Index: /branches/eam_branches/relphot.20210521/Makefile
===================================================================
--- /branches/eam_branches/relphot.20210521/Makefile	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/Makefile	(revision 41603)
@@ -29,5 +29,5 @@
 $(SRC)/StatDataSetOps.$(ARCH).o	 \
 $(SRC)/ConfigInit.$(ARCH).o	 \
-$(SRC)/GridOps.$(ARCH).o	 \
+$(SRC)/GridOpsNew.$(ARCH).o	 \
 $(SRC)/ImageOps.$(ARCH).o	 \
 $(SRC)/ImageSubset.$(ARCH).o	 \
@@ -92,5 +92,5 @@
 $(SRC)/StatDataSetOps.$(ARCH).o	 \
 $(SRC)/ConfigInit.$(ARCH).o	 \
-$(SRC)/GridOps.$(ARCH).o	 \
+$(SRC)/GridOpsNew.$(ARCH).o	 \
 $(SRC)/ImageOps.$(ARCH).o	 \
 $(SRC)/ImageSubset.$(ARCH).o	 \
Index: /branches/eam_branches/relphot.20210521/include/relphot.h
===================================================================
--- /branches/eam_branches/relphot.20210521/include/relphot.h	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/include/relphot.h	(revision 41603)
@@ -4,8 +4,4 @@
 # include <signal.h>
 # include <pthread.h>
-
-/* # define GRID_V1 */
-# define GRID_V2
-# define NO_IMAGE -100
 
 # define ID_SECF_STACK_PRIMARY 0x00004000
@@ -99,5 +95,5 @@
   char skipCal;		      // if TRUE, this mosaic is incomplete and should not be calibrated
   char inTGroup;
-  Coords coords;
+  Coords coords;	      // this is only used to set the mosaic center for assign_images used by region hosts
 } Mosaic; 
 
@@ -413,4 +409,5 @@
 int    MOSAIC_ZEROPT;
 int    TGROUP_ZEROPT;
+int    GRID_ZEROPT;
 
 int    FREEZE_IMAGES;
@@ -424,5 +421,4 @@
 int    MANUAL_ITERATION;
 int    NLOOP;
-int    NGRID;
 int    RESET;
 int    RESET_ZEROPTS;
@@ -444,5 +440,5 @@
 int    IMAGE_BAD;
 int    CALIBRATE_STACKS_AND_WARPS;
-int    USE_GRID;
+
 int    KEEP_UBERCAL;
 char  *OUTROOT;
@@ -457,8 +453,4 @@
 char  *PhotcodeList;
 
-int    RELPHOT_GRID_X;
-int    RELPHOT_GRID_Y;
-int    RELPHOT_GRID_BINNING;
-
 int      *photseclist;
 int      Nphotcodes;
@@ -489,4 +481,7 @@
 SkyRegion UserPatch;
 char *UserCatalog;
+
+enum {GRID_ZPT_MODE_NONE, GRID_ZPT_MODE_ALL, };
+int GRID_ZPT_MODE;
 
 enum {TGROUP_ZPT_MODE_NONE, TGROUP_ZPT_MODE_GOOD_NIGHT, TGROUP_ZPT_MODE_ALL, };
@@ -549,5 +544,5 @@
 int SetSignals (void);
 
-void          freeGridBins        PROTO((int Ncatalog));
+void          freeGridBins        PROTO((void));
 void          freeImageBins       PROTO((int Ncatalog, int doImageList));
 void          freeMosaicBins      PROTO((int Ncatalog, int doMosaicList));
@@ -561,5 +556,6 @@
 float         getMcal             PROTO((off_t meas, int cat, dvoMagClassType class));
 float         getMflat            PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog));
-float         getMgrid            PROTO((off_t meas, int cat));
+float         getMgridTiny        PROTO((MeasureTiny *measure));
+float         getMgrid            PROTO((Measure *measure));
 float         getMmos             PROTO((off_t meas, int cat));
 float         getMgrp             PROTO((off_t meas, int cat, float airmass, float *dZpt));
@@ -571,11 +567,11 @@
 ImageSubset  *getimages_subset    PROTO((off_t *N));
 void          global_stats        PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr, int nloop));
-void          initGrid            PROTO((int dX, int dY));
-void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
+void          initGrid            PROTO(());
+void          initGridBins        PROTO(());
 void          initImageBins       PROTO((Catalog *catalog, int Ncatalog, int doImageList));
 void          initImagesSubset    PROTO((ImageSubset *input, off_t *line_number, off_t N));
 void          initImages          PROTO((Image *input, off_t *LineNumber, off_t N));
 void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
-void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
+void          initMosaicMcal      PROTO((Image *image, off_t Nimage));
 void          initMosaics         PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage));
 void          initTGroups         PROTO((Image *subset, off_t Nsubset));
@@ -621,5 +617,4 @@
 void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
 void          plot_defaults       PROTO((Graphdata *graphdata));
-void          plot_grid           PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
 void          plot_images         PROTO((void));
 void          plot_list           PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *format, ...) OHANA_FORMAT(printf, 6, 7) );
@@ -639,5 +634,7 @@
 int           reload_images       PROTO((FITS_DB *db));
 int           setExclusions       PROTO((Catalog *catalog, int Ncatalog, int verbose));
-void          setMgrid            PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
+void          setMgrid            PROTO((Catalog *catalog, int Ncatalog));
+void          resetMgrid          PROTO((void));
+
 void          setMcalFromMosaics  PROTO((void));
 void          setMcalFromTGroups  PROTO((void));
Index: /branches/eam_branches/relphot.20210521/src/ConfigInit.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/ConfigInit.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/ConfigInit.c	(revision 41603)
@@ -63,8 +63,4 @@
 
   GetConfig (config, "ZERO_PT",                "%lf", 0, &ZERO_POINT);
-
-  GetConfig (config, "RELPHOT_GRID_X",         "%d",  0, &RELPHOT_GRID_X);
-  GetConfig (config, "RELPHOT_GRID_Y",         "%d",  0, &RELPHOT_GRID_Y);
-  GetConfig (config, "RELPHOT_GRID_BINNING",   "%d",  0, &RELPHOT_GRID_BINNING);
   GetConfig (config, "CAMERA_CONFIG",          "%s",  0, CameraConfig);
   GetConfig (config, "MOSAICNAME",             "%s",  0, MOSAICNAME);
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;
+}
Index: /branches/eam_branches/relphot.20210521/src/ImageOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/ImageOps.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/ImageOps.c	(revision 41603)
@@ -319,54 +319,7 @@
 }
 
-int findCCD (off_t idx, off_t meas, int cat, MeasureTiny *measure) {
-
-  int ccdnum;
-  double X, Y;
-  char *pname, *filter, *p, base[256];
-
-  /* identify the ccd on the basis of the photcode name */
-  pname = GetPhotcodeNamebyCode (image[idx].photcode);
-
-  // skip measurements which do not match one of the requested photcodes (
-  // (do we not already exclude in bcatalog -- maybe needed for reload_objects?
-  int Ns = GetActivePhotcodeIndex (measure[0].photcode);
-  if (Ns < 0) return FALSE;
-
-  filter = photcodes[Ns][0].name;
-  sprintf (base, "%s.%s.", MOSAICNAME, filter);
-  if (strncmp (pname, base, strlen (base))) return (FALSE);
-  p = pname + strlen(base);
-
-  /* p is pointing at the DETECTOR ID part of the photcode: CAMERA.FILTER.DETECTOR 
-     for now, we will let this be of the form SSSSnn where SSS is an arbitrary string
-     and nn is an integer -> ccdnum */
-      
-  /// XXX we have an inconsistency here wrt dvo.layout and dvo.photcode
-
-  while (!isdigit(*p) && *p) p++;
-  if (*p == 0) return (FALSE);
-  ccdnum = atoi (p);
-
-  /* ccdnum is an integer, but not necessarily a sequence number.  when the camera layout is
-     constructed, there will be null values for undefined ccdnums */
-
-  // old code to add this measurement to the grid cell for this chip 
-  // ave = measureT[0].averef;
-  // ra  = catalog[cat].averageT[ave].R - measureT[0].dR / 3600.0;
-  // dec = catalog[cat].averageT[ave].D - measureT[0].dD / 3600.0;
-  // RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
-
-  // XXX we can now use these values (but need to be careful about old formats)
-  X = measure[0].Xccd;
-  Y = measure[0].Yccd;
-  setGridMeasure (meas, cat, X, Y, ccdnum);
-
-  return (TRUE);
-}
-
 void matchImage (Catalog *catalog, off_t meas, int cat, int doImageList) {
 
   off_t idx, ID;
-  int status;
   MeasureTiny *measure;
   
@@ -380,12 +333,4 @@
   }
   catalog[cat].measureT[meas].myDet = TRUE;
-
-  if (USE_GRID) {
-    status = findCCD (idx, meas, cat, measure);
-    if (!status) {
-      if (VERBOSE2) fprintf (stderr, "failed to determine CCD for "OFF_T_FMT", %d\n",  meas, cat);
-      return;
-    }
-  }
 
   // index for (catalog, measure) -> image
@@ -690,5 +635,5 @@
 	  continue;
       }
-      float Mgrid = getMgrid (m, c);
+      float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
       if (isnan(Mgrid)) {
 	  Ngrid++;
@@ -889,5 +834,5 @@
 	  continue;
       }
-      float Mgrid = getMgrid (m, c);
+      float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
       if (isnan(Mgrid)) {
 	  Ngrid++;
@@ -1179,5 +1124,5 @@
       Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
       if (isnan(Mgrp)) continue;
-      Mgrid = getMgrid (m, c);
+      Mgrid = getMgridTiny (&catalog[c].measureT[m]);
       if (isnan(Mgrid)) continue;
       N++;
Index: /branches/eam_branches/relphot.20210521/src/MagResidSave.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/MagResidSave.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/MagResidSave.c	(revision 41603)
@@ -192,4 +192,6 @@
   CHECK_STATUS (!status, "ERROR: problem closing meanmags file %s\n", filename);
 
+  FitDataSetFree (&psfStars);
+
   return TRUE;
 }
@@ -225,5 +227,5 @@
     if (isnan(Mgrp)) continue;
 
-    float Mgrid = getMgrid (m, c);
+    float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
     if (isnan(Mgrid)) continue;
 
Index: /branches/eam_branches/relphot.20210521/src/MosaicOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/MosaicOps.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/MosaicOps.c	(revision 41603)
@@ -50,104 +50,4 @@
 # undef COMPARE
 
-}
-
-/* find mosaic frames (unique time periods & photcode name matches mosaic) */
-void initMosaics_old (Image *image, off_t Nimage) {
-
-  off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
-  unsigned int start, stop;
-  char *pname;
-
-  if (!MOSAIC_ZEROPT) return;
-
-  Nmosaic = 0;
-  NMOSAIC = 10;
-  ALLOCATE (mosaic, Mosaic, NMOSAIC);
-
-  ALLOCATE (MosaicToImage, off_t *, NMOSAIC);
-  ALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
-  ALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
-
-  ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs
-
-  /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no 
-   * entry in the image table that represents this mosaic.  Instead, it is an
-   * internal construct that defines a group of related images 
-   */
-
-  /* generate list of unique mosaics */
-  for (i = 0; i < Nimage; i++) {
-    ImageToMosaic[i] = -1;
-
-    /* select valid mosaic images by photcode */
-    pname = GetPhotcodeNamebyCode (image[i].photcode);
-    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
-    if (status) continue;
-
-    /* set image time range */
-    start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
-    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
-
-    /* find existing mosaic with this time range */
-    found = FALSE;
-    for (j = 0; !found && (j < Nmosaic); j++) { 
-      if (stop  < mosaic[j].start) continue;
-      if (start > mosaic[j].stop)  continue;
-      found = TRUE;
-
-      // add reference from image to mosaic
-      ImageToMosaic[i] = j;
-
-      /* add image to mosaic image list */
-      MosaicToImage[j][MosaicN_Image[j]] = i;
-      MosaicN_Image[j] ++;
-      if (MosaicN_Image[j] == MosaicN_IMAGE[j]) {
-	MosaicN_IMAGE[j] += 10;
-	REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]);
-      }
-
-    }
-    if (found) continue;
-    
-    /* a new mosaic, define ranges -- preserve the original values incase this image is not used */
-    mosaic[Nmosaic].start     = start;
-    mosaic[Nmosaic].stop      = stop;
-    mosaic[Nmosaic].McalPSF   = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
-    mosaic[Nmosaic].McalAPER  = 0.0; // note : mosaic stores only offsets relative to the original image values
-    mosaic[Nmosaic].dMcal     = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
-    mosaic[Nmosaic].dMsys     = 0.0;
-    mosaic[Nmosaic].McalChiSq = 0.0;// NAN or 0.0?
-    mosaic[Nmosaic].flags     = image[i].flags;
-    mosaic[Nmosaic].secz      = image[i].secz;
-    mosaic[Nmosaic].photcode  = GetPhotcodeEquivCodebyCode (image[i].photcode);
-
-    // XXX do we need to do something about flag consistency across a mosaic?
-
-    /* add image to mosaic image list */
-    MosaicN_IMAGE[Nmosaic] = 10;
-    MosaicN_Image[Nmosaic] = 1;
-    ALLOCATE (MosaicToImage[Nmosaic], off_t, MosaicN_IMAGE[Nmosaic]);
-    MosaicToImage[Nmosaic][0] = i;
-
-    // add reference from image to mosaic
-    ImageToMosaic[i] = Nmosaic;
-    
-    Nmosaic ++;
-    if (Nmosaic == NMOSAIC) {
-      NMOSAIC += 10;
-      REALLOCATE (mosaic,  Mosaic,  NMOSAIC);
-      REALLOCATE (MosaicToImage,  off_t *, NMOSAIC);
-      REALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
-      REALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
-    }
-  }
-
-  // free this or not?
-  free (MosaicN_IMAGE);
-
-  initMosaicGrid (image, Nimage);
-
-  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
-  return;
 }
 
@@ -318,5 +218,5 @@
   free (startTimesMosaic);
 
-  initMosaicGrid (subset, Nsubset);
+  initMosaicMcal (subset, Nsubset);
 
   fprintf (stderr, "matched %d images to %d mosaics, %d simple chips not matched to mosaics\n", (int) (Nsubset - Nsimple), (int) Nmosaic, (int) Nsimple);
@@ -460,5 +360,5 @@
 
   if (mergeMcal) {
-    initMosaicGrid (image, Nimage);
+    initMosaicMcal (image, Nimage);
   }
 
@@ -511,4 +411,7 @@
 }
 
+// this function sets mosaic->coords to the median of the individual chips.  This
+// coordinate frame is used by the parallel region analysis to assign exposures (mosaics)
+// to specific machines by a single center (rather than individual chips)
 void setMosaicCenters (Image *image, off_t Nimage) {
   OHANA_UNUSED_PARAM(Nimage);
@@ -599,95 +502,24 @@
     mosaic[i].coords.cdelt1 = 1.0 / 3600.0;
     mosaic[i].coords.cdelt2 = 1.0 / 3600.0;
-
-    mosaic[i].McalPSF   = 0.0; 
-    mosaic[i].McalAPER  = 0.0; 
-    mosaic[i].dMcal     = 0.0; 
-    mosaic[i].dMsys     = 0.0;
-    mosaic[i].McalChiSq = 0.0; 
   }
   return;
 }
 
-void initMosaicGrid (Image *image, off_t Nimage) {
+void initMosaicMcal (Image *image, off_t Nimage) {
   OHANA_UNUSED_PARAM(Nimage);
 
-  /* find max dR, dD range for all mosaics */
-  /* define mosaic.coords to cover dR, dD */
-  /* send results to initGridBins */
-
-  off_t i, j, m, NX, NY;
-  int dXmax, dYmax;
-  double dS, dX, dY;
-  double R, D, Rmin, Rmax, Dmin, Dmax;
   double McalPSF, McalAPER, dMcal, McalChiSq;
 
   fprintf (stderr, "*** moving Mcal from image.Mcal to mosaic.Mcal ***\n");
 
-  dXmax = dYmax = 0.0;
-  for (i = 0; i < Nmosaic; i++) {
-    Dmin = Rmin =  360.0;
-    Dmax = Rmax = -360.0;
-    dS = 0.0;
+  for (off_t i = 0; i < Nmosaic; i++) {
     McalPSF = McalAPER = dMcal = McalChiSq = 0;
-    for (j = 0; j < MosaicN_Image[i]; j++) {
-      m = MosaicToImage[i][j];
-      NX = image[m].NX;
-      NY = image[m].NY;
-      dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
-
-      OhanaProjection proj = GetProjection (image[m].coords.ctype);
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-	XY_to_LM (&R, &D, 0.0, 0.0, &image[m].coords);
-      } else {
-	XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-	XY_to_LM (&R, &D, (double) NX, 0.0, &image[m].coords);
-      } else {
-	XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-	XY_to_LM (&R, &D, (double) NX, (double) NY, &image[m].coords);
-      } else {
-	XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-	XY_to_LM (&R, &D, 0.0, (double) NY, &image[m].coords);
-      } else {
-	XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      // XXX : this probably does not handle mosaics at RA = 0,360 well
-
-      /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
-
-      // XXX: how does this work with UBERCAL?  We want to keep the Mcal values supplied by ubercal, but 
-      // solve for a single offset for each exposure (Mosaic.Mcal).  
-      // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not)
-      // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value
+    for (off_t j = 0; j < MosaicN_Image[i]; j++) {
+      off_t m = MosaicToImage[i][j];
 
       if (!isfinite(image[m].McalPSF)) {
-	image[m].McalPSF = 0.0;
-	image[m].McalAPER = 0.0;
-	image[m].dMcal = 0.0;
+	image[m].McalPSF   = 0.0;
+	image[m].McalAPER  = 0.0;
+	image[m].dMcal     = 0.0;
 	image[m].McalChiSq = 0.0;
 	fprintf (stderr, "warning: resetting NAN value for Mcal %s\n", image[m].name);
@@ -703,15 +535,5 @@
       image[m].dMcal     = NAN;
       image[m].McalChiSq = NAN;
-
-    }
-    dS /= MosaicN_Image[i];
-
-    InitCoords (&mosaic[i].coords, "DEC--TAN");
-    mosaic[i].coords.crval1 = Rmin;
-    mosaic[i].coords.crval2 = Dmin;
-    mosaic[i].coords.cdelt1 = dS;
-    mosaic[i].coords.cdelt2 = dS;
-
-    RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
+    }
 
     mosaic[i].McalPSF   = McalPSF / MosaicN_Image[i];
@@ -720,9 +542,4 @@
     mosaic[i].McalChiSq = McalChiSq / MosaicN_Image[i];
   }
-  if (!USE_GRID) return;
-
-  dXmax = MAX (dXmax, dX);
-  dYmax = MAX (dYmax, dY);
-  initGrid (dXmax, dYmax);
   return;
 }
@@ -1231,5 +1048,5 @@
       continue;
     }
-    float Mgrid = getMgrid (m, c);
+    float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
     if (isnan(Mgrid)) {
       info->Ngrid ++;
@@ -1382,5 +1199,5 @@
     float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);
     // float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
-    float Mgrid    = getMgrid (m, c);
+    float Mgrid    = getMgridTiny (&catalog[c].measureT[m]);
     float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
     float Mflat    = getMflat (m, c, flatcorr, catalog);
@@ -1481,4 +1298,5 @@
 	   summary.Nskip);
   free (threadinfo);
+  // XXX SetMmosInfoFree (&summary);
 
   npass_output ++;
@@ -1811,5 +1629,5 @@
       Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
       if (isnan(Mgrp)) continue;
-      Mgrid = getMgrid (m, c);
+      Mgrid = getMgridTiny (&catalog[c].measureT[m]);
       if (isnan(Mgrid)) continue;
       Mrel = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
Index: /branches/eam_branches/relphot.20210521/src/StarOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/StarOps.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/StarOps.c	(revision 41603)
@@ -495,5 +495,5 @@
 	float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
 	if (isnan(Mgrp)) continue;
-	float Mgrid = getMgrid (m, i);
+	float Mgrid = getMgridTiny (&catalog[i].measureT[m]);
 	if (isnan(Mgrid)) continue;
 
@@ -737,5 +737,5 @@
 	  Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
 	  if (isnan(Mgrp)) { Ngrp ++; continue; }
-	  Mgrid = getMgrid (m, i);
+	  Mgrid = getMgridTiny (&catalog[i].measureT[m]);
 	  if (isnan(Mgrid)) { Ngrid ++; continue; }
 
@@ -788,5 +788,5 @@
 	  Mmos  = getMmos  (m, i);
 	  if (isnan(Mmos)) continue;
-	  Mgrid = getMgrid (m, i);
+	  Mgrid = getMgridTiny (&catalog[i].measureT[m]);
 	  if (isnan(Mgrid)) continue;
 
@@ -869,5 +869,5 @@
 	Mmos = getMmos  (m, i);
 	if (isnan(Mmos)) { continue; }
-	Mgrid = getMgrid (m, i);
+	Mgrid = getMgridTiny (&catalog[i].measureT[m]);
 	if (isnan(Mgrid)) { continue;}
 	N++;
@@ -1105,5 +1105,5 @@
 
       float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
-      // float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
+      // float Mgrid    = getMgridTiny (&catalog[c].measureT[m]);                    // camera offset (deprecated?)
       // float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
 
Index: /branches/eam_branches/relphot.20210521/src/TGroupOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/TGroupOps.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/TGroupOps.c	(revision 41603)
@@ -853,5 +853,5 @@
       continue;
     }
-    float Mgrid = getMgrid (m, c);
+    float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
     if (isnan(Mgrid)) {
       info->Ngrid ++;
@@ -1013,5 +1013,5 @@
     float Mmos     = getMmos  (m, c);                    // mosaic zero point
 
-    float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
+    float Mgrid    = getMgridTiny (&catalog[c].measureT[m]);                    // camera offset (deprecated?)
     float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
     float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
@@ -1124,4 +1124,6 @@
   free (threadinfo);
 
+  // XXX not allocated SetMgrpInfoFree (&summary);
+
   npass_output ++;
 
@@ -1268,5 +1270,7 @@
 
 static float MinMaxChiSq = NAN;
+static float MaxMaxChiSq = 6.0;
 static float MinMaxScatter = NAN;
+static float MaxMaxScatter = 0.1;
 
 void clean_tgroups () {
@@ -1301,10 +1305,10 @@
   float ChiSqUpper90 = stats.Upper90; 
   if (isnan (MinMaxChiSq)) MinMaxChiSq = 2.0*stats.median;
-  float MaxChiSq = MAX (MinMaxChiSq, ChiSqUpper90);
+  float MaxChiSq = MIN(MaxMaxChiSq, MAX (MinMaxChiSq, ChiSqUpper90)); // the upper limit should be between MinMaxChiSq and MaxMaxChiSq 
 
   liststats (slist, NULL, NULL, N, &stats);
   float ScatterUpper90 = stats.Upper90; 
   if (isnan (MinMaxScatter)) MinMaxScatter = 2.0*stats.median;
-  float MaxScatter = MAX (MinMaxScatter, ScatterUpper90);
+  float MaxScatter = MIN(MaxMaxScatter, MAX (MinMaxScatter, ScatterUpper90));
 
   fprintf (stderr, "TGROUPS: Max ChiSq: %f, Max Scatter: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", MaxChiSq, MaxScatter, ChiSqUpper90, ScatterUpper90);
@@ -1501,5 +1505,5 @@
 
 	float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
-	// float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
+	// float Mgrid    = getMgridTiny (&catalog[c].measureT[m]); // camera offset (deprecated?)
 	// float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
 
Index: /branches/eam_branches/relphot.20210521/src/ZeroPointModes.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/ZeroPointModes.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/ZeroPointModes.c	(revision 41603)
@@ -49,5 +49,5 @@
       clean_mosaics(); // do this on every pass 
     }
-    if ((CurrentLoop >= 18) && (CurrentLoop <= 999)) {
+    if ((CurrentLoop >= 18) && (CurrentLoop < 21)) {
       // after iterating a few times on the TGroup & Mosaic zero points:
       // * identify the good / bad Mosaics
@@ -60,4 +60,16 @@
       clean_images(); // do this on every pass 
     }
+    if ((CurrentLoop >= 21) && (CurrentLoop <= 999)) {
+      // after iterating a few times on the TGroup & Mosaic zero points:
+      // * identify the good / bad Mosaics
+      // * fit the chips from the bad Mosaics
+      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights
+      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_BAD_NIGHT_GOOD_MOSAIC; // only fit good mosaics on bad nights
+      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_BAD_NIGHT_BAD_MOSAIC; // only fit bad mosaics on bad nights
+      GRID_ZPT_MODE  = GRID_ZPT_MODE_ALL; // only fit bad mosaics on bad nights
+      clean_tgroups(); // do this on every pass 
+      clean_mosaics(); // do this on every pass 
+      clean_images(); // do this on every pass 
+    }
   }
 
@@ -68,5 +80,5 @@
       MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
     }
-    if ((CurrentLoop > 4) && (CurrentLoop <= 999)) {
+    if ((CurrentLoop > 4) && (CurrentLoop <= 8)) {
       // after iterating a few times on the TGroup zero points:
       // * identify the photometric nights
@@ -75,4 +87,14 @@
       MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
       IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_BAD_NIGHT;
+      clean_tgroups(); // do this on every pass or just sometimes?
+    }
+    if ((CurrentLoop > 8) && (CurrentLoop <= 999)) {
+      // after iterating a few times on the TGroup zero points:
+      // * identify the photometric nights
+      // * fit the images from the non-photometric nights
+      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights
+      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
+      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_BAD_NIGHT;
+      GRID_ZPT_MODE  = GRID_ZPT_MODE_ALL; // only fit bad mosaics on bad nights
       clean_tgroups(); // do this on every pass or just sometimes?
     }
Index: /branches/eam_branches/relphot.20210521/src/args.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/args.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/args.c	(revision 41603)
@@ -198,11 +198,4 @@
   }
 
-  NGRID = 8;
-  if ((N = get_argument (argc, argv, "-ngrid"))) {
-    remove_argument (N, &argc, argv);
-    NGRID = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-  }
-
   TEST_IMAGE1 = -1;
   if ((N = get_argument (argc, argv, "-test-image1"))) {
@@ -387,9 +380,10 @@
   }
 
-  // USE_GRID is not valid for all cases, probably should be its own mode...
-  USE_GRID = FALSE;
+  // GRID_ZEROPT is not valid for all cases, probably should be its own mode...
+  GRID_ZEROPT = FALSE;
   if ((N = get_argument (argc, argv, "-grid"))) {
     remove_argument (N, &argc, argv);
-    USE_GRID = TRUE;
+    GRID_ZEROPT = TRUE;
+    GRID_ZPT_MODE  = GRID_ZPT_MODE_NONE; // start with grid off
   }
 
Index: /branches/eam_branches/relphot.20210521/src/assign_images.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/assign_images.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/assign_images.c	(revision 41603)
@@ -36,6 +36,7 @@
   if (MOSAIC_ZEROPT) {
     makeMosaics (image, Nimage, FALSE);
-
-    // center coords and Mcal, dMcal, Mchisq for the mosaics
+    MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime);
+
+    // center coords and zero Mcal, dMcal, Mchisq for the mosaics
     setMosaicCenters (image, Nimage);
     MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime);
@@ -124,9 +125,9 @@
       Mosaic *mosaic = getMosaicForImage (j); 
       if (mosaic) {
-	Rc = mosaic->coords.crval1;
-	Dc = mosaic->coords.crval2;
-	// NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the
-	// chips are located.  but, for host assignment, we rationalize to 0.0 - 360.0
-	Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0);
+    	Rc = mosaic->coords.crval1;
+    	Dc = mosaic->coords.crval2;
+    	// NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the
+    	// chips are located.  but, for host assignment, we rationalize to 0.0 - 360.0
+    	Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0);
       }
     }
Index: /branches/eam_branches/relphot.20210521/src/bcatalog.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/bcatalog.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/bcatalog.c	(revision 41603)
@@ -112,5 +112,5 @@
 
       // require 0x01 in photFlags (fitted with a PSF)
-      if ((catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; }
+// XXX EAM : test is inconsistent -- fix      if ((catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; }
 
       // very loose cut on PSF - Kron 
Index: /branches/eam_branches/relphot.20210521/src/extra.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/extra.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/extra.c	(revision 41603)
@@ -1,3 +1,4 @@
 # include "relphot.h"
+# define DISABLE_CORNER_GPC1 0
 
 int isMosaicChip (int photcode) {
@@ -14,5 +15,5 @@
 int whichGPC1filter (int photcode) {
 
-# if (1)
+# if (DISABLE_CORNER_GPC1)
   // disable the corner chips:
   if (((photcode >  10000) && (photcode <  10077)) || (photcode == 4100)) return PS1_g; // g-band 
@@ -37,5 +38,5 @@
 int isGPC1chip (int photcode) {
 
-# if (1)
+# if (DISABLE_CORNER_GPC1)
   if (((photcode >  10000) && (photcode <  10077)) || (photcode == 4100)) return TRUE; // g-band
   if (((photcode >  10100) && (photcode <  10177)) || (photcode == 4200)) return TRUE; // r-band
Index: /branches/eam_branches/relphot.20210521/src/initialize.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/initialize.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/initialize.c	(revision 41603)
@@ -27,9 +27,4 @@
   }    
     
-  if (USE_GRID && (Nphotcodes > 1)) {
-    fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n");
-    exit (1);
-  }
-
   IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
   STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
@@ -67,9 +62,4 @@
 
     fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF);
-    fprintf (stderr, "GRID_X: %d, GRID_Y: %d, BINNING: %d == Nmx: %d, Nmy: %d\n", 
-	     RELPHOT_GRID_X, 
-	     RELPHOT_GRID_Y, 
-	     RELPHOT_GRID_BINNING, 
-	     (int)(RELPHOT_GRID_X/RELPHOT_GRID_BINNING), (int)(RELPHOT_GRID_Y/RELPHOT_GRID_BINNING));
 
     fprintf (stderr, "STAR_SCATTER           %lf\n", STAR_SCATTER);
Index: /branches/eam_branches/relphot.20210521/src/plot_scatter.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/plot_scatter.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/plot_scatter.c	(revision 41603)
@@ -44,5 +44,5 @@
 		Mmos = getMmos  (m, i);
 		if (isnan(Mmos)) continue;
-		Mgrid = getMgrid (m, i);
+		Mgrid = getMgridTiny (&catalog[i].measureT[m]);
 		if (isnan(Mgrid)) continue;
 
Index: /branches/eam_branches/relphot.20210521/src/reload_catalogs.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/reload_catalogs.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/reload_catalogs.c	(revision 41603)
@@ -126,5 +126,5 @@
     freeMosaicBins (1, FALSE);
     freeTGroupBins (1);
-    freeGridBins (1);
+    freeGridBins ();
     TIMESTAMP(time7);
 
Index: /branches/eam_branches/relphot.20210521/src/relphot_images.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/relphot_images.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/relphot_images.c	(revision 41603)
@@ -54,4 +54,6 @@
   FlatCorrectionTable *flatcorr = NULL;
 
+  initGridBins (); // allocates the empty array of corrections (elements are build below)
+
   if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) {
     /* load catalog data from region files (hostID is 0 since we are not a client */
@@ -68,5 +70,4 @@
     initTGroupBins (catalog, Ncatalog);
 
-    initGridBins   (catalog, Ncatalog);
     initMrel (catalog, Ncatalog);
 
@@ -74,9 +75,11 @@
     MARKTIME("-- set up image indexes: %f sec\n", dtime);
 
-    findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
+    findMosaics (catalog, Ncatalog, TRUE);
     MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
 
     findTGroups (catalog, Ncatalog);
     MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
+
+    initGrid (); // allocate grid correction entries for existing photcodes
 
     SAVEPLOT = FALSE;
@@ -91,22 +94,4 @@
     }
 
-    // if we are measuring the flat-field correction grid, we need to perform a number of iterations first:
-    if (USE_GRID) {
-      int star_toofew;
-
-      // until we finish the grid analysis, do not reject stars out-of-hand based on ID_OBJ_FEW
-      // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel
-      star_toofew = STAR_TOOFEW;
-      STAR_TOOFEW = 0;
-      for (i = 0; i < NGRID; i++) {
-	STAR_BAD = ID_OBJ_POOR;
-	setMrel  (catalog, Ncatalog, flatcorr);
-	STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
-	setMgrid (catalog, flatcorr);
-      }
-      STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
-      STAR_TOOFEW = star_toofew;
-    }
-
     /* determine fit values */
     for (i = 0; i < NLOOP; i++) {
@@ -119,10 +104,9 @@
       setMgrp  (catalog, flatcorr);
 
-      setMgrid (catalog, flatcorr);
+      setMgrid (catalog, Ncatalog);
       MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
     
       if (PLOTSTUFF) {
 	plot_scatter (catalog, Ncatalog, flatcorr); 
-	plot_grid (catalog, flatcorr); 
 	plot_mosaics ();
 	plot_images ();
@@ -140,5 +124,4 @@
     if (PLOTSTUFF) {
       plot_scatter (catalog, Ncatalog, flatcorr); 
-      plot_grid (catalog, flatcorr); 
       plot_mosaics ();
       plot_images ();
@@ -147,5 +130,5 @@
     }
   
-    if (USE_GRID) dump_grid ();
+    // if (GRID_ZEROPT) dump_grid ();
 
     MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
@@ -169,6 +152,4 @@
     freeMosaicBins (Ncatalog, TRUE);
     freeTGroupBins (Ncatalog);
-    freeGridBins (Ncatalog);
-
     // end of if (NLOOP > 0) block : this loop determines the offsets per chip 
   } else { 
@@ -205,4 +186,6 @@
   freeImages (db.ftable.buffer);
   freeMosaics ();
+  freeGridBins ();
+
   gfits_db_free (&db);
 
Index: /branches/eam_branches/relphot.20210521/src/relphot_objects.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/relphot_objects.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/relphot_objects.c	(revision 41603)
@@ -147,5 +147,5 @@
       freeImageBins (1, FALSE);
       freeMosaicBins (1, FALSE);
-      freeGridBins (1);
+      freeGridBins ();
       continue;
     }
@@ -183,5 +183,5 @@
     freeImageBins (1, FALSE);
     freeMosaicBins (1, FALSE);
-    freeGridBins (1);
+    freeGridBins ();
   }
 
Index: /branches/eam_branches/relphot.20210521/src/relphot_parallel_images.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/relphot_parallel_images.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/relphot_parallel_images.c	(revision 41603)
@@ -161,5 +161,5 @@
   freeImageBins(Ncatalog, TRUE);
   freeMosaicBins (Ncatalog, TRUE);
-  freeGridBins (Ncatalog);
+  freeGridBins ();
   freeImages((char *)image);
   free (image);
Index: /branches/eam_branches/relphot.20210521/src/setMrelCatalog.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/setMrelCatalog.c	(revision 41602)
+++ /branches/eam_branches/relphot.20210521/src/setMrelCatalog.c	(revision 41603)
@@ -243,5 +243,5 @@
       Mgrp  = getMgrp  (meas, cat, measureT[k].airmass, &dMgrp);
       if (isnan(Mgrp))  SKIP_THIS_MEAS(Ngrp);
-      Mgrid = getMgrid (meas, cat);
+      Mgrid = getMgridTiny (&measureT[k]);
       if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
       Mflat = isnan (measureT[k].Mflat) ? 0.0 : measureT[k].Mflat;
@@ -732,5 +732,5 @@
       McalAPER  = USE_MCAL_PSF_FOR_STACK_APER ? getMcal  (measSeq, cat, MAG_CLASS_PSF) : getMcal  (measSeq, cat, MAG_CLASS_KRON);
       Mmos      = getMmos  (measSeq, cat);
-      Mgrid     = getMgrid (measSeq, cat);
+      Mgrid     = getMgrid (&measure[meas]);
       // XXX can Mmos and Mgrid exist for stacks?
     }