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 41610)
@@ -0,0 +1,162 @@
+# include "relphot.h"
+
+int GridCorrectionSave (char *filename) {
+
+  /* open file for input */
+  FILE *f = fopen (filename, "w");
+  if (!f) {
+    gprint (GP_ERR, "can't open file for write : %s\n", filename);
+    return FALSE;
+  }
+
+  // generate a PHU 
+
+  Header phu_header;
+  Matrix phu_matrix;
+
+  gfits_init_header   (&phu_header);
+  gfits_create_header (&phu_header);
+  gfits_modify_alt    (&phu_header, "EXTEND", "%t",  1, TRUE);
+  gfits_create_matrix (&phu_header, &phu_matrix);
+  gfits_fwrite_header (f, &phu_header);
+  gfits_fwrite_matrix (f, &phu_matrix);
+  gfits_free_header (&phu_header);
+  gfits_free_matrix (&phu_matrix);
+
+  int Nlast = -1;
+
+  while (TRUE) {
+    // getGridCorr() takes a photcode and returns the *next* valid grid correction pointer (start at -1)
+    GridCorrectionType *GridCorr = getGridCorr (&Nlast);
+    if (!GridCorr) break;
+
+    int filter = GetPhotcodeEquivCodebyCode (GridCorr->photcode);
+
+    Header header;
+    Matrix matrix;
+
+    gfits_init_header (&header);
+    header.bitpix   = -32;
+    header.Naxes    =   2;
+    header.Naxis[0] = GridCorr->Nx;
+    header.Naxis[1] = GridCorr->Ny;
+
+    /* create the appropriate header and matrix */
+    gfits_create_header (&header);
+    gfits_create_matrix (&header, &matrix);
+
+    // 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, "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");
+    if (!gfits_modify (&header, "Y_CHIP",   "%d", 1, -1))                 myAbort ("failed to set Y_CHIP");
+    if (!gfits_modify (&header, "DX_CHIP",  "%f", 1, GridCorr->dX))       myAbort ("failed to set DX_CHIP");
+    if (!gfits_modify (&header, "DY_CHIP",  "%f", 1, GridCorr->dY))       myAbort ("failed to set DY_CHIP");
+    if (!gfits_modify (&header, "NX_CHIP",  "%d", 1, GridCorr->Nx))       myAbort ("failed to set NX_CHIP");
+    if (!gfits_modify (&header, "NY_CHIP",  "%d", 1, GridCorr->Ny))       myAbort ("failed to set NY_CHIP");
+
+    int Nx = GridCorr->Nx;
+    float *value = (float *) 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];
+      }
+    }
+    
+    gfits_primary_to_extended (&header, "IMAGE", "Image Extension");
+    gfits_fwrite_header  (f, &header);
+    gfits_fwrite_matrix  (f, &matrix);
+
+    gfits_free_header (&header);
+    gfits_free_matrix (&matrix);
+  }
+  fclose (f);
+
+  return TRUE;
+}
+
+# if (0) 
+CamPhotomCorrection *CamPhotomCorrectionLoad (char *filename) {
+
+  int i, ix, iy, season, filter;
+
+  Header *phu = gfits_alloc_header();
+
+  /* open file for input */
+  FILE *f = fopen (filename, "r");
+  if (f == (FILE *) NULL) {
+    gprint (GP_ERR, "can't open file for read : %s\n", filename);
+    return FALSE;
+  }
+
+  // read the PHU header
+  if (!gfits_load_header (f, phu)) return FALSE;
+  
+  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;
+
+  CamPhotomCorrection *cam = CamPhotomCorrectionAlloc (Nx, Ny, Nfilter, Nseason);
+  cam->phu = phu;
+
+  cam->dX      = dX;  	   // superpixel sampling, x-dir
+  cam->dY      = dY;  	   // superpixel sampling, y-dir
+
+  cam->NxCCD   = NxCCD;    // pixels per chip
+  cam->NyCCD   = NyCCD;    // pixels per chip
+
+  for (i = 0; i < Nseason; i++) {
+    char field[256];
+    double value;
+    
+    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);
+  }
+
+  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;
+}
+# endif
+
Index: /branches/eam_branches/relphot.20210521/src/GridOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/GridOps.c	(revision 41609)
+++ /branches/eam_branches/relphot.20210521/src/GridOps.c	(revision 41610)
@@ -16,16 +16,4 @@
  */
 
-typedef struct {
-  unsigned short photcode;
-  float **Mgrid; // grid of average corrections
-  float **dMgrid; // grid of stdev of corrections
-    int **nMgrid; // grid of number of stars for corrections
-  int Nx; // bin = int(Xchip * (Nx / NxChip))
-  int Ny; // bin = int(Ychip * (Ny / NyChip))
-  float dX; // bin = int(Xchip * dX), dX = Nx / NxChip
-  float dY; // bin = int(Ychip * dY), dY = Ny / NyChip
-  // NxChip, NyChip = 4900,4900 for now
-} GridCorrectionType;
-
 static GridCorrectionType **GridCorr = NULL;
 static int                 NGridCorr = 0;
@@ -95,4 +83,21 @@
 
   return;
+}
+
+GridCorrectionType *getGridCorr (int *Nlast) {
+
+  if (GridCorr == NULL) return NULL;
+
+  if (*Nlast >= NGridCorr) return NULL;
+
+  if (*Nlast < 0) *Nlast = -1;
+
+  GridCorrectionType *result = NULL;
+  for (int i = *Nlast + 1; i < NGridCorr; i++) {
+    if (GridCorr[i] == NULL) continue;
+    result = GridCorr[i];
+    *Nlast = i;
+  }
+  return result;
 }
 
Index: /branches/eam_branches/relphot.20210521/src/relphot_images.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/relphot_images.c	(revision 41609)
+++ /branches/eam_branches/relphot.20210521/src/relphot_images.c	(revision 41610)
@@ -143,4 +143,6 @@
     freeMosaicBins (Ncatalog, TRUE);
     freeTGroupBins (Ncatalog);
+
+    GridCorrectionSave ("test.grid.fits");
     // end of if (NLOOP > 0) block : this loop determines the offsets per chip 
   } else {