Index: branches/eam_branches/relphot.20210521/src/GridIO.c =================================================================== --- branches/eam_branches/relphot.20210521/src/GridIO.c (revision 41610) +++ branches/eam_branches/relphot.20210521/src/GridIO.c (revision 41615) @@ -1,5 +1,28 @@ # include "relphot.h" -int GridCorrectionSave (char *filename) { +int GridCorrectionSave () { + + if (!GRID_ZEROPT) return TRUE; + + char filename[1024], uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + + snprintf_nowarn (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + + snprintf (filename, 1024, "%s/gridcorr.mean.%s.fits", CATDIR, uniquer); + if (!GridCorrectionSaveFile (filename, GRID_MEAN)) return FALSE; + GRID_MEANFILE = strcreate (filename); // save in the global so we can pass to e.g., reload_catalogs + + snprintf (filename, 1024, "%s/gridcorr.stdev.%s.fits", CATDIR, uniquer); + if (!GridCorrectionSaveFile (filename, GRID_STDEV)) return FALSE; + + snprintf (filename, 1024, "%s/gridcorr.npts.%s.fits", CATDIR, uniquer); + if (!GridCorrectionSaveFile (filename, GRID_NPTS)) return FALSE; + + return TRUE; +} + +int GridCorrectionSaveFile (char *filename, GridOutputMode mode) { /* open file for input */ @@ -27,9 +50,11 @@ while (TRUE) { - // getGridCorr() takes a photcode and returns the *next* valid grid correction pointer (start at -1) - GridCorrectionType *GridCorr = getGridCorr (&Nlast); + // getGridCorreNext() takes a photcode and returns the *next* valid grid correction pointer (start at -1) + GridCorrectionType *GridCorr = getGridCorrNext (&Nlast); if (!GridCorr) break; int filter = GetPhotcodeEquivCodebyCode (GridCorr->photcode); + char *filtname = GetPhotcodeNamebyCode (filter); // reference, do not free + char *photname = GetPhotcodeNamebyCode (GridCorr->photcode); // reference, do not free Header header; @@ -37,6 +62,6 @@ gfits_init_header (&header); - header.bitpix = -32; - header.Naxes = 2; + header.bitpix = (mode == GRID_NPTS) ? 32 : -32; + header.Naxes = 2; header.Naxis[0] = GridCorr->Nx; header.Naxis[1] = GridCorr->Ny; @@ -48,5 +73,5 @@ // create header & matrix if (!gfits_modify (&header, "PHOTCODE", "%d", 1, GridCorr->photcode)) myAbort ("failed to set PHOTCODE"); - if (!gfits_modify (&header, "FILTER", "%d", 1, filter)) myAbort ("failed to set FILTER"); + if (!gfits_modify (&header, "FILTER", "%s", 1, filtname)) myAbort ("failed to set FILTER"); if (!gfits_modify (&header, "SEASON", "%d", 1, 0)) myAbort ("failed to set SEASON"); // XXX hard-wired for now if (!gfits_modify (&header, "X_CHIP", "%d", 1, -1)) myAbort ("failed to set X_CHIP"); @@ -58,12 +83,27 @@ int Nx = GridCorr->Nx; - float *value = (float *) matrix.buffer; + float *fvalue = (float *) matrix.buffer; + int *ivalue = (int *) matrix.buffer; + for (int iy = 0; iy < GridCorr->Ny; iy++) { for (int ix = 0; ix < GridCorr->Nx; ix++) { - value[ix + iy*Nx] = GridCorr->Mgrid[ix][iy]; + switch (mode) { + case GRID_MEAN: fvalue[ix + iy*Nx] = GridCorr-> Mgrid[ix][iy]; break; + case GRID_STDEV: fvalue[ix + iy*Nx] = GridCorr->dMgrid[ix][iy]; break; + case GRID_NPTS: ivalue[ix + iy*Nx] = GridCorr->nMgrid[ix][iy]; break; + } } } + char extname[85]; + switch (mode) { + // XXX If we add seasons, add here: + case GRID_MEAN: snprintf (extname, 85, "%s.%s", photname, "MEAN"); break; + case GRID_STDEV: snprintf (extname, 85, "%s.%s", photname, "STDEV"); break; + case GRID_NPTS: snprintf (extname, 85, "%s.%s", photname, "NPTS"); break; + } + gfits_modify (&header, "EXTNAME", "%s", 1, extname); + gfits_primary_to_extended (&header, "IMAGE", "Image Extension"); gfits_fwrite_header (f, &header); @@ -78,10 +118,16 @@ } -# if (0) -CamPhotomCorrection *CamPhotomCorrectionLoad (char *filename) { +void GridCorrectionLoad (char *filename) { - int i, ix, iy, season, filter; + int photcode; - Header *phu = gfits_alloc_header(); + Header phu; + Header header; + Matrix matrix; + + if (!GRID_ZEROPT) return; + + // allocate the basic containing structures: + initGridBins (); /* open file for input */ @@ -89,74 +135,47 @@ if (f == (FILE *) NULL) { gprint (GP_ERR, "can't open file for read : %s\n", filename); - return FALSE; + return; } // read the PHU header - if (!gfits_load_header (f, phu)) return FALSE; + if (!gfits_fread_header (f, &phu)) return; - int Nbytes = gfits_data_size (phu); + // move to the first extension + int Nbytes = gfits_data_size (&phu); fseeko (f, Nbytes, SEEK_CUR); - int Nseason, Nfilter, Nx, Ny, dX, dY, NxCCD, NyCCD; - if (!gfits_scan (phu, "NSEASON", "%d", 1, &Nseason)) return FALSE; - if (!gfits_scan (phu, "NFILTER", "%d", 1, &Nfilter)) return FALSE; - if (!gfits_scan (phu, "NX", "%d", 1, &Nx)) return FALSE; - if (!gfits_scan (phu, "NY", "%d", 1, &Ny)) return FALSE; - if (!gfits_scan (phu, "DX", "%d", 1, &dX)) return FALSE; - if (!gfits_scan (phu, "DY", "%d", 1, &dY)) return FALSE; - if (!gfits_scan (phu, "NX_CHIP", "%d", 1, &NxCCD)) return FALSE; - if (!gfits_scan (phu, "NY_CHIP", "%d", 1, &NyCCD)) return FALSE; + gfits_free_header (&phu); - CamPhotomCorrection *cam = CamPhotomCorrectionAlloc (Nx, Ny, Nfilter, Nseason); - cam->phu = phu; + while (gfits_fread_header (f, &header)) { + if (!gfits_fread_matrix (f, &matrix, &header)) myAbort ("failed to read matrix"); - cam->dX = dX; // superpixel sampling, x-dir - cam->dY = dY; // superpixel sampling, y-dir + if (!gfits_scan (&header, "PHOTCODE", "%d", 1, &photcode)) myAbort ("failed to get PHOTCODE"); + GridCorrectionType *GridCorr = getGridCorrByCode (photcode); + GridCorr->photcode = photcode; - cam->NxCCD = NxCCD; // pixels per chip - cam->NyCCD = NyCCD; // pixels per chip + if (!gfits_scan (&header, "DX_CHIP", "%f", 1, &GridCorr->dX)) myAbort ("failed to get DX_CHIP"); + if (!gfits_scan (&header, "DY_CHIP", "%f", 1, &GridCorr->dY)) myAbort ("failed to get DY_CHIP"); + if (!gfits_scan (&header, "NX_CHIP", "%d", 1, &GridCorr->Nx)) myAbort ("failed to get NX_CHIP"); + if (!gfits_scan (&header, "NY_CHIP", "%d", 1, &GridCorr->Ny)) myAbort ("failed to get NY_CHIP"); + + GridCorr->dMgrid = NULL; + GridCorr->nMgrid = NULL; + ALLOCATE (GridCorr->Mgrid, float *, GridCorr->Nx); + for (int ix = 0; ix < GridCorr->Nx; ix++) { + ALLOCATE (GridCorr->Mgrid[ix], float, GridCorr->Ny); + } - for (i = 0; i < Nseason; i++) { - char field[256]; - double value; + int Nx = GridCorr->Nx; + + float *fvalue = (float *) matrix.buffer; - snprintf (field, 256, "TS_%03d", i); - if (!gfits_scan (phu, field, "%lf", 1, &value)) return FALSE; - cam->tstart[i] = ohana_mjd_to_sec(value); - - snprintf (field, 256, "TE_%03d", i); - if (!gfits_scan (phu, field, "%lf", 1, &value)) return FALSE; - cam->tstop[i] = ohana_mjd_to_sec(value); + for (int iy = 0; iy < GridCorr->Ny; iy++) { + for (int ix = 0; ix < GridCorr->Nx; ix++) { + GridCorr-> Mgrid[ix][iy] = fvalue[ix + iy*Nx]; + } + } + gfits_free_header (&header); + gfits_free_matrix (&matrix); } - - while (TRUE) { - // load data for this header : if not found, assume we hit the end of the file - Header *header = gfits_alloc_header(); - if (!gfits_load_header (f, header)) { - free (header); - break; - } - - if (!gfits_scan (header, "FILTER", "%d", 1, &filter)) myAbort ("failed to find FILTER"); - if (!gfits_scan (header, "SEASON", "%d", 1, &season)) myAbort ("failed to find SEASON"); - if (!gfits_scan (header, "X_CHIP", "%d", 1, &ix)) myAbort ("failed to find X_CHIP"); - if (!gfits_scan (header, "Y_CHIP", "%d", 1, &iy)) myAbort ("failed to find Y_CHIP"); - // XXX NOTE: astroflat.20150209.fits had ix and iy backwards in header - // double-check that the new flats are OK - - Matrix *matrix = gfits_alloc_matrix(); - if (!gfits_load_matrix (f, matrix, header)) myAbort ("failed to read matrix"); - - int index = ix + iy*cam->Nx + filter*cam->Nchips + season*cam->Nflats; - myAssert (index >= 0, "index too small"); - myAssert (index < cam->Nvalues, "index too big"); - - myAssert (!cam->matrix[index], "entry already assigned?"); - - cam->matrix[index] = matrix; - cam->header[index] = header; - } - return cam; + gfits_free_header (&header); } -# endif -