Index: /branches/eam_branches/relphot.20210521/include/relphot.h =================================================================== --- /branches/eam_branches/relphot.20210521/include/relphot.h (revision 41614) +++ /branches/eam_branches/relphot.20210521/include/relphot.h (revision 41615) @@ -66,4 +66,10 @@ PS1_w = 6, } PS1filters; + +typedef enum { + GRID_MEAN = 0, + GRID_STDEV = 1, + GRID_NPTS = 2, +} GridOutputMode; typedef struct { @@ -374,4 +380,5 @@ int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ char *SYNTH_ZERO_POINTS; +char *GRID_MEANFILE; // globals for parallel region operations @@ -574,6 +581,9 @@ void setMflatFromGrid PROTO((Catalog *catalog)); void initGridBins PROTO(()); -GridCorrectionType *getGridCorr PROTO((int *Nlast)); -int GridCorrectionSave PROTO((char *filename)); +GridCorrectionType *getGridCorrByCode PROTO((int code)); +GridCorrectionType *getGridCorrNext PROTO((int *Nlast)); +int GridCorrectionSave PROTO((void)); +int GridCorrectionSaveFile PROTO((char *filename, GridOutputMode mode)); +void GridCorrectionLoad PROTO((char *filename)); void initImageBins PROTO((Catalog *catalog, int Ncatalog, int doImageList)); void initImagesSubset PROTO((ImageSubset *input, off_t *line_number, off_t N)); Index: /branches/eam_branches/relphot.20210521/src/GridIO.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/GridIO.c (revision 41614) +++ /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 - Index: /branches/eam_branches/relphot.20210521/src/GridOps.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/GridOps.c (revision 41614) +++ /branches/eam_branches/relphot.20210521/src/GridOps.c (revision 41615) @@ -68,6 +68,6 @@ for (int ix = 0; ix < GridCorr[code]->Nx; ix++) { FREE (GridCorr[code]-> Mgrid[ix]); - FREE (GridCorr[code]->dMgrid[ix]); - FREE (GridCorr[code]->nMgrid[ix]); + if (GridCorr[code]->dMgrid) FREE (GridCorr[code]->dMgrid[ix]); + if (GridCorr[code]->nMgrid) FREE (GridCorr[code]->nMgrid[ix]); } @@ -85,5 +85,5 @@ } -GridCorrectionType *getGridCorr (int *Nlast) { +GridCorrectionType *getGridCorrNext (int *Nlast) { if (GridCorr == NULL) return NULL; @@ -99,4 +99,16 @@ *Nlast = i; } + return result; +} + +GridCorrectionType *getGridCorrByCode (int code) { + + if (GridCorr == NULL) return NULL; + + if (code >= NGridCorr) return NULL; + if (code < 0) return NULL; + + ALLOCATE(GridCorr[code], GridCorrectionType, 1) + GridCorrectionType *result = GridCorr[code]; return result; } @@ -135,7 +147,7 @@ ALLOCATE (GridCorr[code]->nMgrid, int *, NX_BIN); for (int ix = 0; ix < NX_BIN; ix++) { - ALLOCATE (GridCorr[code]-> Mgrid[ix], float, NX_BIN); - ALLOCATE (GridCorr[code]->dMgrid[ix], float, NX_BIN); - ALLOCATE (GridCorr[code]->nMgrid[ix], int, NX_BIN); + ALLOCATE (GridCorr[code]-> Mgrid[ix], float, NY_BIN); + ALLOCATE (GridCorr[code]->dMgrid[ix], float, NY_BIN); + ALLOCATE (GridCorr[code]->nMgrid[ix], int, NY_BIN); } } Index: /branches/eam_branches/relphot.20210521/src/args.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/args.c (revision 41614) +++ /branches/eam_branches/relphot.20210521/src/args.c (revision 41615) @@ -382,4 +382,5 @@ // GRID_ZEROPT is not valid for all cases, probably should be its own mode... GRID_ZEROPT = FALSE; + GRID_MEANFILE = NULL; // this is set by GridCorrectionSave() if ((N = get_argument (argc, argv, "-grid"))) { remove_argument (N, &argc, argv); @@ -602,4 +603,5 @@ FREE (REGION_FILE); FREE (SYNTH_ZERO_POINTS); + FREE (GRID_MEANFILE); FREE (IMAGE_TABLE); FREE (CATDIR); @@ -858,4 +860,14 @@ } + // load the GridCorrection file + GRID_ZEROPT = FALSE; + GRID_MEANFILE = NULL; + if ((N = get_argument (argc, argv, "-grid"))) { + GRID_ZEROPT = TRUE; + remove_argument (N, &argc, argv); + GRID_MEANFILE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + ImagSelect = FALSE; if ((N = get_argument (argc, argv, "-instmag"))) { @@ -916,4 +928,5 @@ FREE (SYNTH_ZERO_POINTS); + FREE (GRID_MEANFILE); FREE (BOUNDARY_TREE); FREE (UPDATE_CATFORMAT); Index: /branches/eam_branches/relphot.20210521/src/reload_catalogs.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/reload_catalogs.c (revision 41614) +++ /branches/eam_branches/relphot.20210521/src/reload_catalogs.c (revision 41615) @@ -264,4 +264,5 @@ if (BOUNDARY_TREE) { strextend (&command, "-boundary-tree %s", BOUNDARY_TREE); } if (SYNTH_ZERO_POINTS) { strextend (&command, "-synthphot-zpts %s", SYNTH_ZERO_POINTS); } + if (GRID_ZEROPT) { strextend (&command, "-grid %s", GRID_MEANFILE); } if (USE_BASIC_CHECK) { strextend (&command, "-basic-image-search"); } if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); } Index: /branches/eam_branches/relphot.20210521/src/relphot_client.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/relphot_client.c (revision 41614) +++ /branches/eam_branches/relphot.20210521/src/relphot_client.c (revision 41615) @@ -80,7 +80,11 @@ initImagesSubset (image, NULL, Nimage); + // load grid corrections here + GridCorrectionLoad (GRID_MEANFILE); + reload_catalogs (skylist, HOST_ID, HOSTDIR); freeImages ((char *)image); free (image); + freeGridBins (); client_logger_message ("updated catalogs\n"); relphot_client_free (sky, skylist); Index: /branches/eam_branches/relphot.20210521/src/relphot_images.c =================================================================== --- /branches/eam_branches/relphot.20210521/src/relphot_images.c (revision 41614) +++ /branches/eam_branches/relphot.20210521/src/relphot_images.c (revision 41615) @@ -144,5 +144,6 @@ freeTGroupBins (Ncatalog); - GridCorrectionSave ("test.grid.fits"); + GridCorrectionSave (); + // end of if (NLOOP > 0) block : this loop determines the offsets per chip } else {