IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41695


Ignore:
Timestamp:
Jul 7, 2021, 4:54:07 PM (5 years ago)
Author:
eugene
Message:

merge changes from trunk: relphot improvements, ASCII FITS table I/O fixes, mpcorb_predict

Location:
tags/ipp-ps1-20210510/Ohana
Files:
24 edited
3 copied

Legend:

Unmodified
Added
Removed

  • tags/ipp-ps1-20210510/Ohana

  • tags/ipp-ps1-20210510/Ohana/src/getstar

  • tags/ipp-ps1-20210510/Ohana/src/kapa2/src/Center.c

    r25757 r41695  
    1717
    1818  // enforce integer center here?
    19   image[0].picture.Xc = X;
    20   image[0].picture.Yc = Y;
     19  // NAN value means retain existing center
     20  if (isfinite(X)) image[0].picture.Xc = X;
     21  if (isfinite(Y)) image[0].picture.Yc = Y;
    2122  if ((zoom != 0) && (zoom != -1)) {
    2223    image[0].picture.expand = zoom;

  • tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_get_column.c

    r41474 r41695  
    310310   C - complex float
    311311   M - complex double
    312    P - var lenght array descrpt (64 bit)
     312   P - var length array descrpt (64 bit)
    313313   
    314314   all can be preceeded by integer N which specified number
    … …  
    341341
    342342  if (!gfits_table_format (format, type, Nval, &Nbytes)) return (FALSE);
    343   *Nval = 1;
    344343
    345344  return (TRUE);
    … …  
    379378  gfits_scan (header, "NAXIS2",  OFF_T_FMT, 1,  &Ny);
    380379
    381   /* scan columns to find insert point */
    382   Nstart = 0;
    383   for (i = 1; i < N; i++) {
    384     snprintf (field, 256, "TFORM%d", i);
    385     gfits_scan (header, field, "%s", 1, format);
    386     gfits_table_format (format, tmp, &Nv, &Nb);
    387     Nstart += Nv*Nb;
     380  // find the starting byte for this column
     381  // FITS ASCII table is supposed to have TBCOLn to specify the starting point
     382  // but if it is missing, we can try to find by counting
     383  snprintf (field, 256, "TBCOL%d", N);
     384  int status = gfits_scan (header, field, "%d", 1, &Nstart);
     385  if (!status) {
     386    /* scan columns to find insert point */
     387    Nstart = 0;
     388    for (i = 1; i < N; i++) {
     389      snprintf (field, 256, "TFORM%d", i);
     390      gfits_scan (header, field, "%s", 1, format);
     391      gfits_table_format (format, tmp, &Nv, &Nb);
     392      Nstart += Nv*Nb;
     393    }
     394  } else {
     395    Nstart --;
    388396  }
    389397
    390398  /* allocate temporary line, init pointers */
    391   ALLOCATE (line, char, Nval+1);
    392   bzero (line, Nval+1);
     399  ALLOCATE (line, char, Nval*Nbytes+1);
     400  bzero (line, Nval*Nbytes+1);
    393401  Pin  = table[0].buffer + Nstart;
    394402

  • tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_set_column.c

    r41422 r41695  
    478478  /* print line with appropriate formatting */
    479479  ALLOCATE (array, char, Nbytes*Nval*Nrow);
    480   ALLOCATE (line, char, Nval+1);
     480  ALLOCATE (line, char, Nbytes+1);
    481481  Pin = data;
    482482  Pout = array;
    483483  if (!strcmp (type, "char")) {
    484     for (i = 0; i < Nval*Nrow; i++, Pin++, Pout++) {
     484    for (i = 0; i < Nbytes*Nval*Nrow; i++, Pin++, Pout++) {
    485485      *Pout = *Pin;
    486486    }
    487487  }
    488488  if (!strcmp (type, "int")) {
    489     for (i = 0; i < Nrow; i++, Pin+=4, Pout+=Nval) {
    490       snprintf (line, Nval + 1, cformat, *(int *)Pin);
    491       memcpy (Pout, line, Nval);
     489    for (i = 0; i < Nrow; i++, Pin+=4, Pout+=Nbytes) {
     490      snprintf (line, Nbytes + 1, cformat, *(int *)Pin);
     491      memcpy (Pout, line, Nbytes);
    492492    }
    493493  }
    494494  if (!strcmp (type, "float")) {
    495     for (i = 0; i < Nrow; i++, Pin+=4, Pout+=Nval) {
    496       snprintf (line, Nval + 1, cformat, *(float *)Pin);
    497       memcpy (Pout, line, Nval);
     495    for (i = 0; i < Nrow; i++, Pin+=4, Pout+=Nbytes) {
     496      snprintf (line, Nbytes + 1, cformat, *(float *)Pin);
     497      memcpy (Pout, line, Nbytes);
    498498    }
    499499  }

  • tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_table_format.c

    r41474 r41695  
    55// get the format of a table column
    66// Nval : number of joined columns
    7 // Nbytes : width of column
     7// Nbytes : width of column, determined by the type (e.g., float = 4 bytes)
    88int gfits_bintable_format (char *format, char *type, int *Nval, int *Nbytes) {
    99
    … …  
    9090
    9191/***********************/
     92// get the format of a table column
     93// Nval : number of fields (1 for all types except string)
     94// Nbytes : width of FITS table column, ie number of bytes in ASCII table
    9295int gfits_table_format (char *format, char *type, int *Nval, int *Nbytes) {
    9396
    … …  
    109112  int isLong = FALSE;
    110113  char Type = 'x';
    111   if (Fchar == 'D') { *Nbytes = 1;  strcpy (type, "double"); Type = 'e'; *Nval = Nv; isLong = TRUE; }
    112   if (Fchar == 'E') { *Nbytes = 1;  strcpy (type, "float");  Type = 'e'; *Nval = Nv; }
    113   if (Fchar == 'F') { *Nbytes = 1;  strcpy (type, "float");  Type = 'f'; *Nval = Nv; }
    114   if (Fchar == 'I') { *Nbytes = 1;  strcpy (type, "int");    Type = 'd'; *Nval = Nv; }
    115   if (Fchar == 'A') { *Nbytes = 1;  strcpy (type, "char");   Type = 's'; *Nval = Nv; }
     114  if (Fchar == 'D') { *Nbytes = Nv;  strcpy (type, "double"); Type = 'e'; *Nval =  1; isLong = TRUE; }
     115  if (Fchar == 'E') { *Nbytes = Nv;  strcpy (type, "float");  Type = 'e'; *Nval =  1; }
     116  if (Fchar == 'F') { *Nbytes = Nv;  strcpy (type, "float");  Type = 'f'; *Nval =  1; }
     117  if (Fchar == 'I') { *Nbytes = Nv;  strcpy (type, "int");    Type = 'd'; *Nval =  1; }
     118  if (Fchar == 'A') { *Nbytes =  1;  strcpy (type, "char");   Type = 's'; *Nval = Nv; }
    116119  if (!*Nbytes) { return (FALSE); }
    117120 
    … …  
    133136/*
    134137  valid TABLE column formats:
    135   FN.N  - floating point
    136   INN   - integer
     138  FN.N  - floating point taking NN characters in the table
     139  INN   - integer taking NN characters in the table
    137140  ANN   - NN char string
    138141 

  • tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/Makefile

    r40376 r41695  
    4242$(SRC)/fixcols.$(ARCH).o           \
    4343$(SRC)/fiximage.$(ARCH).o          \
     44$(SRC)/forcedphot.$(ARCH).o        \
    4445$(SRC)/gauss.$(ARCH).o             \
    4546$(SRC)/getvel.$(ARCH).o            \
    … …  
    9091$(SRC)/imfit-qgauss-psf.$(ARCH).o          \
    9192$(SRC)/imfit-sgauss.$(ARCH).o      \
    92 $(SRC)/imfit-sgauss-psf.$(ARCH).o          \
     93$(SRC)/imfit-sgauss-psf.$(ARCH).o  \
    9394$(SRC)/imfit-qfgauss.$(ARCH).o     \
    94 $(SRC)/imfit-qrgauss.$(ARCH).o     
     95$(SRC)/imfit-qrgauss.$(ARCH).o     \
     96$(SRC)/imfit-trail.$(ARCH).o
    9597
    9698# dependancy rules for include files ########################

  • tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/imfit.c

    r39233 r41695  
    6565    qfgauss_setup (argv[N]);
    6666    qrgauss_setup (argv[N]);
     67    trail_setup (argv[N]);
    6768    if (fitfunc == NULL) {
    6869      gprint (GP_ERR, "unknown function %s\n", argv[N]);
    … …  
    124125  dchisq = ochisq;
    125126  chisq  = ochisq;
    126   for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
     127
     128//for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
     129
     130  for (i = 0; (i < 25); i++) {
    127131    chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
    128132    dchisq = ochisq - chisq;
    129133    ochisq = chisq;
     134    fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq);
    130135  } 
    131136  set_int_variable ("Niter",  i);

  • tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/init.c

    r41379 r41695  
    2727int fixcols                 PROTO((int, char **));
    2828int fixrows                 PROTO((int, char **));
     29int forcedphot              PROTO((int, char **));
    2930int gauss                   PROTO((int, char **));
    3031int gaussfit                PROTO((int, char **));
    … …  
    100101  {1, "fixcols",     fixcols,      "fix bad columns by comparing with others"},
    101102  {1, "fixrows",     fixrows,      "fix bad rows by comparing with others"},
     103  {1, "forcedphot",  forcedphot,   "forced photometry on a star or set of stars, assuming gaussian profile"},
    102104  {1, "gauss",       gauss,        "get statistics on a star, assuming gaussian profile"},
    103105  {1, "getvel",      getvel,       "rotcurve to velocities"},

  • tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.data/read_vectors.c

    r41432 r41695  
    706706    char name[80];
    707707     
     708    // for both BINARY and ASCII tables, Nval is the number of fields of the given type
     709    // e.g., there may be 3 related floating point fields fields or 21 connected ascii characters.
     710
    708711    Nval = 0;
    709712    if (Binary) {

  • tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.data/wd.c

    r38441 r41695  
    11# include "data.h"
     2FILE *gfits_open_and_extend (char *filename);
    23
    34int wd (int argc, char **argv) {
    … …  
    9899  // it updates NEXTEND and set EXTEND to TRUE, and modifies the PHU header (neither should happen)
    99100  // if those keywords do not exist...
    100   if (Extend) {
    101     Header Xhead;
    102     FILE *f;
    103     off_t nbytes;
    104     int status, Nextend;
    105 
    106     /* assume failure means non-existent file */
    107     if (!gfits_read_header (argv[2], &Xhead)) {
    108 
    109       gfits_init_header (&Xhead);
    110       Xhead.bitpix = 16;
    111       Xhead.extend = TRUE;
    112       gfits_create_header (&Xhead);
    113 
    114       gfits_modify (&Xhead, "NEXTEND", "%d", 1, 0);
    115       f = fopen (argv[2], "w");
    116       fclose (f);
    117     }
    118 
    119     gfits_modify_alt (&Xhead, "EXTEND", "%t", 1, TRUE);
    120 
    121     Nextend = 0;
    122     gfits_scan (&Xhead, "NEXTEND", "%d", 1, &Nextend);
    123     Nextend ++;
    124     gfits_modify (&Xhead, "NEXTEND", "%d", 1, Nextend);
    125 
    126     /* write the main header to the start of the file */
    127     f = fopen (argv[2], "r+");
    128     if (f == NULL) {
    129       gprint (GP_ERR, "failed to write file\n");
     101  if (Extend && !CompressMode) {
     102    int status = TRUE;
     103
     104    FILE *f = gfits_open_and_extend (argv[2]);
     105    if (!f) {
     106      gprint (GP_ERR, "failed to extend file %s\n", argv[2]);
    130107      status = FALSE;
    131108      goto done1;
    132109    }
    133    
    134     /* position to begining of file to write header */
    135     fseeko (f, 0LL, SEEK_SET);
    136     nbytes = fwrite (Xhead.buffer, 1, Xhead.datasize, f);
    137     if (nbytes != Xhead.datasize) {
    138       gprint (GP_ERR, "ERROR: failed writing data to image header\n");
    139       status = FALSE;
    140       goto done1;
    141     }
    142    
     110
    143111    /* fix up header */
    144112    {
    … …  
    153121    /* position to end of file to write new extend */
    154122    fseeko (f, 0LL, SEEK_END);
    155     nbytes = fwrite (temp_header.buffer, 1, temp_header.datasize, f);
    156     fclose (f);
     123    off_t nbytes = fwrite (temp_header.buffer, 1, temp_header.datasize, f);
    157124    if (nbytes != temp_header.datasize) {
     125      fclose (f);
    158126      gprint (GP_ERR, "failed to write file\n");
    159127      status = FALSE;
    160128      goto done1;
    161129    }
     130
    162131    /* write the matrix buffer (automatically goes to end of file */
    163     if (!gfits_write_matrix (argv[2], &temp_matrix)) {
     132    if (!gfits_fwrite_matrix (f, &temp_matrix)) {
     133      fclose (f);
    164134      gprint (GP_ERR, "failed to write file\n");
    165135      status = FALSE;
    166136      goto done1;
    167137    }
     138    fclose (f);
    168139    status = TRUE;
     140
    169141  done1:
    170     gfits_free_header (&Xhead);
    171142    gfits_free_header (&temp_header);
    172143    gfits_free_matrix (&temp_matrix);
    … …  
    180151    Matrix myMatrix;
    181152
    182     FILE *f = fopen (argv[2], "w");
    183     if (!f) {
    184       fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", argv[2]);
    185       FREE (CompressMode);
    186       return FALSE;
    187     }
    188    
    189     gfits_init_header (&myHeader);
    190     myHeader.extend = TRUE;
    191     gfits_create_header (&myHeader);
    192     gfits_create_matrix (&myHeader, &myMatrix);
    193     gfits_fwrite_header  (f, &myHeader);
    194     gfits_fwrite_matrix  (f, &myMatrix);
    195     gfits_free_header (&myHeader);
    196     gfits_free_matrix (&myMatrix);
     153    FILE *f;
     154
     155    if (Extend) {
     156      // keeps an existing file
     157      f = gfits_open_and_extend (argv[2]);
     158      fseeko (f, 0LL, SEEK_END);
     159    } else {
     160      // replaces an existing file
     161      f = fopen (argv[2], "w");
     162      if (!f) {
     163        fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", argv[2]);
     164        FREE (CompressMode);
     165        return FALSE;
     166      }
     167   
     168      gfits_init_header (&myHeader);
     169      myHeader.extend = TRUE;
     170      gfits_create_header (&myHeader);
     171      gfits_create_matrix (&myHeader, &myMatrix);
     172      gfits_fwrite_header  (f, &myHeader);
     173      gfits_fwrite_matrix  (f, &myMatrix);
     174      gfits_free_header (&myHeader);
     175      gfits_free_matrix (&myMatrix);
     176    }
    197177
    198178    FTable ftable;
    … …  
    236216  return (TRUE);
    237217}
     218
     219// prepare the file to add an extension:
     220// read the PHU, update NEXTEND, write back header
     221// return the file pointer open and ready for write
     222FILE *gfits_open_and_extend (char *filename) {
     223
     224  Header Xhead;
     225  FILE *f = NULL;
     226  struct stat filestat;
     227
     228  // does the file already exist?
     229  int status = stat (filename, &filestat);
     230  if (status == 0) { /* file exists, are permissions OK? */
     231    f = fopen (filename, "r+");
     232    if (f == NULL) {
     233      gprint (GP_ERR, "failed to open file %s for write\n", filename);
     234      return NULL;
     235    }
     236
     237    status = gfits_fread_header (f, &Xhead);
     238    if (!status) {
     239      gprint (GP_ERR, "failed to read header for existing file %s\n", filename);
     240      return NULL;
     241    }
     242  } else {
     243    f = fopen (filename, "w");
     244    if (f == NULL) {
     245      gprint (GP_ERR, "failed to create file %s for write\n", filename);
     246      return NULL;
     247    }
     248
     249    // create an empty header
     250    gfits_init_header (&Xhead);
     251    Xhead.bitpix = 16;
     252    Xhead.extend = TRUE;
     253    gfits_create_header (&Xhead);
     254
     255    gfits_modify (&Xhead, "NEXTEND", "%d", 1, 0);
     256  }
     257
     258  gfits_modify_alt (&Xhead, "EXTEND", "%t", 1, TRUE);
     259
     260  int Nextend = 0;
     261  gfits_scan (&Xhead, "NEXTEND", "%d", 1, &Nextend);
     262  Nextend ++;
     263  gfits_modify (&Xhead, "NEXTEND", "%d", 1, Nextend);
     264
     265  /* position to begining of file to write header */
     266  fseeko (f, 0LL, SEEK_SET);
     267  off_t nbytes = fwrite (Xhead.buffer, 1, Xhead.datasize, f);
     268  if (nbytes != Xhead.datasize) {
     269    gfits_free_header (&Xhead);
     270    gprint (GP_ERR, "ERROR: failed writing data to image header\n");
     271    return NULL;
     272  }
     273  gfits_free_header (&Xhead);
     274   
     275  return (f);
     276}
     277
     278
     279///    /* fix up header */
     280///    {
     281///      static char simple[] = "XTENSION= 'IMAGE  '            / Image extension";
     282///      int Ns, No;
     283///      Ns = strlen (simple);
     284///      No = 80 - Ns;
     285///      strncpy (temp_header.buffer, simple, Ns);
     286///      memset (&temp_header.buffer[Ns], ' ', No);
     287///    }
     288///
     289///    /* position to end of file to write new extend */
     290///    fseeko (f, 0LL, SEEK_END);
     291///    nbytes = fwrite (temp_header.buffer, 1, temp_header.datasize, f);
     292///    fclose (f);
     293///    if (nbytes != temp_header.datasize) {
     294///      gprint (GP_ERR, "failed to write file\n");
     295///      status = FALSE;
     296///      goto done1;
     297///    }
     298///
     299///  }

  • tags/ipp-ps1-20210510/Ohana/src/opihi/include/imfit.h

    r21153 r41695  
    1919void qgauss_psf_setup (char *name);
    2020void sgauss_psf_setup (char *name);
     21
     22void trail_setup (char *name);

  • tags/ipp-ps1-20210510/Ohana/src/relphot

  • tags/ipp-ps1-20210510/Ohana/src/relphot/include/relphot.h

    r41649 r41695  
    163163
    164164typedef enum {
    165   STATS_NONE,
    166   STATS_MEAN,
    167   STATS_MEDIAN,
    168   STATS_WT_MEAN,
    169   STATS_INNER_MEAN,
    170   STATS_INNER_WTMEAN,
    171   STATS_CHI_INNER_MEAN,
    172   STATS_CHI_INNER_WTMEAN
     165  // these modes are primary and mutually exclusive:
     166  STATS_NONE               = 0x0000,
     167  STATS_MEAN               = 0x0001,
     168  STATS_MEDIAN             = 0x0002,
     169  STATS_WT_MEAN            = 0x0003,
     170  STATS_INNER_MEAN         = 0x0004,
     171  STATS_INNER_WTMEAN       = 0x0005,
     172  STATS_CHI_INNER_MEAN     = 0x0006,
     173  STATS_CHI_INNER_WTMEAN   = 0x0007,
     174  // these modes are additional options
     175  STATS_PRIMARY            = 0x0007, // use this to mask out the optional bits
     176  STATS_VARSTATS           = 0x0010,
    173177} ListStatsMode;
    174178
    … …  
    185189  double Lower20;
    186190  double Lower10;
     191  double Upper90Nsig;
     192  double Upper80Nsig;
     193  double Lower20Nsig;
     194  double Lower10Nsig;
    187195  double total;
    188196  int    Nmeas;
    … …  
    419427double IMFIT_SYS_SIGMA_LIM;
    420428double CLOUD_TOLERANCE;
     429int    VARIABILITY_STATS;
     430int    USE_OLS_FOR_AVERAGES;
    421431
    422432int    PARALLEL;
    … …  
    874884void ResetAverageActivePhotcodes (SecFilt *secfilt);
    875885
    876 
     886void   rationalize_zeropoints (int Niter);
     887double get_median_zpt_images (short photcode);
     888void   set_median_zpt_images (short photcode, double zpt);
     889double get_median_zpt_mosaics (short photcode);
     890void   set_median_zpt_mosaics (short photcode, double zpt);
     891double get_median_zpt_tgroups (short photcode);
     892void   set_median_zpt_tgroups (short photcode, double zpt);

  • tags/ipp-ps1-20210510/Ohana/src/relphot/src/GridOps.c

    r41649 r41695  
    274274    if (!GridCorr[code]) continue;
    275275
     276    // float GridSum = 0.0;
     277    // int   GridCnt =   0;
    276278    for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
    277279      for (int iy = 0; iy < GridCorr[code]->Ny; iy++) {
    … …  
    292294        GridCorr[code]->dMgrid[ix][iy] = sqrt(r*(Mgrid2 - Mgrid*Mgrid)); // sample stdev
    293295        // fprintf (stderr, "grid code %d, %d x %d : %f +/- %f : %d\n", code, ix, iy, GridCorr[code]-> Mgrid[ix][iy], GridCorr[code]->dMgrid[ix][iy], GridCorr[code]->nMgrid[ix][iy]);
     296        // GridSum += Mgrid;
     297        // GridCnt ++;
    294298      }
    295299    }
     300    // float GridAve = GridSum / GridCnt;
     301    // fprintf (stderr, "grid average: %f\n", GridAve);
    296302  }
    297303  return;

  • tags/ipp-ps1-20210510/Ohana/src/relphot/src/ImageOps.c

    r41649 r41695  
    951951}
    952952
     953// find the median zero point and subtract it (for each active average photcode)
     954void rationalize_zeropoints (int Niter) {
     955
     956  if (VERBOSE) fprintf (stderr, "rationalize zero points\n");
     957
     958  // if we are calculating zero points for tgroups,
     959  // rationalize based on the tgroups
     960  if (TGROUP_ZEROPT) {
     961    for (int ic = 0; ic < Nphotcodes; ic++) {
     962      double zpt = get_median_zpt_tgroups (photcodes[ic][0].code);
     963      fprintf (stderr, "rationalize zero points by tgroup for %s (%d) : %f\n", photcodes[ic][0].name, photcodes[ic][0].code, zpt);
     964      if (isnan(zpt)) continue;
     965      set_median_zpt_tgroups (photcodes[ic][0].code, zpt);
     966      set_median_zpt_mosaics (photcodes[ic][0].code, zpt);
     967      set_median_zpt_images  (photcodes[ic][0].code, zpt);
     968    }
     969    return;
     970  }
     971
     972  // if we are calculating zero points for mosaics, but not tgroups,
     973  // rationalize based on the mosaics
     974  if (MOSAIC_ZEROPT) {
     975    for (int ic = 0; ic < Nphotcodes; ic++) {
     976      double zpt = get_median_zpt_mosaics (photcodes[ic][0].code);
     977      fprintf (stderr, "rationalize zero points by mosaic for %s (%d) : %f\n", photcodes[ic][0].name, photcodes[ic][0].code, zpt);
     978      if (isnan(zpt)) continue;
     979      set_median_zpt_tgroups (photcodes[ic][0].code, zpt);
     980      set_median_zpt_mosaics (photcodes[ic][0].code, zpt);
     981      set_median_zpt_images  (photcodes[ic][0].code, zpt);
     982    }
     983    return;
     984  }   
     985
     986  // otherwise, rationalize based on the images
     987  for (int ic = 0; ic < Nphotcodes; ic++) {
     988    double zpt = get_median_zpt_images (photcodes[ic][0].code);
     989    fprintf (stderr, "rationalize zero points by image for %s (%d) : %f\n", photcodes[ic][0].name, photcodes[ic][0].code, zpt);
     990    if (isnan(zpt)) continue;
     991    set_median_zpt_tgroups (photcodes[ic][0].code, zpt);
     992    set_median_zpt_mosaics (photcodes[ic][0].code, zpt);
     993    set_median_zpt_images  (photcodes[ic][0].code, zpt);
     994  }
     995}
     996
     997double get_median_zpt_images (short photcode) {
     998
     999  double *mlist;
     1000
     1001  ALLOCATE (mlist, double, Nimage);
     1002
     1003  int N = 0;
     1004  for (int i = 0; i < Nimage; i++) {
     1005    int mycode = GetPhotcodeEquivCodebyCode (image[i].photcode);
     1006    if (mycode != photcode) continue;
     1007    if (!(image[i].flags & ID_IMAGE_IMAGE_PHOTCAL)) continue;
     1008    mlist[N] = image[i].McalPSF;
     1009    N++;
     1010  }
     1011
     1012  if (N == 0) {
     1013    free (mlist);
     1014    return NAN;
     1015  }
     1016
     1017  StatType stats;
     1018  liststats_setmode (&stats, "MEAN");
     1019
     1020  liststats (mlist, NULL, NULL, N, &stats);
     1021  free (mlist);
     1022
     1023  fprintf (stderr, "rationalize by image using %d pts\n", N);
     1024  return stats.median;
     1025}
     1026
     1027void set_median_zpt_images (short photcode, double zpt) {
     1028
     1029  if (!isfinite(zpt)) return; // do not break the zero points
     1030
     1031  for (int i = 0; i < Nimage; i++) {
     1032    int mycode = GetPhotcodeEquivCodebyCode (image[i].photcode);
     1033    if (mycode != photcode) continue;
     1034    // fprintf (stderr, "IMAGE %d zpt %f -> ", i, image[i].McalPSF);
     1035
     1036    int applyOffset = TRUE;
     1037    TGroup *mygrp = getTGroupForImage (i);
     1038    if (mygrp && (mygrp->flags & ID_IMAGE_TGROUP_PHOTCAL)) applyOffset = FALSE;
     1039
     1040    Mosaic *mymos = getMosaicForImage (i);
     1041    if (mymos && (mymos->flags & ID_IMAGE_MOSAIC_PHOTCAL)) applyOffset = FALSE;
     1042
     1043    if (applyOffset) {
     1044      image[i].McalPSF -= zpt;
     1045      image[i].McalAPER -= zpt;
     1046    }
     1047    // fprintf (stderr, "%f (%d)\n", image[i].McalPSF, applyOffset);
     1048  }
     1049
     1050  return;
     1051}
     1052
    9531053static int setMcal_init_done = FALSE;
    9541054void plot_setMcal (double *list, int Npts) {

  • tags/ipp-ps1-20210510/Ohana/src/relphot/src/MosaicOps.c

    r41649 r41695  
    15741574}
    15751575
     1576TGroup *getTGroupForMosaic (int mos) {
     1577  if (mos >= Nmosaic) return NULL;
     1578  if (mos < 0) return NULL;
     1579
     1580  if (!mosaic[mos].inTGroup) return NULL;
     1581
     1582  int imageIdx = MosaicToImage[mos][0];
     1583
     1584  TGroup *mygrp = getTGroupForImage (imageIdx);
     1585  return mygrp;
     1586}
     1587
     1588double get_median_zpt_mosaics (short photcode) {
     1589
     1590  double *mlist;
     1591
     1592  if (!MOSAIC_ZEROPT) return NAN;
     1593
     1594  ALLOCATE (mlist, double, Nmosaic);
     1595
     1596  int N = 0;
     1597  for (int i = 0; i < Nmosaic; i++) {
     1598    if (mosaic[i].photcode != photcode) continue;
     1599    if (!(mosaic[i].flags & ID_IMAGE_MOSAIC_PHOTCAL)) continue;
     1600    mlist[N] = mosaic[i].McalPSF;
     1601    N++;
     1602  }
     1603
     1604  if (N == 0) {
     1605    free (mlist);
     1606    return NAN;
     1607  }
     1608
     1609  StatType stats;
     1610  liststats_setmode (&stats, "MEAN");
     1611
     1612  liststats (mlist, NULL, NULL, N, &stats);
     1613  free (mlist);
     1614
     1615  fprintf (stderr, "rationalize by mosaic using %d pts\n", N);
     1616  return stats.median;
     1617}
     1618
     1619void set_median_zpt_mosaics (short photcode, double zpt) {
     1620
     1621  if (!MOSAIC_ZEROPT) return;
     1622  if (!isfinite(zpt)) return; // do not break the zero points
     1623
     1624  for (int i = 0; i < Nmosaic; i++) {
     1625    if (mosaic[i].photcode != photcode) continue;
     1626    // fprintf (stderr, "MOSAIC %d zpt %f -> ", i, mosaic[i].McalPSF);
     1627
     1628    int applyOffset = TRUE;
     1629    TGroup *mygrp = getTGroupForMosaic (i);
     1630    if (mygrp && (mygrp->flags & ID_IMAGE_TGROUP_PHOTCAL)) applyOffset = FALSE;
     1631    if (!(mosaic[i].flags & ID_IMAGE_MOSAIC_PHOTCAL)) applyOffset = FALSE;
     1632
     1633    if (applyOffset) {
     1634      mosaic[i].McalPSF -= zpt;
     1635      mosaic[i].McalAPER -= zpt;
     1636    }
     1637    // fprintf (stderr, "%f (%d)\n", mosaic[i].McalPSF, applyOffset);
     1638  }
     1639
     1640  return;
     1641}
     1642
    15761643void plot_mosaic_fields (Catalog *catalog) {
    15771644

  • tags/ipp-ps1-20210510/Ohana/src/relphot/src/TGroupOps.c

    r41649 r41695  
    13351335      }
    13361336      // too few exposures (configure)
    1337       if (tgroup[j].Nmosaic < 4) {
     1337      if ((tgroup[j].Nmosaic < 4) && (tgroup[j].Nimage < 4)) {
    13381338        mark = TRUE;
    13391339        NfewExp ++;
    … …  
    13651365}
    13661366
     1367double get_median_zpt_tgroups (short photcode) {
     1368
     1369  double *mlist;
     1370
     1371  if (!TGROUP_ZEROPT) return NAN;
     1372
     1373  ALLOCATE (mlist, double, NtgroupTimes);
     1374
     1375  int N = 0;
     1376  for (int i = 0; i < NtgroupTimes; i++) {
     1377    TGroup *tgroup = tgroupTimes[i][0].byCode;
     1378    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
     1379      if (tgroup[j].photcode != photcode) continue;
     1380      if (!(tgroup[j].flags & ID_IMAGE_TGROUP_PHOTCAL)) continue;
     1381      mlist[N] = tgroup[j].McalPSF;
     1382      N++;
     1383    }
     1384  }
     1385
     1386  if (N == 0) {
     1387    free (mlist);
     1388    return NAN;
     1389  }
     1390
     1391  StatType stats;
     1392  liststats_setmode (&stats, "MEAN");
     1393
     1394  liststats (mlist, NULL, NULL, N, &stats);
     1395  free (mlist);
     1396
     1397  fprintf (stderr, "rationalize by tgroup using %d pts\n", N);
     1398  return stats.median;
     1399}
     1400
     1401void set_median_zpt_tgroups (short photcode, double zpt) {
     1402
     1403  if (!TGROUP_ZEROPT) return;
     1404  if (!isfinite(zpt)) return; // do not break the zero points
     1405
     1406  for (int i = 0; i < NtgroupTimes; i++) {
     1407    TGroup *tgroup = tgroupTimes[i][0].byCode;
     1408    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
     1409      if (tgroup[j].photcode != photcode) continue;
     1410      if (!(tgroup[j].flags & ID_IMAGE_TGROUP_PHOTCAL)) continue;
     1411      tgroup[j].McalPSF -= zpt;
     1412      tgroup[j].McalAPER -= zpt;
     1413    }
     1414  }
     1415  return;
     1416}
     1417
    13671418void plot_tgroup_fields (Catalog *catalog) {
    13681419

  • tags/ipp-ps1-20210510/Ohana/src/relphot/src/args.c

    r41649 r41695  
    434434    remove_argument (N, &argc, argv);
    435435    KEEP_UBERCAL = FALSE;
     436  }
     437
     438  VARIABILITY_STATS = FALSE;
     439  if ((N = get_argument (argc, argv, "-varstats"))) {
     440    remove_argument (N, &argc, argv);
     441    VARIABILITY_STATS = TRUE;
     442  }
     443  USE_OLS_FOR_AVERAGES = FALSE;
     444  if ((N = get_argument (argc, argv, "-use-ols-for-averages"))) {
     445    remove_argument (N, &argc, argv);
     446    USE_OLS_FOR_AVERAGES = TRUE;
    436447  }
    437448
    … …  
    865876
    866877
     878  VARIABILITY_STATS = FALSE;
     879  if ((N = get_argument (argc, argv, "-varstats"))) {
     880    remove_argument (N, &argc, argv);
     881    VARIABILITY_STATS = TRUE;
     882  }
     883  USE_OLS_FOR_AVERAGES = FALSE;
     884  if ((N = get_argument (argc, argv, "-use-ols-for-averages"))) {
     885    remove_argument (N, &argc, argv);
     886    USE_OLS_FOR_AVERAGES = TRUE;
     887  }
     888
    867889  MIN_ERROR = 0.001;
    868890  if ((N = get_argument (argc, argv, "-minerror"))) {

  • tags/ipp-ps1-20210510/Ohana/src/relphot/src/liststats.c

    r41555 r41695  
    33void liststats_setmode (StatType *stats, char *strmode) {
    44
    5   stats->statmode = -1;
     5  stats->statmode = STATS_NONE;
    66  if (!strcmp (strmode, "MEAN"))             { stats->statmode = STATS_MEAN; return; }
    77  if (!strcmp (strmode, "MEDIAN"))           { stats->statmode = STATS_MEDIAN; return; }
    … …  
    1717
    1818int liststats_init (StatType *stats) {
    19   stats->median  = NAN;
    20   stats->mean    = NAN;
    21   stats->sigma   = NAN;
    22   stats->error   = NAN;
    23   stats->chisq   = NAN;
    24   stats->min     = NAN;
    25   stats->max     = NAN;
    26   stats->Upper80 = NAN;
    27   stats->Lower20 = NAN;
    28   stats->Upper90 = NAN;
    29   stats->Lower10 = NAN;
    30   stats->total   = NAN;
    31   stats->Nmeas = 0;
     19  stats->median      = NAN;
     20  stats->mean        = NAN;
     21  stats->sigma       = NAN;
     22  stats->error       = NAN;
     23  stats->chisq       = NAN;
     24  stats->min         = NAN;
     25  stats->max         = NAN;
     26  stats->Upper80     = NAN;
     27  stats->Lower20     = NAN;
     28  stats->Upper90     = NAN;
     29  stats->Lower10     = NAN;
     30  stats->Upper80Nsig = NAN;
     31  stats->Lower20Nsig = NAN;
     32  stats->Upper90Nsig = NAN;
     33  stats->Lower10Nsig = NAN;
     34  stats->total       = NAN;
     35  stats->Nmeas       = 0;
    3236  return TRUE;
    3337}
    … …  
    3943
    4044  myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set");
     45  ListStatsMode myMode = stats->statmode & STATS_PRIMARY; // exclude the option bits
    4146
    4247  liststats_init (stats);
    … …  
    7580  stats[0].Lower10 = value[N10];
    7681
    77   switch (stats->statmode) {
     82  switch (myMode) {
    7883    case STATS_MEDIAN:
    7984      ks = 0;
    … …  
    100105      break;
    101106    case STATS_NONE:
     107    default:
    102108      myAbort ("undefined stats");
    103109  }   
    … …  
    105111  // for these two modes, I need a vector of the chi-square contribution
    106112  // I'm actually just using chisq to get the correct sorting order
    107   if ((stats->statmode == STATS_CHI_INNER_MEAN) || (stats->statmode == STATS_CHI_INNER_WTMEAN)) {
     113  if ((myMode == STATS_CHI_INNER_MEAN) || (myMode == STATS_CHI_INNER_WTMEAN)) {
    108114    ALLOCATE (chi, double, N);
    109115    for (i = 0; i < N; i++) {
    … …  
    119125
    120126  int WeightedMean = FALSE;
    121   WeightedMean |= (stats->statmode == STATS_WT_MEAN);
    122   WeightedMean |= (stats->statmode == STATS_INNER_WTMEAN);
    123   WeightedMean |= (stats->statmode == STATS_CHI_INNER_WTMEAN);
     127  WeightedMean |= (myMode == STATS_WT_MEAN);
     128  WeightedMean |= (myMode == STATS_INNER_WTMEAN);
     129  WeightedMean |= (myMode == STATS_CHI_INNER_WTMEAN);
    124130  if (!dvalue) WeightedMean = FALSE;  // warn the user?
    125131
    … …  
    278284
    279285  // make the initial guess based on the median (not weighted mean)
    280   // ALLOCATE_PTR (values, double, Npoints);
    281286  for (int i = 0; i < Npoints; i++) {
    282287    dataset->values[i] = dataset->flxlist[i];
    … …  
    288293  // if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, NULL, Npoints)) return FALSE;
    289294 
    290   // XXX add to dataset elements?
    291   // ALLOCATE_PTR (wtvals, double, Npoints);
    292   // ALLOCATE_PTR (wtlist, double, Npoints);
    293 
    294295  int converged = FALSE;
    295296  for (int iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations++) {
    … …  
    315316     
    316317  // calculate the weight thresholds to mask the bad points:
    317   // double Sum_W = 0.0;
    318318  for (int i = 0; i < Npoints; i++) {
    319319    dataset->wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
    320320    dataset->wtlist[i] = dataset->wtvals[i];
    321     // Sum_W += dataset->wtvals[i];
    322321  }
    323322  dsort (dataset->wtlist, Npoints);
    324323  double WtMedian = (Npoints % 2) ? dataset->wtlist[midpt] : 0.5*(dataset->wtlist[midpt] + dataset->wtlist[midpt-1]);
    325   // double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * Npoints);
    326324  double WtThreshold = WEIGHT_THRESHOLD * WtMedian;
    327 
    328   // generate the unmasked subset
    329   // XXX add these to the dataset elements?
    330   // ALLOCATE_PTR ( ykeep, double, Npoints);
    331   // ALLOCATE_PTR (dykeep, double, Npoints);
    332   // ALLOCATE_PTR (wtkeep, double, Npoints);
    333325
    334326  // save unmasked points
    … …  
    343335      continue;
    344336    }
    345      dataset->ykeep[Nkeep] = dataset->flxlist[i];
     337    dataset-> ykeep[Nkeep] = dataset->flxlist[i];
    346338    dataset->dykeep[Nkeep] = dataset->errlist[i];
    347339    dataset->wtkeep[Nkeep] = dataset->wgtlist[i]; // externally-supplied weight
    … …  
    359351  stats->sigma = sqrt (dSig / (Nkeep - 1));
    360352
    361   // bootstrap resampling to generate the errorbars
    362   // XXX add these to the dataset elements?
    363   // ALLOCATE_PTR (ysample,  double, Nkeep);
    364   // ALLOCATE_PTR (dysample, double, Nkeep);
    365   // ALLOCATE_PTR (wtsample, double, Nkeep);
    366   // ALLOCATE_PTR (bvalue,   double, NBOOTSTRAP); // vector to save the bootstrap values
     353  // if percentile ranges are desired, calculate them
     354  if (stats->statmode & STATS_VARSTATS) {
     355    // save ykeep values in a vector to be sorted
     356    for (int i = 0; i < Nkeep; i++) {
     357      dataset->wtvals[i] = dataset->ykeep[i];
     358    }
     359    dsort (dataset->wtvals, Nkeep);
     360
     361    int N90 = MIN (Nkeep-1, 0.9*Nkeep);
     362    int N80 = MIN (Nkeep-1, 0.8*Nkeep);
     363    int N20 = MAX (0, 0.2*Nkeep);
     364    int N10 = MAX (0, 0.1*Nkeep);
     365
     366    stats->Upper90 = dataset->wtvals[N90];
     367    stats->Upper80 = dataset->wtvals[N80];
     368    stats->Lower20 = dataset->wtvals[N20];
     369    stats->Lower10 = dataset->wtvals[N10];
     370   
     371    int NkeepMid = 0.5 * Nkeep;
     372    double ClipMedian = (Nkeep % 2) ? dataset->wtlist[NkeepMid] : 0.5*(dataset->wtlist[NkeepMid] + dataset->wtlist[NkeepMid-1]);
     373
     374    // save ykeep values in a vector to be sorted
     375    for (int i = 0; i < Nkeep; i++) {
     376      dataset->wtvals[i] = (dataset->ykeep[i] - ClipMedian) / dataset->dykeep[i];
     377    }
     378    dsort (dataset->wtvals, Nkeep);
     379
     380    stats->Upper90Nsig = dataset->wtvals[N90];
     381    stats->Upper80Nsig = dataset->wtvals[N80];
     382    stats->Lower20Nsig = dataset->wtvals[N20];
     383    stats->Lower10Nsig = dataset->wtvals[N10];
     384  }
    367385
    368386  int Nboot = 0;

  • tags/ipp-ps1-20210510/Ohana/src/relphot/src/reload_catalogs.c

    r41649 r41695  
    237237    if (USE_BASIC_CHECK)   {  strextend (&command, "-basic-image-search"); }
    238238    // if (USE_ALL_IMAGES)    { strextend (&command, "-use-all-images"); }
     239    if (VARIABILITY_STATS) { strextend (&command, "-varstats"); }
    239240    if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); }
    240241

  • tags/ipp-ps1-20210510/Ohana/src/relphot/src/relphot_images.c

    r41649 r41695  
    9696    for (i = 0; i < NLOOP; i++) {
    9797      SetZptIteration (i);
     98
     99      rationalize_zeropoints (i);
    98100
    99101      setMrel  (catalog, Ncatalog); // threaded (calls setMrelCatalog)

  • tags/ipp-ps1-20210510/Ohana/src/relphot/src/setMrelCatalog.c

    r41649 r41695  
    11# include "relphot.h"
     2
     3# define UPPER90_VAL MpsfStk
     4# define UPPER80_VAL FpsfStk
     5# define LOWER20_VAL MkronStk
     6# define LOWER10_VAL FkronStk
     7# define UPPER90_SIG MpsfWrp
     8# define UPPER80_SIG FpsfWrp
     9# define LOWER20_SIG MkronWrp
     10# define LOWER10_SIG FkronWrp
     11
    212
    313# if (0)
    … …  
    7181  SetMrelInfoReset (results); // reset the counters
    7282
    73   useOLS = UseStandardOLS (ZPT_STARS);
     83  // XXX if we are iterating on the image zero points, we should use OLS
     84  // until we have performed a few iterations.  UseStandardOLS() returns
     85  // true (use OLS) or false depending on the zero point mode (TGROUP, MOSAIC, IMAGE) and the
     86  // number of iterations.
     87
     88  // On the final pass (isSetMrelFinal), we should use OLS or IRLS depending on
     89  // the user choice (USE_OLS_FOR_AVERAGES)
     90
     91  useOLS = isSetMrelFinal ? USE_OLS_FOR_AVERAGES : UseStandardOLS (ZPT_STARS);
    7492  for (j = 0; j < catalog[Nc].Naverage; j++) {
    7593
    … …  
    435453    StatType *psfstats = &results->psfstats;
    436454
     455    // if isSetMrelFinal and varstats are requested, request percentile ranges
     456    if (isSetMrelFinal && VARIABILITY_STATS) psfstats->statmode |= STATS_VARSTATS;
     457
    437458    int NrankingPSF = useOLS ? magStatsByRanking (&results->psfData[Nsec], psfstats) : magStatsByRankingIRLS (&results->psfData[Nsec], psfstats);
    438459    if (NrankingPSF < Nminmeas) {
    … …  
    469490        secfilt[Nsec].Mmax     = psfstats->max;
    470491        secfilt[Nsec].Mmin     = psfstats->min;
     492      }
     493      if (VARIABILITY_STATS) {
     494        // XXX NOTE: these are overloaded on warp and stack mags and fluxes
     495        secfilt[Nsec].UPPER90_VAL = psfstats->Upper90;
     496        secfilt[Nsec].UPPER80_VAL = psfstats->Upper80;
     497        secfilt[Nsec].LOWER20_VAL = psfstats->Lower20;
     498        secfilt[Nsec].LOWER10_VAL = psfstats->Lower10;
     499        secfilt[Nsec].UPPER90_SIG = psfstats->Upper90Nsig;
     500        secfilt[Nsec].UPPER80_SIG = psfstats->Upper80Nsig;
     501        secfilt[Nsec].LOWER20_SIG = psfstats->Lower20Nsig;
     502        secfilt[Nsec].LOWER10_SIG = psfstats->Lower10Nsig;
    471503      }
    472504      secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];

  • tags/ipp-ps1-20210510/Ohana/src/tools/Makefile

    r40547 r41695  
    6161$(PROGRAMS) : % : $(BIN)/%.$(ARCH)
    6262
     63$(BIN)/mpcorb_predict.$(ARCH) : $(SRC)/mpcorb_predict.$(ARCH).o
     64        @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi
     65        $(CC) $(FULL_CFLAGS) -o $@ $^ $(FULL_LDFLAGS) -lsla
     66        @echo "compiled $*"
     67        @echo ""
     68
     69$(DESTBIN)/mpcorb_predict: $(BIN)/mpcorb_predict.$(ARCH)
     70        @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi
     71        rm -f $(DESTBIN)/mpcorb_predict
     72        cp $(BIN)/mpcorb_predict.$(ARCH) $(DESTBIN)/mpcorb_predict
     73        @echo "installed mpcorb_predict"
     74        @echo ""
     75
     76$(DESTLIB)/libsla.a:
     77        echo "missing SLALIB libsla.a: install in $(DESTLIB)"
     78
     79$(DESTINC)/slalib.h:
     80        echo "missing SLALIB slalib.h: install in $(DESTINC)"
     81
     82mpcorb_predict: $(DESTBIN)/mpcorb_predict $(DESTLIB)/libsla.a $(DESTINC)/slalib.h
     83
    6384%.clean:
    6485        rm -f $(SRC)/$*.$(ARCH).o
Note: See TracChangeset for help on using the changeset viewer.