IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33180


Ignore:
Timestamp:
Jan 30, 2012, 7:43:27 AM (14 years ago)
Author:
eugene
Message:

merge the NO_METADATA code into load_zpt_ubercal

Location:
branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/load_zpt_table.c

    r33156 r33180  
    5353ZptTable *load_zpt_ubercal(char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable) {
    5454
    55   int i, nfilter, nseason, ix, iy, ixc, iyc, Ncol, seq;
     55  int i, nfilter, nseason, ix, iy, ixc, iyc, Ncol;
    5656  off_t Nrow;
    5757  char type[16], filter[80];
    … …  
    101101  // for current testing, make fake smfs that correspond to specific chips, filter, and mjd ranges?
    102102
    103   if (!gfits_scan (&header, "NSEASON", "%d", 1, &NSEASON)) {
    104     fprintf (stderr, "cannot find NSEASON in header of %s\n", filename);
    105     fclose (f);
    106     return NULL;
    107   }
    108   if (!gfits_scan (&header, "NFILTER", "%d", 1, &NFILTER)) {
    109     fprintf (stderr, "cannot find NFILTER in header of %s\n", filename);
    110     fclose (f);
    111     return NULL;
    112   }
    113   if (!gfits_scan (&header, "NCHIP_X", "%d", 1, &NCHIP_X)) {
    114     fprintf (stderr, "cannot find NCHIP_X in header of %s\n", filename);
    115     fclose (f);
    116     return NULL;
    117   }
    118   if (!gfits_scan (&header, "NCHIP_Y", "%d", 1, &NCHIP_Y)) {
    119     fprintf (stderr, "cannot find NCHIP_Y in header of %s\n", filename);
    120     fclose (f);
    121     return NULL;
    122   }
    123   if (!gfits_scan (&header, "NCELL_X", "%d", 1, &NCELL_X)) {
    124     fprintf (stderr, "cannot find NCELL_X in header of %s\n", filename);
    125     fclose (f);
    126     return NULL;
    127   }
    128   if (!gfits_scan (&header, "NCELL_Y", "%d", 1, &NCELL_Y)) {
    129     fprintf (stderr, "cannot find NCELL_Y in header of %s\n", filename);
    130     fclose (f);
    131     return NULL;
    132   }
    133   if (!gfits_scan (&header, "CHIP_DX", "%d", 1, &CHIP_DX)) {
    134     fprintf (stderr, "cannot find CHIP_DX in header of %s\n", filename);
    135     fclose (f);
    136     return NULL;
    137   }
    138   if (!gfits_scan (&header, "CHIP_DY", "%d", 1, &CHIP_DY)) {
    139     fprintf (stderr, "cannot find CHIP_DY in header of %s\n", filename);
    140     fclose (f);
    141     return NULL;
     103  // the simple files from Eddie have no internal metadata describing the corrections,
     104  // so they must be manually encoded
     105
     106  // hard-wired values which describe the ubercal analysis
     107  if (NO_METADATA) {
     108    NFILTER = 5;
     109    NSEASON = 4;
     110    NCHIP_X = 8;
     111    NCHIP_Y = 8;
     112    NCELL_X = 2;
     113    NCELL_Y = 2;
     114    CHIP_DX = 4880;
     115    CHIP_DY = 4864;
     116  } else {
     117    if (!gfits_scan (&header, "NSEASON", "%d", 1, &NSEASON)) {
     118      fprintf (stderr, "cannot find NSEASON in header of %s\n", filename);
     119      fclose (f);
     120      return NULL;
     121    }
     122    if (!gfits_scan (&header, "NFILTER", "%d", 1, &NFILTER)) {
     123      fprintf (stderr, "cannot find NFILTER in header of %s\n", filename);
     124      fclose (f);
     125      return NULL;
     126    }
     127    if (!gfits_scan (&header, "NCHIP_X", "%d", 1, &NCHIP_X)) {
     128      fprintf (stderr, "cannot find NCHIP_X in header of %s\n", filename);
     129      fclose (f);
     130      return NULL;
     131    }
     132    if (!gfits_scan (&header, "NCHIP_Y", "%d", 1, &NCHIP_Y)) {
     133      fprintf (stderr, "cannot find NCHIP_Y in header of %s\n", filename);
     134      fclose (f);
     135      return NULL;
     136    }
     137    if (!gfits_scan (&header, "NCELL_X", "%d", 1, &NCELL_X)) {
     138      fprintf (stderr, "cannot find NCELL_X in header of %s\n", filename);
     139      fclose (f);
     140      return NULL;
     141    }
     142    if (!gfits_scan (&header, "NCELL_Y", "%d", 1, &NCELL_Y)) {
     143      fprintf (stderr, "cannot find NCELL_Y in header of %s\n", filename);
     144      fclose (f);
     145      return NULL;
     146    }
     147    if (!gfits_scan (&header, "CHIP_DX", "%d", 1, &CHIP_DX)) {
     148      fprintf (stderr, "cannot find CHIP_DX in header of %s\n", filename);
     149      fclose (f);
     150      return NULL;
     151    }
     152    if (!gfits_scan (&header, "CHIP_DY", "%d", 1, &CHIP_DY)) {
     153      fprintf (stderr, "cannot find CHIP_DY in header of %s\n", filename);
     154      fclose (f);
     155      return NULL;
     156    }
    142157  }
    143158
    … …  
    146161  ALLOCATE (flatcorrTable->tstop,  e_time, NSEASON);
    147162
     163  // ubercal hard-coded values
     164  char filters_uc[5][3] = {"g", "r", "i", "z", "y"};
     165  double tstart_uc[] = {55000.0, 55296.0, 55327.0, 55662.0};
     166  double tstop_uc[]  = {55296.0, 55327.0, 55662.0, 60000.0};
     167
     168  double mjdstart, mjdstop;
    148169  for (i = 0; i < NSEASON; i++) {
    149     double mjdvalue;
    150     char name[9];
    151     snprintf (name, 9, "TS0_%04d", i);
    152     if (!gfits_scan (&header, name, "%lf", 1, &mjdvalue)) {
    153       fprintf (stderr, "cannot find %s in header of %s\n", name, filename);
    154       fclose (f);
    155       return NULL;
    156     }
    157     flatcorrTable->tstart[i] = ohana_mjd_to_sec(mjdvalue);
    158 
    159     snprintf (name, 9, "TS1_%04d", i);
    160     if (!gfits_scan (&header, name, "%lf", 1, &mjdvalue)) {
    161       fprintf (stderr, "cannot find %s in header of %s\n", name, filename);
    162       fclose (f);
    163       return NULL;
    164     }
    165     flatcorrTable->tstop[i]  = ohana_mjd_to_sec(mjdvalue);
     170    if (NO_METADATA) {
     171      mjdstart = tstart_uc[i];
     172      mjdstop  = tstop_uc[i];
     173  } else {
     174      char name[9];
     175      snprintf (name, 9, "TS0_%04d", i);
     176      if (!gfits_scan (&header, name, "%lf", 1, &mjdstart)) {
     177        fprintf (stderr, "cannot find %s in header of %s\n", name, filename);
     178        fclose (f);
     179        return NULL;
     180      }
     181
     182      snprintf (name, 9, "TS1_%04d", i);
     183      if (!gfits_scan (&header, name, "%lf", 1, &mjdstop)) {
     184        fprintf (stderr, "cannot find %s in header of %s\n", name, filename);
     185        fclose (f);
     186        return NULL;
     187      }
     188    }
     189    flatcorrTable->tstart[i] = ohana_mjd_to_sec(mjdstart);
     190    flatcorrTable->tstop[i]  = ohana_mjd_to_sec(mjdstop);
    166191  }
    167192
    … …  
    216241    if (!gfits_load_header (f, &header)) return (NULL);
    217242
    218     if (!gfits_scan (&header, "FILTER", "%s", 1, filter)) {
    219       fprintf (stderr, "cannot find FILTER in header of %s\n", filename);
    220       fclose (f);
    221       return NULL;
    222     }
    223 
     243    if (NO_METADATA) {
     244      strcpy (filter, filters_uc[nfilter]);
     245    } else {
     246      if (!gfits_scan (&header, "FILTER", "%s", 1, filter)) {
     247        fprintf (stderr, "cannot find FILTER in header of %s\n", filename);
     248        fclose (f);
     249        return NULL;
     250      }
     251    }
     252   
    224253    // read the fits table bytes
    225254    double *offset64 = NULL;
    226255    float *offset32 = NULL;
    227256    int use32 = TRUE;
    228 
     257   
    229258    if (!gfits_fread_matrix (f, &matrix, &header)) return (NULL);
    230259    switch (header.bitpix) {
    … …  
    242271    }
    243272
     273    // XXX the initial hacked-together table from Eddie is missing the last 2 elements.  they should be zero
     274    if (NO_METADATA) {
     275      matrix.Naxis[0] += 2;
     276      REALLOCATE (offset64, double, matrix.Naxis[0]);
     277      matrix.buffer = (char *) offset64;
     278    }
     279
    244280    for (nseason = 0; nseason < NSEASON; nseason++) { // seasons
    245281      for (iy = 0; iy < NCHIP_Y; iy++) { // y-chip
    … …  
    260296          flatcorrTable->image[Nimage].tstop  = flatcorrTable->tstop[nseason];
    261297         
     298          int seq_full = -1; // sequence number within the full table (all filters concatenated together)
     299          int seq_filt = -1; // sequence number for just this filter
     300
     301          // XXX we should have a NCHIP_X * NCHIP_Y array of x and y parity values
     302
    262303          // This enforces a 180 chip rotation for XY3n - XY7n & is only known to be valid for GPC1 (the XYnn names as well)
    263304          for (iyc = 0; iyc < NCELL_Y; iyc++) {
    264305            for (ixc = 0; ixc < NCELL_X; ixc++) {
    265               if (ix > 3) {
    266                 // XXX I think this is backwards : chips ix < 4 should be flipped
    267                 seq = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + NCELL_Y - 1 - iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + NCELL_X - 1 - ixc);
     306              if (ix < 4) {
     307                // chips ix < 4 should be flipped
     308                seq_full = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + NCELL_Y - 1 - iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + NCELL_X - 1 - ixc);
     309                seq_filt =                                                   nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + NCELL_Y - 1 - iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + NCELL_X - 1 - ixc);
     310                // = .... + y_parity*iyc - (y_parity - 1)*(NCELL_Y - 1) / 2
    268311              } else {
    269                 seq = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + ixc);
     312                seq_full = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + ixc);
     313                seq_filt =                                                   nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + ixc);
    270314              }
    271               assert (!flatcorrTable->corr[seq].ID);
    272               flatcorrTable->corr[seq].x = ixc;
    273               flatcorrTable->corr[seq].y = iyc;
     315              assert (seq_full > -1);
     316              assert (seq_filt > -1);
     317              assert (seq_full < flatcorrTable->Ncorr);
     318              if (matrix.Naxes == 1) {
     319                assert (seq_filt < matrix.Naxis[0]);
     320              }
     321              if (matrix.Naxes == 2) {
     322                assert (seq_filt < matrix.Naxis[0]*matrix.Naxis[1]);
     323              }
     324              assert (!flatcorrTable->corr[seq_full].ID);
     325              flatcorrTable->corr[seq_full].x = ixc;
     326              flatcorrTable->corr[seq_full].y = iyc;
    274327              if (use32) {
    275                 flatcorrTable->corr[seq].offset = offset32[seq];
     328                flatcorrTable->corr[seq_full].offset = offset32[seq_filt];
    276329              } else {
    277                 flatcorrTable->corr[seq].offset = offset64[seq];
     330                flatcorrTable->corr[seq_full].offset = offset64[seq_filt];
    278331              }
    279               flatcorrTable->corr[seq].ID = corrID;
     332              flatcorrTable->corr[seq_full].ID = corrID;
    280333            }
    281334          }
    … …  
    295348  return zpts;
    296349}
    297 
    298 ZptTable *load_zpt_ubercal_nometadata (char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable) {
    299 
    300   int i, nfilter, nseason, ix, iy, ixc, iyc, Ncol, seq;
    301   off_t Nrow;
    302   char type[16];
    303   int Nzpts, NZPTS;
    304   ZptTable *zpts;
    305   Header header;
    306   Header theader;
    307   Matrix matrix;
    308   FTable ftable;
    309 
    310   // hard-wired values which describe the ubercal analysis
    311   int NFILTER = 5;
    312   int NSEASON = 4;
    313   int NCHIP_X = 8;
    314   int NCHIP_Y = 8;
    315   int NCELL_X = 2;
    316   int NCELL_Y = 2;
    317   int CHIP_DX = 4880;
    318   int CHIP_DY = 4864;
    319 
    320   flatcorrTable->Nseason = NSEASON;
    321   ALLOCATE (flatcorrTable->tstart, e_time, NSEASON);
    322   ALLOCATE (flatcorrTable->tstop,  e_time, NSEASON);
    323 
    324   char filter[] = {'g', 'r', 'i', 'z', 'y'};
    325 
    326   double tstart_mjd[] = {55000.0, 55296.0, 55327.0, 55662.0};
    327   double tstop_mjd[]  = {55296.0, 55327.0, 55662.0, 60000.0};
    328 
    329   for (i = 0; i < NSEASON; i++) {
    330     flatcorrTable->tstart[i] = ohana_mjd_to_sec(tstart_mjd[i]);
    331     flatcorrTable->tstop[i]  = ohana_mjd_to_sec(tstop_mjd[i]);
    332   }
    333 
    334   FILE *f = fopen (filename, "r");
    335   if (!f) {
    336     fprintf (stderr, "ERROR: cannot open zpt table file %s\n", filename);
    337     exit (1);
    338   }
    339 
    340   /* load in PHU segment (ignore) */
    341   if (!gfits_fread_header (f, &header)) {
    342     if (VERBOSE) fprintf (stderr, "can't read Flat Correction header\n");
    343     fclose (f);
    344     return (NULL);
    345   }
    346   if (!gfits_fread_matrix (f, &matrix, &header)) {
    347     if (VERBOSE) fprintf (stderr, "can't read Flat Correction matrix\n");
    348     gfits_free_header (&header);
    349     fclose (f);
    350     return (NULL);
    351   }
    352 
    353   Nzpts = 0;
    354   NZPTS = 0;
    355   ALLOCATE (zpts, ZptTable, NZPTS);
    356 
    357   // this function would be better if we read the list of filters, seasons, and the dimensions from the header
    358   // for current testing, make fake smfs that correspond to specific chips, filter, and mjd ranges?
    359 
    360   // we have 5 filters, and 4 flat-field correction sets for each
    361   flatcorrTable->Ncorr = NFILTER*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y;
    362   flatcorrTable->Nimage = NFILTER*NSEASON*NCHIP_X*NCHIP_Y;
    363 
    364   ALLOCATE (flatcorrTable->corr, FlatCorrection, flatcorrTable->Ncorr);
    365   ALLOCATE (flatcorrTable->image, FlatCorrectionImage, flatcorrTable->Nimage);
    366   memset (flatcorrTable->corr, 0, flatcorrTable->Ncorr*sizeof(FlatCorrection));
    367 
    368   int corrID = 1;
    369   int Nimage = 0;
    370   ftable.header = &theader;
    371   for (nfilter = 0; nfilter < NFILTER; nfilter++) {
    372     // load data for this header
    373     if (!gfits_load_header (f, &theader)) return (NULL);
    374 
    375     // read the fits table bytes
    376     if (!gfits_fread_ftable_data (f, &ftable, FALSE)) return (NULL);
    377  
    378     // skip over remaining bytes in data segment
    379     fseeko (f, ftable.datasize - ftable.validsize, SEEK_CUR);
    380 
    381     // need to create and assign to flat-field correction
    382     double *mjd = gfits_get_bintable_column_data (&theader, &ftable, "mjd_obs", type, &Nrow, &Ncol);
    383     assert (!strcmp(type, "double"));
    384 
    385     double *zp = gfits_get_bintable_column_data (&theader, &ftable, "zp", type, &Nrow, &Ncol);
    386     assert (!strcmp(type, "double"));
    387      
    388     // float *zperr = gfits_get_bintable_column_data (&theader, &ftable, "resid", type, &Nrow, &Ncol);
    389     // assert (!strcmp(type, "float"));
    390      
    391     NZPTS += Nrow;
    392     REALLOCATE (zpts, ZptTable, NZPTS);
    393     for (i = 0; i < Nrow; i++) {
    394       zpts[i+Nzpts].time = ohana_mjd_to_sec(mjd[i]);
    395       zpts[i+Nzpts].zpt = zp[i];
    396       // zpts[i+Nzpts].zpt_err = zperr[i];
    397     }
    398     Nzpts += Nrow;
    399 
    400     // the image contains the flat-field corrections
    401 
    402     // load data for this header
    403     if (!gfits_load_header (f, &header)) return (NULL);
    404 
    405     // read the fits table bytes
    406     if (!gfits_fread_matrix (f, &matrix, &header)) return (NULL);
    407     assert (header.bitpix == -64); // hardwired as a double
    408     double *offset = (double *) matrix.buffer;
    409 
    410     for (nseason = 0; nseason < NSEASON; nseason++) { // seasons
    411       for (iy = 0; iy < NCHIP_Y; iy++) { // y-chip
    412         for (ix = 0; ix < NCHIP_X; ix++) { // x-chip
    413 
    414           // photcode name
    415           char photname[64];
    416           snprintf (photname, 64, "GPC1.%c.XY%d%d", filter[nfilter], ix, iy);
    417           // note that the XY00, XY07, etc, chips will have photcode values of 0
    418 
    419           flatcorrTable->image[Nimage].photcode = GetPhotcodeCodebyName(photname);
    420           flatcorrTable->image[Nimage].Nx = NCELL_X;
    421           flatcorrTable->image[Nimage].Ny = NCELL_Y;
    422           flatcorrTable->image[Nimage].ID = corrID;
    423           flatcorrTable->image[Nimage].DX = CHIP_DX;
    424           flatcorrTable->image[Nimage].DY = CHIP_DY;
    425           flatcorrTable->image[Nimage].tstart = flatcorrTable->tstart[nseason];
    426           flatcorrTable->image[Nimage].tstop  = flatcorrTable->tstop[nseason];
    427          
    428           for (iyc = 0; iyc < NCELL_Y; iyc++) {
    429             for (ixc = 0; ixc < NCELL_X; ixc++) {
    430               if (ix > 3) {
    431                 seq = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + NCELL_Y - 1 - iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + NCELL_X - 1 - ixc);
    432               } else {
    433                 seq = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + ixc);
    434               }
    435               assert (!flatcorrTable->corr[seq].ID);
    436               flatcorrTable->corr[seq].x = ixc;
    437               flatcorrTable->corr[seq].y = iyc;
    438               flatcorrTable->corr[seq].offset = offset[seq];
    439               flatcorrTable->corr[seq].ID = corrID;
    440             }
    441           }
    442           corrID ++;
    443           Nimage ++;
    444         }
    445       }
    446     }         
    447   }
    448 
    449   /*** convert from corr,image format to offsets ***/
    450   FlatCorrectionInternal (flatcorrTable);
    451 
    452   fprintf (stderr, "loaded %d zero points\n", Nzpts);
    453 
    454   *nzpts = Nzpts;
    455   return zpts;
    456 }

  • branches/eam_branches/ipp-20111122/Ohana/src/uniphot/src/setphot.c

    r33156 r33180  
    1919
    2020  if (UBERCAL) {
    21     if (NO_METADATA) {
    22       // the simple files from Eddie have no internal metadata describing the corrections,
    23       // so they must be manually encoded
    24       zpts = load_zpt_ubercal_nometadata (argv[1], &Nzpts, &flatcorrTable);
    25     } else {
    26       zpts = load_zpt_ubercal (argv[1], &Nzpts, &flatcorrTable);
    27     }
     21    zpts = load_zpt_ubercal (argv[1], &Nzpts, &flatcorrTable);
    2822
    2923    char flatcorrfile[64];
Note: See TracChangeset for help on using the changeset viewer.