IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 25, 2021, 5:00:52 PM (5 years ago)
Author:
eugene
Message:

save grid correction into well-defined location; pass to relphot_client if needed

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/relphot.20210521/src/GridIO.c

    r41610 r41615  
    11# include "relphot.h"
    22
    3 int GridCorrectionSave (char *filename) {
     3int GridCorrectionSave () {
     4
     5  if (!GRID_ZEROPT) return TRUE;
     6
     7  char filename[1024], uniquer[12];
     8  int TIME = time(NULL);
     9  int PID = getpid();
     10
     11  snprintf_nowarn (uniquer, 12, "%05d.%05d", PID, TIME % 100000);
     12
     13  snprintf (filename, 1024, "%s/gridcorr.mean.%s.fits", CATDIR, uniquer);
     14  if (!GridCorrectionSaveFile (filename, GRID_MEAN)) return FALSE;
     15  GRID_MEANFILE = strcreate (filename);  // save in the global so we can pass to e.g., reload_catalogs
     16
     17  snprintf (filename, 1024, "%s/gridcorr.stdev.%s.fits", CATDIR, uniquer);
     18  if (!GridCorrectionSaveFile (filename, GRID_STDEV)) return FALSE;
     19
     20  snprintf (filename, 1024, "%s/gridcorr.npts.%s.fits", CATDIR, uniquer);
     21  if (!GridCorrectionSaveFile (filename, GRID_NPTS)) return FALSE;
     22
     23  return TRUE;
     24}
     25
     26int GridCorrectionSaveFile (char *filename, GridOutputMode mode) {
    427
    528  /* open file for input */
    … …  
    2750
    2851  while (TRUE) {
    29     // getGridCorr() takes a photcode and returns the *next* valid grid correction pointer (start at -1)
    30     GridCorrectionType *GridCorr = getGridCorr (&Nlast);
     52    // getGridCorreNext() takes a photcode and returns the *next* valid grid correction pointer (start at -1)
     53    GridCorrectionType *GridCorr = getGridCorrNext (&Nlast);
    3154    if (!GridCorr) break;
    3255
    3356    int filter = GetPhotcodeEquivCodebyCode (GridCorr->photcode);
     57    char *filtname = GetPhotcodeNamebyCode (filter); // reference, do not free
     58    char *photname = GetPhotcodeNamebyCode (GridCorr->photcode); // reference, do not free
    3459
    3560    Header header;
    … …  
    3762
    3863    gfits_init_header (&header);
    39     header.bitpix   = -32;
    40     header.Naxes    =   2;
     64    header.bitpix   = (mode == GRID_NPTS) ? 32 : -32;
     65    header.Naxes    = 2;
    4166    header.Naxis[0] = GridCorr->Nx;
    4267    header.Naxis[1] = GridCorr->Ny;
    … …  
    4873    // create header & matrix
    4974    if (!gfits_modify (&header, "PHOTCODE", "%d", 1, GridCorr->photcode)) myAbort ("failed to set PHOTCODE");
    50     if (!gfits_modify (&header, "FILTER",   "%d", 1, filter))             myAbort ("failed to set FILTER");
     75    if (!gfits_modify (&header, "FILTER",   "%s", 1, filtname))           myAbort ("failed to set FILTER");
    5176    if (!gfits_modify (&header, "SEASON",   "%d", 1, 0))                  myAbort ("failed to set SEASON"); // XXX hard-wired for now
    5277    if (!gfits_modify (&header, "X_CHIP",   "%d", 1, -1))                 myAbort ("failed to set X_CHIP");
    … …  
    5883
    5984    int Nx = GridCorr->Nx;
    60     float *value = (float *) matrix.buffer;
    6185
     86    float *fvalue = (float *) matrix.buffer;
     87    int   *ivalue = (int   *) matrix.buffer;
     88   
    6289    for (int iy = 0; iy < GridCorr->Ny; iy++) {
    6390      for (int ix = 0; ix < GridCorr->Nx; ix++) {
    64         value[ix + iy*Nx] = GridCorr->Mgrid[ix][iy];
     91        switch (mode) {
     92          case GRID_MEAN:  fvalue[ix + iy*Nx] = GridCorr-> Mgrid[ix][iy]; break;
     93          case GRID_STDEV: fvalue[ix + iy*Nx] = GridCorr->dMgrid[ix][iy]; break;
     94          case GRID_NPTS:  ivalue[ix + iy*Nx] = GridCorr->nMgrid[ix][iy]; break;
     95        }
    6596      }
    6697    }
    6798   
     99    char extname[85];
     100    switch (mode) {
     101      // XXX If we add seasons, add here:
     102      case GRID_MEAN:  snprintf (extname, 85, "%s.%s", photname, "MEAN");  break;
     103      case GRID_STDEV: snprintf (extname, 85, "%s.%s", photname, "STDEV"); break;
     104      case GRID_NPTS:  snprintf (extname, 85, "%s.%s", photname, "NPTS");  break;
     105    }
     106    gfits_modify (&header, "EXTNAME", "%s", 1, extname);
     107
    68108    gfits_primary_to_extended (&header, "IMAGE", "Image Extension");
    69109    gfits_fwrite_header  (f, &header);
    … …  
    78118}
    79119
    80 # if (0)
    81 CamPhotomCorrection *CamPhotomCorrectionLoad (char *filename) {
     120void GridCorrectionLoad (char *filename) {
    82121
    83   int i, ix, iy, season, filter;
     122  int photcode;
    84123
    85   Header *phu = gfits_alloc_header();
     124  Header phu;
     125  Header header;
     126  Matrix matrix;
     127
     128  if (!GRID_ZEROPT) return;
     129
     130  // allocate the basic containing structures:
     131  initGridBins ();
    86132
    87133  /* open file for input */
    … …  
    89135  if (f == (FILE *) NULL) {
    90136    gprint (GP_ERR, "can't open file for read : %s\n", filename);
    91     return FALSE;
     137    return;
    92138  }
    93139
    94140  // read the PHU header
    95   if (!gfits_load_header (f, phu)) return FALSE;
     141  if (!gfits_fread_header (f, &phu)) return;
    96142 
    97   int Nbytes = gfits_data_size (phu);
     143  // move to the first extension
     144  int Nbytes = gfits_data_size (&phu);
    98145  fseeko (f, Nbytes, SEEK_CUR);
    99146
    100   int Nseason, Nfilter, Nx, Ny, dX, dY, NxCCD, NyCCD;
    101   if (!gfits_scan (phu, "NSEASON", "%d", 1, &Nseason)) return FALSE;
    102   if (!gfits_scan (phu, "NFILTER", "%d", 1, &Nfilter)) return FALSE;
    103   if (!gfits_scan (phu, "NX",      "%d", 1, &Nx))      return FALSE;
    104   if (!gfits_scan (phu, "NY",      "%d", 1, &Ny))      return FALSE;
    105   if (!gfits_scan (phu, "DX",      "%d", 1, &dX))      return FALSE;
    106   if (!gfits_scan (phu, "DY",      "%d", 1, &dY))      return FALSE;
    107   if (!gfits_scan (phu, "NX_CHIP", "%d", 1, &NxCCD))   return FALSE;
    108   if (!gfits_scan (phu, "NY_CHIP", "%d", 1, &NyCCD))   return FALSE;
     147  gfits_free_header (&phu);
    109148
    110   CamPhotomCorrection *cam = CamPhotomCorrectionAlloc (Nx, Ny, Nfilter, Nseason);
    111   cam->phu = phu;
     149  while (gfits_fread_header (f, &header)) {
     150    if (!gfits_fread_matrix (f, &matrix, &header))                         myAbort ("failed to read matrix");
    112151
    113   cam->dX      = dX;       // superpixel sampling, x-dir
    114   cam->dY      = dY;       // superpixel sampling, y-dir
     152    if (!gfits_scan (&header, "PHOTCODE", "%d", 1, &photcode))    myAbort ("failed to get PHOTCODE");
     153    GridCorrectionType *GridCorr = getGridCorrByCode (photcode);
     154    GridCorr->photcode = photcode;
    115155
    116   cam->NxCCD   = NxCCD;    // pixels per chip
    117   cam->NyCCD   = NyCCD;    // pixels per chip
     156    if (!gfits_scan (&header, "DX_CHIP",  "%f", 1, &GridCorr->dX)) myAbort ("failed to get DX_CHIP");
     157    if (!gfits_scan (&header, "DY_CHIP",  "%f", 1, &GridCorr->dY)) myAbort ("failed to get DY_CHIP");
     158    if (!gfits_scan (&header, "NX_CHIP",  "%d", 1, &GridCorr->Nx)) myAbort ("failed to get NX_CHIP");
     159    if (!gfits_scan (&header, "NY_CHIP",  "%d", 1, &GridCorr->Ny)) myAbort ("failed to get NY_CHIP");
     160   
     161    GridCorr->dMgrid = NULL;
     162    GridCorr->nMgrid = NULL;
     163    ALLOCATE (GridCorr->Mgrid, float *, GridCorr->Nx);
     164    for (int ix = 0; ix < GridCorr->Nx; ix++) {
     165      ALLOCATE (GridCorr->Mgrid[ix], float, GridCorr->Ny);
     166    }     
    118167
    119   for (i = 0; i < Nseason; i++) {
    120     char field[256];
    121     double value;
     168    int Nx = GridCorr->Nx;
     169
     170    float *fvalue = (float *) matrix.buffer;
    122171   
    123     snprintf (field, 256, "TS_%03d", i);
    124     if (!gfits_scan (phu, field, "%lf", 1, &value)) return FALSE;
    125     cam->tstart[i] = ohana_mjd_to_sec(value);
    126    
    127     snprintf (field, 256, "TE_%03d", i);
    128     if (!gfits_scan (phu, field, "%lf", 1, &value)) return FALSE;
    129     cam->tstop[i] = ohana_mjd_to_sec(value);
     172    for (int iy = 0; iy < GridCorr->Ny; iy++) {
     173      for (int ix = 0; ix < GridCorr->Nx; ix++) {
     174        GridCorr-> Mgrid[ix][iy] = fvalue[ix + iy*Nx];
     175      }
     176    }
     177    gfits_free_header (&header);
     178    gfits_free_matrix (&matrix);
    130179  }
    131 
    132   while (TRUE) {
    133     // load data for this header : if not found, assume we hit the end of the file
    134     Header *header = gfits_alloc_header();
    135     if (!gfits_load_header (f, header)) {
    136       free (header);
    137       break;
    138     }
    139    
    140     if (!gfits_scan (header, "FILTER", "%d", 1, &filter)) myAbort ("failed to find FILTER");
    141     if (!gfits_scan (header, "SEASON", "%d", 1, &season)) myAbort ("failed to find SEASON");
    142     if (!gfits_scan (header, "X_CHIP", "%d", 1, &ix))     myAbort ("failed to find X_CHIP");
    143     if (!gfits_scan (header, "Y_CHIP", "%d", 1, &iy))     myAbort ("failed to find Y_CHIP");
    144     // XXX NOTE: astroflat.20150209.fits had ix and iy backwards in header
    145     // double-check that the new flats are OK
    146 
    147     Matrix *matrix = gfits_alloc_matrix();
    148     if (!gfits_load_matrix (f, matrix, header)) myAbort ("failed to read matrix");
    149 
    150     int index = ix + iy*cam->Nx + filter*cam->Nchips + season*cam->Nflats;
    151     myAssert (index >= 0, "index too small");
    152     myAssert (index < cam->Nvalues, "index too big");
    153 
    154     myAssert (!cam->matrix[index], "entry already assigned?");
    155 
    156     cam->matrix[index] = matrix;
    157     cam->header[index] = header;
    158   }
    159   return cam;
     180  gfits_free_header (&header);
    160181}
    161 # endif
    162 
Note: See TracChangeset for help on using the changeset viewer.