Index: branches/eam_branches/relphot.20210521/src/GridOps.c
===================================================================
--- branches/eam_branches/relphot.20210521/src/GridOps.c	(revision 41602)
+++ branches/eam_branches/relphot.20210521/src/GridOps.c	(revision 41603)
@@ -1,898 +1,302 @@
 # include "relphot.h"
 
-enum {
-    GRID_FITTED,
-    GRID_FROZEN,
-    GRID_REFERENCE,
-};
-
-static int     Ngrid;  // number of grid elements (gridX * gridY)
-static float   *gridM; // magnitude offset for this grid cell 
-static float   *gridS; // stdev of the magnitude offset for this grid cell 
-static int     *gridN; // number of stars used to measure the magnitude offset for this grid cell 
-static int     *gridV; // data mode for this cell: fitted, frozen, reference
-static int      gridX; // number of grid elements in X direction
-static int      gridY; // number of grid elements in Y direction
-
-static int    **bin;   // link from catalog, measure to grid element
-static int    **Xmeas; // grid x-coordinate for a measurement
-static int    **Ymeas; // grid y-coordinate for a measurement
-
-static int    **clist; // link from measurement on a cell to catalog containing measurement
-static off_t  **mlist; // link from measurement on a cell to measurement in a catalog
-static off_t   *Nlist; // number of measurements for each grid cell
-static off_t   *NLIST; // allocated number of measurements for each grid cell
-
-static struct {
-  int Nchip;
-  int Mx, My;  /* mosaic size in chips */
-  int Nx, Ny;  /* chip size in pixels */
-  int *Fx, *Fy;  /* chip flip */
-  int *Ox, *Oy;  /* chip offset */
-  int *valid;
-  char **ccdname;
-} camera;    
-
-static int *ccdnum;
-static char *config;
-
-void initGrid (int dX, int dY) {
-  OHANA_UNUSED_PARAM(dX);
-  OHANA_UNUSED_PARAM(dY);
-
-  int i, N, ccdnum_max, refX, refY, refBin;
-  char *p, field[64], line[256];
-  int *Fx, *Fy;  /* chip flip */
-  int *Ox, *Oy;  /* chip offset */
-  char **ccdname;
-
-  /* load camera config file */
-  config = LoadConfigFile (CameraConfig);
-  if (config == (char *) NULL) {
-    fprintf (stderr, "ERROR: can't find camera config file %s\n", CameraConfig);
-    exit (1);
-  }
-
-  /* load basic mosaic parameters */
-  ScanConfig (config, "NCCD", "%d", 1, &camera.Nchip);
-  ScanConfig (config, "MOSAIC_X", "%d", 1, &camera.Mx);
-  ScanConfig (config, "MOSAIC_Y", "%d", 1, &camera.My);
-  ScanConfig (config, "NAXIS1", "%d", 1, &camera.Nx);
-  ScanConfig (config, "NAXIS2", "%d", 1, &camera.Ny);
-
-  ScanConfig (config, "REFCELL.X", "%d", 1, &refX);
-  ScanConfig (config, "REFCELL.Y", "%d", 1, &refY);
-
-  // temporary storage
-  ALLOCATE (Fx,      int,    camera.Nchip);
-  ALLOCATE (Fy,      int,    camera.Nchip);
-  ALLOCATE (Ox,      int,    camera.Nchip);
-  ALLOCATE (Oy,      int,    camera.Nchip);
-  ALLOCATE (ccdname, char *, camera.Nchip);
-  ALLOCATE (ccdnum,  int,    camera.Nchip);
-
-  /* load per-chip parameters */
-  for (i = 0; i < camera.Nchip; i++) {
-    ALLOCATE (ccdname[i], char, 256);
-    sprintf (field, "CCD.%d", i);
-    ScanConfig (config, field, "%s", 1, line);
-    // XXX get error status!
-
-    sscanf (line, "%s %d %d %d %d", ccdname[i], &Ox[i], &Oy[i], &Fx[i], &Fy[i]);
+/* 
+   We define a 'grid correction', essentially the flat-field correction, as a correction per 
+   photcode subdivided into an NxN array.  The dimensions of the chips corresponding to a
+   photcode will need to be added to the photcode table.  This means a schema update,
+   which I detest.  For now (2021.05.16), I will hard-wire the GPC1 / GPC2 chip size and
+   worry about HSC & Megacam in the future.
+
+   We have a collection of photcodes, with photcodeID limited by design to 64k (unsigned short).
+   
+   Thus we can generate an array of pointers to the grid correction structures and access them 
+   by photcode.
+
+
+ */
+
+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;
+
+// PS1: # define NX_CHIP 4900
+// PS1: # define NY_CHIP 4900
+// PS1: # define NX_BIN 16
+// PS1: # define NY_BIN 16
+
+# define NX_CHIP 1000
+# define NY_CHIP 1000
+# define NX_BIN 2
+# define NY_BIN 2
+
+void initGridBins () {
+
+  // allocate the full possible range of GridCorrectionType pointers, Nphotcode
+  // loop over all images to find actual existing photcodes
+  // generate initial grid values for each existing photcode
+
+  if (!GRID_ZEROPT) return; // skip if we are ignoring the grid correction
+
+  if (GridCorr) return;
+
+  PhotCodeData *photcodes = GetPhotcodeTable();
+  if (!photcodes) return;
+
+  // we have photcodes->Ncodes actually loaded.  loop over them and allocate
+  // array only as large as the max photcodes->code[i].code
+
+  int maxCode = 0;
+  for (int i = 0; i < photcodes->Ncode; i++) {
+    maxCode = MAX(maxCode, photcodes->code[i].code);
+    myAssert (maxCode < 0x10000, "oops");
+  }
+
+  NGridCorr = maxCode + 1;
+  ALLOCATE (GridCorr, GridCorrectionType *, NGridCorr);
+  for (int i = 0; i < NGridCorr; i++) {
+    GridCorr[i] = NULL;
+  }
+  return;
+}
+
+void freeGridBins() {
+
+  if (!GridCorr) return;
+
+  for (int code = 0; code < NGridCorr; code++) {
+    if (!GridCorr[code]) continue;
+
+    for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
+      FREE (GridCorr[code]-> Mgrid[ix]);
+      FREE (GridCorr[code]->dMgrid[ix]);
+      FREE (GridCorr[code]->nMgrid[ix]);
+    }      
+
+    FREE (GridCorr[code]-> Mgrid);
+    FREE (GridCorr[code]->dMgrid);
+    FREE (GridCorr[code]->nMgrid);
+
+    FREE(GridCorr[code]);
+  }
+
+  FREE (GridCorr);
+  GridCorr = NULL;
+
+  return;
+}
+
+/* for GPC1, we have 60 chips, 5 filters, 16x16 grid cells = 2MB of memory for this stuff
+   if we go to 64x64 grid cells (~75 pixels), then it is still only 29MB */
+void initGrid () {
+
+  if (!GRID_ZEROPT) return;
+
+  off_t Nimages = 0;
+  Image *images = getimages (&Nimages, NULL);
+
+  for (int i = 0; i < Nimages; i++) {
+    int code = images[0].photcode;
+    myAssert (code >= 0, "oops");
+    myAssert (code < NGridCorr, "oops");
     
-    p = ccdname[i];
-    while (!isdigit(*p) && *p) p++;
-    if (*p == 0) continue;
-    ccdnum[i] = atoi (p);
-  }
-
-  /* we now have the parameters loaded into a minimal length list; reshuffle to the a list of length 0 - MAX(ccdnum) */
-  ccdnum_max = 0;
-  for (i = 0; i < camera.Nchip; i++) {
-    ccdnum_max = MAX(ccdnum_max, ccdnum[i]);
-  }
-  ccdnum_max ++;
-
-  if (ccdnum_max < camera.Nchip) {
-    fprintf (stderr, "problem with camera config: duplicate ccd IDs\n");
-    exit (1);
-  }
-  
-  if (ccdnum_max > 0x1000) {
-    fprintf (stderr, "problem with camera config: absurd max ccd ID number %d\n", ccdnum_max);
-    exit (1);
-  }
-
-  ALLOCATE (camera.Fx, 	    int,    ccdnum_max);
-  ALLOCATE (camera.Fy, 	    int,    ccdnum_max);
-  ALLOCATE (camera.Ox, 	    int,    ccdnum_max);
-  ALLOCATE (camera.Oy, 	    int,    ccdnum_max);
-  ALLOCATE (camera.valid,   int,    ccdnum_max);
-  ALLOCATE (camera.ccdname, char *, ccdnum_max);
-
-  for (i = 0; i < ccdnum_max; i++) {
-    camera.valid[i] = FALSE;
-  }
-
-  for (i = 0; i < camera.Nchip; i++) {
-    N = ccdnum[i];
-    camera.Fx[N] = Fx[i];
-    camera.Fy[N] = Fy[i];
-    camera.Ox[N] = Ox[i];
-    camera.Oy[N] = Oy[i];
-    camera.ccdname[N] = ccdname[i];
-    camera.valid[N] = TRUE;
-  }
-
-  /* define mosaic 2d correction grid: 
-   * GRID_X is the number of grid pixels per chip in the X direction */
-  gridX = RELPHOT_GRID_X * camera.Mx;
-  gridY = RELPHOT_GRID_Y * camera.My;
-  Ngrid = gridX * gridY;
-
-  ALLOCATE (gridM, float, Ngrid);
-  ALLOCATE (gridS, float, Ngrid);
-  ALLOCATE (gridN, int,   Ngrid);
-  ALLOCATE (gridV, int,   Ngrid);
-
-  // the grid bins may have one of three possible states: fitted, frozen, reference
-  // set the initial values to indicate that the bins are frozen 
-  for (i = 0; i < Ngrid; i++) {
-    gridM[i] = 0.0;
-    gridS[i] = 0.0;
-    gridN[i] = 0;
-    gridV[i] = GRID_FROZEN;
-  }
-
-  // refBin is the index of the grid cell which is kept at 0.0 (all others are relative to this)
-  refBin = refX + refY*gridX;
-  gridV[refBin] = GRID_REFERENCE;
-}
-
-void initGridBins (Catalog *catalog, int Ncatalog) {
-
-  int i, j;
-
-  if (!USE_GRID) return;
-
-  /* define cat,meas -> grid pointers */
-  ALLOCATE (bin, int *, Ncatalog);
-  ALLOCATE (Xmeas, int *, Ncatalog);
-  ALLOCATE (Ymeas, int *, Ncatalog);
-  for (i = 0; i < Ncatalog; i++) {
-    ALLOCATE (bin[i],   int, MAX (catalog[i].Nmeasure, 1));
-    ALLOCATE (Xmeas[i], int, MAX (catalog[i].Nmeasure, 1));
-    ALLOCATE (Ymeas[i], int, MAX (catalog[i].Nmeasure, 1));
-    for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
-  }
-
-  /* define grid -> cat,meas pointers */
-  ALLOCATE (Nlist, off_t,   Ngrid);
-  ALLOCATE (NLIST, off_t,   Ngrid);
-  ALLOCATE (clist, int *,   Ngrid);
-  ALLOCATE (mlist, off_t *, Ngrid);
-
-  for (i = 0; i < Ngrid; i++) {
-    Nlist[i] = 0;
-    NLIST[i] = 100;
-    ALLOCATE (clist[i], int,   NLIST[i]);
-    ALLOCATE (mlist[i], off_t, NLIST[i]);
-  }
-}
-
-void freeGridBins (int Ncatalog) {
-
-  int i;
-
-  if (!USE_GRID) return;
-
-  /* define cat,meas -> grid pointers */
-  for (i = 0; i < Ncatalog; i++) {
-    free (bin[i]);
-    free (Xmeas[i]);
-    free (Ymeas[i]);
-  }
-  free (bin);
-  free (Xmeas);
-  free (Ymeas);
-
-  /* define grid -> cat,meas pointers */
-  for (i = 0; i < Ngrid; i++) {
-    free (clist[i]);
-    free (mlist[i]);
-  }
-  free (Nlist);
-  free (NLIST);
-  free (clist);
-  free (mlist);
-}
-
-int setGridMeasure (off_t meas, int cat, double X, double Y, int ccdnum) {
-
-  int ix, iy, Cx, Cy, i;
-  double x, y;
-
-  /* X, Y are chip coords on chip ccdnum */
-
-  /* normalize X & Y */
-  x = X;
-  if (camera.Fx[ccdnum]) x = camera.Nx - X;
-  y = Y;
-  if (camera.Fy[ccdnum]) y = camera.Ny - Y;
-
-  /* grid coords on the chip */
-  Cx = MIN (MAX ((x / camera.Nx) * RELPHOT_GRID_X, 0), RELPHOT_GRID_X - 1);
-  Cy = MIN (MAX ((y / camera.Ny) * RELPHOT_GRID_Y, 0), RELPHOT_GRID_Y - 1);
-  
-  /* coordinates in the grid */
-  ix = Cx + camera.Ox[ccdnum]*RELPHOT_GRID_X;
-  iy = Cy + camera.Oy[ccdnum]*RELPHOT_GRID_Y;
-
-  i = ix + iy*gridX;
-
-  bin[cat][meas] = i;
-  Xmeas[cat][meas] = x + camera.Ox[ccdnum]*camera.Nx;
-  Ymeas[cat][meas] = y + camera.Oy[ccdnum]*camera.Ny;
-  clist[i][Nlist[i]] = cat;
-  mlist[i][Nlist[i]] = meas;
-
-  // for the moment, add up the total grid count
-  gridN[i]++;
-
-  Nlist[i] ++;
-  if (Nlist[i] == NLIST[i]) {
-    NLIST[i] += 100;
-    REALLOCATE (clist[i], int,   NLIST[i]);
-    REALLOCATE (mlist[i], off_t, NLIST[i]);
-  }	
-  return (TRUE);
-
-  fprintf (stderr, "error: star out of grid\n");
-  exit (1);
-}
-
-int showGridCount() {
-
-    int ix, iy, i;
-
-    for (iy = 0; iy < gridY; iy++) {
-	for (ix = 0; ix < gridX; ix++) {
-	    i = ix + iy*gridX;
-	    fprintf (stderr, "%3d  ", gridN[i]);
-	}
-	fprintf (stderr, "\n");
+    // valid photcodes values (code) are in range 1 <= code < 0x10000
+    // photcode == 0 are e.g., PHU (mosaic) images, and should be ignored here
+    if (!code) continue;
+
+    if (GridCorr[code]) continue; // already created this one
+
+    ALLOCATE(GridCorr[code], GridCorrectionType, 1); 
+    GridCorr[code]->photcode = code;
+    GridCorr[code]->Nx = NX_BIN;
+    GridCorr[code]->Ny = NY_BIN;
+    GridCorr[code]->dX = NX_BIN / (float) NX_CHIP;
+    GridCorr[code]->dY = NY_BIN / (float) NX_CHIP;
+
+    // we are normally accessing this array randomly, so there is no advantage to
+    // doing Mgrid[y][x] vs Mgrid[x][y]
+    ALLOCATE (GridCorr[code]-> Mgrid, float *, NX_BIN);
+    ALLOCATE (GridCorr[code]->dMgrid, float *, NX_BIN);
+    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);
+    }      
+  }
+  resetMgrid(); // start with values of 0
+}
+
+// reset the values in the arrays to 0
+void resetMgrid () {
+
+  if (!GRID_ZEROPT) return;
+  if (!GridCorr) return;
+
+  for (int code = 0; code < NGridCorr; code++) {
+    if (!GridCorr[code]) continue;
+
+    for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
+      for (int iy = 0; iy < GridCorr[code]->Ny; iy++) {
+	GridCorr[code]-> Mgrid[ix][iy] = 0.0;
+	GridCorr[code]->dMgrid[ix][iy] = 0.0;
+	GridCorr[code]->nMgrid[ix][iy] =   0;
+      }
+    }      
+  }
+}
+
+void setMgrid (Catalog *catalog, int Ncatalog) {
+
+  // check if we are actually doing this step
+  if (!GRID_ZEROPT) return;
+  if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return;
+
+  resetMgrid(); // start with values of 0
+
+  // loop over all measurements, accumulate Sum (in Mgrid), Sum2 (in dMgrid), and Npts
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  for (int nc = 0; nc < Ncatalog; nc++) {
+    for (int na = 0; na < catalog[nc].Naverage; na++) {
+
+      int nm = catalog[nc].averageT[na].measureOffset;
+      for (int k = 0; k < catalog[nc].averageT[na].Nmeasure; k++, nm++) {
+
+	// XXX I need to skip bad measurements, but I am not certain this
+	// flag is set correctly.  see note below 
+	if (catalog[nc].measureT[nm].dbFlags & MEAS_BAD) continue;
+
+	float Mcal = getMcal  (nm, nc, MAG_CLASS_PSF);
+	if (isnan(Mcal)) continue;
+
+	float Mgrp = getMgrp  (nm, nc, catalog[nc].measureT[nm].airmass, NULL);
+	if (isnan(Mgrp)) continue;
+
+	float Mmos  = getMmos  (nm, nc);
+	if (isnan(Mmos)) continue;
+
+	// Mrel* is the average magnitude for this star.  For PS1 stacks, we have too much
+	// PSF variability.  We need to calibrate the PSF magnitudes separately from the
+	// Aperture-like magnitues.  (We have an option to use the kron magnitudes or the
+	// other apertures here).  I basically need to do this analysis separately for each
+	// magnitude type
+    
+	float MrelPSF = getMrel  (catalog, nm, nc, MAG_CLASS_PSF, MAG_SRC_CHP);
+	if (isnan(MrelPSF)) continue;
+      
+	float MsysPSF = PhotSysTiny (&catalog[nc].measureT[nm], &catalog[nc].averageT[na], &catalog[nc].secfilt[na*Nsecfilt], MAG_CLASS_PSF);
+	if (isnan(MsysPSF)) continue;
+
+	// what about Mflat?
+	float Moff =  Mcal + Mgrp + Mmos;
+
+	// Msys = Mrel + Moff + Mgrid
+	// thus Mgrid = Msys - Mrel - Moff
+
+	int code = catalog[nc].measureT[nm].photcode;
+	if (code <= 0) continue;
+	if (code >= NGridCorr) continue; // does not match one of our image, skip
+	
+	GridCorrectionType *grid = GridCorr[code];
+	if (!grid) continue; // does not match one of our images, skip
+
+	// edge effects could cause some positions to be slightly out of range
+	// probably should trap extreme outliers
+	int ix = MIN(MAX(0, (int)(catalog[nc].measureT[nm].Xccd * grid->dX)), grid->Nx - 1);
+	int iy = MIN(MAX(0, (int)(catalog[nc].measureT[nm].Yccd * grid->dY)), grid->Ny - 1);
+
+	float dM = MsysPSF - MrelPSF - Moff;
+	
+	// XXX by the time we get here, we should have already mostly fixed up the zero
+	// points.  reject measurements which are way off
+	if (fabs(dM) > 0.5) continue;
+
+	grid-> Mgrid[ix][iy] += dM;
+	grid->dMgrid[ix][iy] += dM*dM;
+	grid->nMgrid[ix][iy] ++;
+      }
     }
-    return (TRUE);
-}
-
-float getMgrid (off_t meas, int cat) {
-
-  int i;
-  float value;
-
-  if (!USE_GRID) return (0);
-  i = bin[cat][meas];
-  if (i == -1) return (NAN);
-
-  // during the grid annealing process, we skip over grid cells until they have enough
-  // valid stars to be fitted
-  if (gridV[i] == GRID_FROZEN) {
-    return (NAN);
-  }
-
-  value = gridM[i];
-  return (value);
-}
-
-/* direct (non-iterative) solution for Mgrid values for all grid bins */
-void setMgridDirect (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
-
-  int **gotstar, **gridmeas;
-  int i, j, k, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys, Ngrp;
-  double **A, **B, *Mjx, *Wjx;
-  float Msys, Mcal, Mmos, Mgrp, Merr, Wsys;
-  double Mj, Wj;
-  
-  if (!USE_GRID) return;
-
-  ALLOCATE (A, double *, Ngrid);
-  ALLOCATE (B, double *, Ngrid);
-  for (i = 0; i < Ngrid; i++) {
-    ALLOCATE (A[i], double, Ngrid);
-    ALLOCATE (B[i], double, 1);
-    memset (A[i], 0, Ngrid*sizeof(double));
-    memset (B[i], 0, sizeof(double));
-  }
-
-  Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = Ngrp = 0;
-
-  ALLOCATE (gotstar, int *, Ncatalog);
-  for (i = 0; i < Ncatalog; i++) {
-    ALLOCATE (gotstar[i], int, catalog[i].Naverage);
-  }  
-
-  // set up gridmeas table : grid index for each measurement
-  ALLOCATE (gridmeas, int *, Ncatalog);
-  for (i = 0; i < Ncatalog; i++) {
-    ALLOCATE (gridmeas[i], int, catalog[i].Nmeasure);
-    for (j = 0; j < catalog[i].Nmeasure; j++) {
-      gridmeas[i][j] = -1;
-    }
-  }  
-  for (i = 0; i < Ngrid; i++) {
-    for (j = 0; j < Nlist[i]; j++) {
-      int m, c;
-      m = mlist[i][j];
-      c = clist[i][j];
-      gridmeas[c][m] = i;
-    }
-  }
-
-  // as we loop over the grid cells, we need to accumulate the values of Wjx for each star
-  ALLOCATE (Wjx, double, Ngrid);
-  ALLOCATE (Mjx, double, Ngrid);
-
-  // accumulate the elements of the matrix equation.  We have an equation of the form: Ax = B
-  // where x is the vector of grid cell values G_x (x = 0 - Ngrid), A is an Ngrid x Ngrid matrix,
-  // and B is an Ngrid vector.  For a cell A(x,y), we need the following elements from each
-  // star which touches grid cell (x):
-  // 
-  // Mj  : sum over all measurements of Msys / dMsys^2 
-  // Wj  : sum over all measurements of 1.0 / dMsys^2 
-  // Mjx : sum over all measurements which touch cell (x) of Msys / dMsys^2
-  // Wjx : sum over all measurements which touch cell (x) of 1.0 / dMsys^2
-
-  // Mjx and Wjx can be calculated by summing over all measurements which touch the cell 
-  // Mj requires looping over stars which touch (x)
-
-  // this is tricky because we need to know both the measurements which touch a cell
-  // and the stars for which any measurement touches a cell.  we need to accumulate
-  // sums for each star which touches as cell on both bases.
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-  int thisCode = photcodes[0][0].code;
-  int Nsec = GetPhotcodeNsec(thisCode);
-
-  for (i = 0; i < Ngrid; i++) {
-    
-    for (j = 0; j < Ncatalog; j++) {
-      memset (gotstar[j], 0, catalog[j].Naverage*sizeof(int));
-    }  
-
-    // we are looping over the stars, but doing so by looping over the set of measurements:
-    // every star which touches this grid cell has a measurement in Nlist[i]
-    for (j = 0; j < Nlist[i]; j++) {
-      
-      int mx, c, n, m0, Npts;
-
-      mx = mlist[i][j];
-      c  = clist[i][j];
-      n  = catalog[c].measureT[mx].averef;
-      
-      // if we have already visited this star, skip the stuff below
-      if (gotstar[c][n]) continue;
-      gotstar[c][n] = TRUE;
-
-      // skip stars marked as BAD
-      if (catalog[c].secfilt[n*Nsecfilt+Nsec].flags & STAR_BAD) {
-	Nrel ++;
-	continue;
-      }
-
-      m0 = catalog[c].average[n].measureOffset;
-
-      // we accumuate an entry for each cell
-      memset (Wjx, 0, Ngrid*sizeof(double));
-      memset (Mjx, 0, Ngrid*sizeof(double));
-      Npts = Mj = Wj = 0.0;
-
-      // if we have not yet visited this star, accumulate the Mj, Wj entries
-      for (k = 0; k < catalog[c].average[n].Nmeasure; k++) {
-
-	int m, Ng;
-
-	m = m0 + k;
-
-	// skip measurements marked as BAD
-	if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
-	  Nbad ++;
+  }
+
+  // now calculate Mgrid, dMgrid, nMgrid from Sum, Sum, Npt
+  for (int code = 0; code < NGridCorr; code++) {
+    if (!GridCorr[code]) continue;
+
+    for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
+      for (int iy = 0; iy < GridCorr[code]->Ny; iy++) {
+
+	// cells without sufficient coverage stay at 0.0
+	if (GridCorr[code]->nMgrid[ix][iy] < 5) {
+	  GridCorr[code]-> Mgrid[ix][iy] = 0.0;
+	  GridCorr[code]->dMgrid[ix][iy] = 0.0;
 	  continue;
 	}
 
-	// skip images marked as BAD
-	Mcal = getMcal  (m, c, MAG_CLASS_PSF);
-	if (isnan(Mcal)) {
-	  Ncal ++;
-	  continue;
-	}
-
-	// skip mosaics marked as BAD
-	Mmos = getMmos  (m, c);
-	if (isnan(Mmos)) {
-	  Nmos ++;
-	  continue;
-	}
-
-	// skip mosaics marked as BAD
-	Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
-	if (isnan(Mgrp)) {
-	  Ngrp ++;
-	  continue;
-	}
-
-	// select the color- and airmass-corrected observed magnitude for this star
-	// XXX need to be able to turn off the color-correction until initial average mags are found
-	Msys = PhotCatTiny (&catalog[c].measureT[m], MAG_CLASS_PSF);
-	if (isnan(Msys)) {
-	  Nsys++;
-	  continue;
-	}
-
-	// mag-error for this measurement
-	Merr =  MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-
-	// disable Wsys for now
-	Wsys = TRUE ? 1.0 : 1.0 / SQ(Merr);
-
-	Ng = gridmeas[c][m];
-	if (Ng == -1) continue;  // skip measurements which do not touch any cell
-
-	Mj += Msys * Wsys;  // we are only including measurements touching this cell
-	Wj += Wsys;  // we are only including measurements touching this cell
-	Npts ++;
-	Ngood ++;
-
-	Mjx[Ng] += Msys * Wsys;  // we are only including measurements touching cell (x)
-	Wjx[Ng] += Wsys;  // we are only including measurements touching cell (x)
-      }
-
-      // some stars will not have any valid measurements, skip these
-      if (Npts == 0) continue;
-
-      B[i][0] += Mj*Wjx[i]/Wj - Mjx[i];
-      A[i][i] -= Wjx[i];
-      for (k = 0; k < Ngrid; k++) {
-	A[i][k] += Wjx[i]*Wjx[k]/Wj;
-	// fprintf (stderr, "%3.0f ", Wjx[k]);
+	float Mgrid  = GridCorr[code]-> Mgrid[ix][iy] / GridCorr[code]->nMgrid[ix][iy]; // average Mgrid
+	float Mgrid2 = GridCorr[code]-> dMgrid[ix][iy] / GridCorr[code]->nMgrid[ix][iy]; // average Mgrid^2
+
+	float r = GridCorr[code]->nMgrid[ix][iy] / (float) (GridCorr[code]->nMgrid[ix][iy] - 1.0); // pop -> sample stdev
+
+	GridCorr[code]-> Mgrid[ix][iy] = Mgrid;
+	GridCorr[code]->dMgrid[ix][iy] = sqrt(r*(Mgrid2 - Mgrid*Mgrid)); // sample stdev
+	fprintf (stderr, "grid code %d, %d x %d : %f +/- %f : %d\n", code, ix, iy,
+		 GridCorr[code]-> Mgrid[ix][iy], GridCorr[code]->dMgrid[ix][iy], GridCorr[code]->nMgrid[ix][iy]);
       }
     }
   }
-
-  if (1) {
-
-    FILE *f;
-    Header theader;
-    Matrix matrix;
-
-    /* we are writing to this file */
-    f = fopen ("matrix.fits", "w");
-    if (f == (FILE *) NULL) { 
-      fprintf (stderr, "cannot open matrix.fits for output\n");
-      return;
-    }
-
-    /* save grid mag values */
-    gfits_init_header (&theader);
-    theader.Naxes = 2;
-    theader.Naxis[0] = Ngrid;
-    theader.Naxis[1] = Ngrid;
-    theader.bitpix   = -32;
-    gfits_create_Theader (&theader, "IMAGE");
-    gfits_modify (&theader, "EXTNAME", "%s", 1, "MATRIX");
-    gfits_create_matrix  (&theader, &matrix);
-    for (i = 0; i < Ngrid; i++) {
-      for (j = 0; j < Ngrid; j++) {
-	gfits_set_matrix_value (&matrix, i, j, (double) A[i][j]);
-      }
-    }
-    gfits_fwrite_header (f, &theader);
-    gfits_fwrite_matrix (f, &matrix);
-    gfits_free_matrix (&matrix);
-
-    gfits_modify (&theader, "EXTNAME", "%s", 1, "TRPOSE");
-    gfits_create_matrix  (&theader, &matrix);
-    for (i = 0; i < Ngrid; i++) {
-      for (j = 0; j < Ngrid; j++) {
-	gfits_set_matrix_value (&matrix, i, j, (double) A[j][i]);
-      }
-    }
-    gfits_fwrite_header (f, &theader);
-    gfits_fwrite_matrix (f, &matrix);
-    gfits_free_matrix (&matrix);
-    fclose (f);
-  }
-
-  dgaussjordan (A, B, Ngrid, 1);
-
-  fprintf (stderr, "grid cells fitted (Ngood: %d, Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Ngood, Nbad, Nmos, Ncal, Nrel, Nsys);
-
-  for (i = 0; i < Ngrid; i++) {
-    gridM[i] = B[i][0];
-    gridS[i] = sqrt(A[i][i]);
-    gridN[i] = Ngood;
-  }
-
-  free (Wjx);
-  free (Mjx);
-
-  for (i = 0; i < Ngrid; i++) {
-    free (A[i]);
-  }
-  free (A);
-  free (B);
-
-  for (i = 0; i < Ncatalog; i++) {
-    free (gotstar[i]);
-    free (gridmeas[i]);
-  }  
-  free (gotstar);
-  free (gridmeas);
-}
-
-/* determine Mgrid values for all grid bins */
-void setMgrid (Catalog *catalog, FlatCorrectionTable *flatcorr) {
-
-  int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ngrp, Ncal, Nrel, Nsys, Nfit;
-  double *list, *dlist;
-  float Msys, Mrel, Mcal, Mmos, Mgrp;
-
-  StatType stats;
-  liststats_setmode (&stats, "INNER_WTMEAN");
+  return;
+}
+
+float getMgrid (Measure *measure) {
+
+  if (!GRID_ZEROPT) return 0.0;
+  if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return 0.0;
+
+  int code = measure->photcode;
+  if (code <= 0) return 0.0;
+  if (code >= NGridCorr) return 0.0; // does not match one of our image, skip
+	
+  GridCorrectionType *grid = GridCorr[code];
+  if (!grid) return 0.0; // does not match one of our images, skip
   
-  if (!USE_GRID) return;
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  Nmax = Nlist[0];
-  for (i = 0; i < Ngrid; i++) {
-    Nmax = MAX (Nmax, Nlist[i]);
-  }
-  ALLOCATE (list, double, Nmax);
-  ALLOCATE (dlist, double, Nmax);
-
-  Nbad = Ncal = Nmos = Ngrp = Nrel = Nsys = Nfit = 0;
-
-  for (i = 0; i < Ngrid; i++) {
-
-    N = 0;
-    for (j = 0; j < Nlist[i]; j++) {
-      
-      m = mlist[i][j];
-      c = clist[i][j];
-      
-      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
-	Nbad ++;
-	continue;
-      }
-      Mcal = getMcal  (m, c, MAG_CLASS_PSF);
-      if (isnan(Mcal)) {
-	Ncal ++;
-	continue;
-      }
-      Mmos = getMmos  (m, c);
-      if (isnan(Mmos)) {
-	Nmos ++;
-	continue;
-      }
-      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
-      if (isnan(Mgrp)) {
-	Ngrp ++;
-	continue;
-      }
-      Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
-      if (isnan(Mrel)) {
-	Nrel ++;
-	continue;
-      }
-      
-      n = catalog[c].measureT[m].averef;
-      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
-      if (isnan(Msys)) {
-	Nsys++;
-	continue;
-      }
-      list[N] = Msys - Mrel - Mcal - Mmos - Mgrp;
-      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-      N++;
-    }
-
-    // the reference Cell is forced to have a value of 0.0, and is never changed to GRID_FITTED
-    if (gridV[i] == GRID_REFERENCE) {
-      gridM[i] = 0.0;
-      gridS[i] = 0.0;
-      gridN[i] = N;
-      continue;
-    }
-
-    // until we have enough valid measurements on this grid cell, skip it
-    if (N < GRID_TOOFEW) {
-      gridV[i] = GRID_FROZEN;
-      continue;
-    }
-
-    liststats (list, dlist, NULL, N, &stats);
-    gridM[i] = stats.mean;
-    gridS[i] = stats.sigma;
-    gridN[i] = N;
-    gridV[i] = GRID_FITTED;
-    Nfit++;
-  }
-
-  fprintf (stderr, "%d of %d grid cells fitted (+ reference cell) (Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Nfit, Ngrid, Nbad, Nmos, Ncal, Nrel, Nsys);
-
-  free (list);
-  free (dlist);
-}
-
-void plot_grid (Catalog *catalog, FlatCorrectionTable *flatcorr) {
-
-  int i, j, m, c, n, N, Narea;
-  float Msys, Mrel, Mcal, Mmos, Mgrp;
-  double *xlist, *Mlist, *dlist, *ylist;
-  Graphdata graphdata;
-
-  if (!USE_GRID) return;
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  N = 0;
-  for (i = 0; i < Ngrid; i++) 
-    N += Nlist[i];
-
-  ALLOCATE (xlist, double, N);
-  ALLOCATE (ylist, double, N);
-  ALLOCATE (Mlist, double, N);
-  ALLOCATE (dlist, double, N);
-
-  Narea = 0;
-  N = 0;
-  for (i = 0; i < Ngrid; i++) {
-    for (j = 0; j < Nlist[i]; j++) {
-      
-      m = mlist[i][j];
-      c = clist[i][j];
-      
-      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
-	Narea ++;
-	continue;
-      }
-      Mcal  = getMcal  (m, c, MAG_CLASS_PSF);
-      if (isnan(Mcal)) continue;
-      Mmos  = getMmos  (m, c);
-      if (isnan(Mmos)) continue;
-      Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
-      if (isnan(Mgrp)) continue;
-      Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
-      if (isnan(Mrel)) continue;
-
-      n = catalog[c].measureT[m].averef;
-      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
-
-      xlist[N] = Xmeas[c][m];
-      ylist[N] = Ymeas[c][m];
-      Mlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp;
-      dlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp - gridM[i];
-      N++;
-    }
-  }
-
-  fprintf (stderr, "skipped %d meas for area\n", Narea);
-
-  plot_defaults (&graphdata);
-  graphdata.ymin = PlotdMmin;
-  graphdata.ymax = PlotdMmax;
-  plot_list (&graphdata, xlist, Mlist, N, "X vs dM raw", "%s.XdM.png", OUTROOT);
-  plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "%s.XdMf.png", OUTROOT);
-  plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "%s.YdMf.png", OUTROOT);
-
-  plot_defaults (&graphdata);
-  plot_list (&graphdata, xlist, ylist, N, "X vs Y", "%s.XY.png", OUTROOT);
-
-  free (ylist);
-  free (xlist);
-  free (Mlist);
-  free (dlist);
-
-}
-
-void dump_grid () { 
-
-  off_t i, Nimage;
-  int j, Nbytes, Nformat;
-  FILE *f;
-  Header header, theader;
-  Matrix matrix;
-  Mosaic *refmosaic;
-  char *filename;
-  char formatline[32], key[32], value[64];
-    
-  Nbytes = strlen (OUTROOT) + 6;
-  ALLOCATE (filename, char, Nbytes);
-  snprintf (filename, Nbytes, "%s.fits", OUTROOT);
-
-  /* select reference mosaic image */
-  // off_t *imlist = SelectRefMosaic (&refmosaic, &Nimage); return value ignored
-  SelectRefMosaic (&refmosaic, &Nimage);
-
-  /* we are writing to this file */
-  f = fopen (filename, "w");
-  if (f == (FILE *) NULL) { 
-    fprintf (stderr, "cannot open %s for output\n", filename);
-    free (filename);
-    return;
-  }
-
-  /* create empty phu */
-  gfits_init_header (&header);
-  header.extend = TRUE;
-  gfits_create_header (&header);
-  gfits_create_matrix (&header, &matrix);
-  gfits_modify (&header, "NEXTEND", OFF_T_FMT, 1,  Nimage + 3);
-  gfits_modify (&header, "FILTER", "%s", 1, photcodes[0][0].name);  // XXXX note that this expects a single photcode, enforced in initialize.d
-  gfits_modify_alt (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis");
-
-  // we need to add lines to the PHU to identify the camera and format; these are used by the ipp config system
-  // Note that config must have been loaded (and not freed) above.
-  ScanConfig (config, "NFORMAT", "%d", 1, &Nformat);
-  for (i = 1; i <= Nformat; i++) {
-      ScanConfig (config, "FORMAT", "%s", i, formatline);
-      sscanf (formatline, "%s %s", key, value);
-      gfits_modify (&header, key, "%s", 1, value);
-  }
-
-  gfits_fwrite_header (f, &header);
-  gfits_fwrite_matrix (f, &matrix);
-  gfits_free_matrix (&matrix);
-
-  /* save grid mag values */
-  gfits_init_header (&theader);
-  theader.Naxes = 2;
-  theader.Naxis[0] = gridX;
-  theader.Naxis[1] = gridY;
-  theader.bitpix   = -32;
-  gfits_create_Theader (&theader, "IMAGE");
-  gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
-  gfits_modify (&theader, "EXTNAME", "%s", 1, "MAG_OFFSET");
-  gfits_create_matrix  (&theader, &matrix);
-  for (i = 0; i < gridX; i++) {
-    for (j = 0; j < gridY; j++) {
-      gfits_set_matrix_value (&matrix, i, j, (double) gridM[i + j*gridX]);
-    }
-  }
-  write_coords (&theader, &refmosaic[0].coords);
-  gfits_fwrite_header (f, &theader);
-  gfits_fwrite_matrix (f, &matrix);
-  gfits_free_matrix (&matrix);
-
-  /* save grid Nmeas values */
-  gfits_modify (&theader, "EXTNAME", "%s", 1, "NMEAS");
-  gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
-  gfits_create_matrix  (&theader, &matrix);
-  for (i = 0; i < gridX; i++) {
-    for (j = 0; j < gridY; j++) {
-      gfits_set_matrix_value (&matrix, i, j, (double) gridN[i + j*gridX]);
-    }
-  }
-  write_coords (&theader, &refmosaic[0].coords);
-  gfits_fwrite_header (f, &theader);
-  gfits_fwrite_matrix (f, &matrix);
-  gfits_free_matrix (&matrix);
-
-  /* save grid sigma values */
-  gfits_modify (&theader, "EXTNAME", "%s", 1, "SIGMA");
-  gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
-  gfits_create_matrix  (&theader, &matrix);
-  for (i = 0; i < gridX; i++) {
-    for (j = 0; j < gridY; j++) {
-      gfits_set_matrix_value (&matrix, i, j, (double) gridS[i + j*gridX]);
-    }
-  }
-  write_coords (&theader, &refmosaic[0].coords);
-  gfits_fwrite_header (f, &theader);
-  gfits_fwrite_matrix (f, &matrix);
-  gfits_free_matrix (&matrix);
-
-  /* calculate value for each CCD pixel, write out CCD images */
-  /* grid pixels are tied to detector pixels, but are flipped to match focal plane */
-  for (i = 0; i < camera.Nchip; i++) {
-    int N, ix, iy, x, y, X, Y, bin;
-
-    N = ccdnum[i];
-
-    gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[N]);
-    gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
-    gfits_modify (&theader, "NX", "%d", 1, camera.Nx);
-    gfits_modify (&theader, "NY", "%d", 1, camera.Ny);
-      
-    theader.Naxis[0] = RELPHOT_GRID_X;
-    theader.Naxis[1] = RELPHOT_GRID_Y;
-    gfits_modify (&theader, "NAXIS1", "%d", 1, RELPHOT_GRID_X);
-    gfits_modify (&theader, "NAXIS2", "%d", 1, RELPHOT_GRID_Y);
-    gfits_create_matrix  (&theader, &matrix);
-
-    for (Y = 0; Y < RELPHOT_GRID_Y; Y++) {
-      for (X = 0; X < RELPHOT_GRID_X; X++) {
-	      
-	/* normalize X & Y */
-	x = X;
-	if (camera.Fx[N]) x = RELPHOT_GRID_X - X - 1;
-	y = Y;
-	if (camera.Fy[N]) y = RELPHOT_GRID_Y - Y - 1;
-	      
-	/* coordinates in the grid */
-	ix = x + camera.Ox[N]*RELPHOT_GRID_X;
-	iy = y + camera.Oy[N]*RELPHOT_GRID_Y;
-	      
-	bin = ix + iy*gridX;
-	gfits_set_matrix_value (&matrix, X, Y, (double) gridM[bin]);
-      }
-    }
-    gfits_fwrite_header (f, &theader);
-    gfits_fwrite_matrix (f, &matrix);
-    gfits_free_matrix (&matrix);
-  }
-
-  free (filename);
-}
-
-void InterpolateGrid (float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords) {
-
-  int i, j;
-  double x, y, r, d, X, Y, dx, dy;
-  double V00, V01, V10, V11;
-  double wV00, wV01, wV10, wV11;
-  double dV00, dV01, dV10, dV11;
-  double v1, v2, value;
-  int ix, iy, N;
-
-  for (i = 0; i < Nx; i++) {
-    for (j = 0; j < Ny; j++) {
-      x = i * RELPHOT_GRID_BINNING / 2;
-      y = j * RELPHOT_GRID_BINNING / 2;
-      XY_to_RD (&r, &d, x, y, ccd);
-      RD_to_XY (&X, &Y, r, d, gcoords);
-
-      X = X / RELPHOT_GRID_BINNING;
-      Y = Y / RELPHOT_GRID_BINNING;
-
-      ix = (int) X;
-      dx = X - ix;
-      iy = (int) Y;
-      dy = Y - iy;
-
-      if (ix < 0) continue;
-      if (iy < 0) continue;
-      if (ix >= gridX) continue;
-      if (iy >= gridY) continue;
-
-      N = ix + iy*gridX;
-      V00 = gridM[N];
-      V10 = gridM[N + 1];
-      V01 = gridM[N + gridX];
-      V11 = gridM[N + gridX + 1];
-
-      dV00 = gridS[N];
-      dV10 = gridS[N + 1];
-      dV01 = gridS[N + gridX];
-      dV11 = gridS[N + gridX + 1];
-
-      wV00 = (dV00 == 0) ? 0.0 : 1 / SQ(dV00);
-      wV01 = (dV01 == 0) ? 0.0 : 1 / SQ(dV01);
-      wV10 = (dV10 == 0) ? 0.0 : 1 / SQ(dV10);
-      wV11 = (dV11 == 0) ? 0.0 : 1 / SQ(dV11);
-
-      v1 = wV00*V00*(1 + dx*dy - dx - dy) +
-	wV10*V10*(dx - dx*dy) +
-	wV01*V01*(dy - dx*dy) +
-	wV11*V11*(dx*dy);
-
-      v2 = wV00*(1 + dx*dy - dx - dy) +
-	wV10*(dx - dx*dy) +
-	wV01*(dy - dx*dy) +
-	wV11*(dx*dy);
-
-      value = v1 / v2;
-      buffer[j*Nx + i] = value;
-    }
-  }
-}
+  // edge effects could cause some positions to be slightly out of range
+  // probably should trap extreme outliers
+  int ix = MIN(MAX(0, (int)(measure->Xccd * grid->dX)), grid->Nx - 1);
+  int iy = MIN(MAX(0, (int)(measure->Yccd * grid->dY)), grid->Ny - 1);
+
+  float Mgrid = grid-> Mgrid[ix][iy];
+  return Mgrid;
+}
+
+float getMgridTiny (MeasureTiny *measure) {
+
+  if (!GRID_ZEROPT) return 0.0;
+  if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return 0.0;
+
+  int code = measure->photcode;
+  if (code <= 0) return 0.0;
+  if (code >= NGridCorr) return 0.0; // does not match one of our image, skip
+	
+  GridCorrectionType *grid = GridCorr[code];
+  if (!grid) return 0.0; // does not match one of our images, skip
+  
+  // edge effects could cause some positions to be slightly out of range
+  // probably should trap extreme outliers
+  int ix = MIN(MAX(0, (int)(measure->Xccd * grid->dX)), grid->Nx - 1);
+  int iy = MIN(MAX(0, (int)(measure->Yccd * grid->dY)), grid->Ny - 1);
+
+  float Mgrid = grid-> Mgrid[ix][iy];
+  return Mgrid;
+}