Index: /trunk/Ohana/src/relphot/include/relphot.h
===================================================================
--- /trunk/Ohana/src/relphot/include/relphot.h	(revision 20322)
+++ /trunk/Ohana/src/relphot/include/relphot.h	(revision 20323)
@@ -63,4 +63,5 @@
 
 int    NLOOP;
+int    NGRID;
 int    RESET;
 int    UPDATE;
@@ -73,8 +74,10 @@
 int    MEAS_BAD;
 int    STAR_TOOFEW;
+int    GRID_TOOFEW;
 int    IMAGE_TOOFEW;
 double IMAGE_GOOD_FRACTION;
 int    IMAGE_BAD;
 int    FREEZE_IMAGES;
+int    FREEZE_MOSAICS;
 int    USE_GRID;
 char  *OUTROOT;
Index: /trunk/Ohana/src/relphot/src/ConfigInit.c
===================================================================
--- /trunk/Ohana/src/relphot/src/ConfigInit.c	(revision 20322)
+++ /trunk/Ohana/src/relphot/src/ConfigInit.c	(revision 20323)
@@ -25,4 +25,5 @@
 
   GetConfig (config, "STAR_CHISQ",             "%lf", 0, &STAR_CHISQ);
+  GetConfig (config, "GRID_TOOFEW",            "%d",  0, &GRID_TOOFEW);
   GetConfig (config, "STAR_TOOFEW",            "%d",  0, &STAR_TOOFEW);
   GetConfig (config, "IMAGE_TOOFEW",           "%d",  0, &IMAGE_TOOFEW);
Index: /trunk/Ohana/src/relphot/src/GridOps.c
===================================================================
--- /trunk/Ohana/src/relphot/src/GridOps.c	(revision 20322)
+++ /trunk/Ohana/src/relphot/src/GridOps.c	(revision 20323)
@@ -1,8 +1,15 @@
 # include "relphot.h"
 
+typedef enum {
+  GRID_FITTED,
+  GRID_FROZEN,
+  GRID_REFERENCE,
+};
+
 static int     Ngrid;
-static float   *gridM;
-static float   *gridS;
-static int     *gridN;
+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;
 static int      gridY;
@@ -30,28 +37,7 @@
 static char *config;
 
-# ifdef GRID_V1
 void initGrid (int dX, int dY) {
 
-  int i;
-
-  /* define mosaic 2d correction grid */
-  gridX = dX / RELPHOT_GRID_BINNING + 1;
-  gridY = dY / RELPHOT_GRID_BINNING + 1;
-  Ngrid = gridX * gridY;
-
-  ALLOCATE (gridM, float, Ngrid);
-  ALLOCATE (gridS, float, Ngrid);
-  ALLOCATE (gridN, int,   Ngrid);
-  bzero (gridM, Ngrid*sizeof(float));
-  bzero (gridS, Ngrid*sizeof(float));
-  bzero (gridN, Ngrid*sizeof(int));
-
-}
-# endif
-
-# ifdef GRID_V2
-void initGrid (int dX, int dY) {
-
-  int i, N, ccdnum_max;
+  int i, N, ccdnum_max, refX, refY, refBin;
   char *p, field[64], line[256];
   int *Fx, *Fy;  /* chip flip */
@@ -72,4 +58,7 @@
   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
@@ -143,10 +132,19 @@
   ALLOCATE (gridS, float, Ngrid);
   ALLOCATE (gridN, int,   Ngrid);
-  bzero (gridM, Ngrid*sizeof(float));
-  bzero (gridS, Ngrid*sizeof(float));
-  bzero (gridN, Ngrid*sizeof(int));
-
-}
-# endif
+  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) {
@@ -208,21 +206,30 @@
 }
 
-# ifdef GRID_V1
-int setGridMeasure (int meas, int cat, double X, double Y) {
-
-  int ix, iy, i;
-
-  ix = X / RELPHOT_GRID_BINNING;
-  iy = Y / RELPHOT_GRID_BINNING;
-  if (ix < 0) goto escape;
-  if (iy < 0) goto escape;
-  if (ix >= gridX) goto escape;
-  if (iy >= gridY) goto escape;
+int setGridMeasure (int 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;
-  Ymeas[cat][meas] = Y;
+  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;
@@ -236,52 +243,7 @@
   return (TRUE);
 
-escape:
   fprintf (stderr, "error: star out of grid\n");
   exit (1);
 }
-# endif
-
-# ifdef GRID_V2
-int setGridMeasure (int 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;
-
-  Nlist[i] ++;
-  if (Nlist[i] == NLIST[i]) {
-    NLIST[i] += 100;
-    REALLOCATE (clist[i], int, NLIST[i]);
-    REALLOCATE (mlist[i], int, NLIST[i]);
-  }	
-  return (TRUE);
-
-  fprintf (stderr, "error: star out of grid\n");
-  exit (1);
-}
-# endif
 
 float getMgrid (int meas, int cat) {
@@ -294,4 +256,8 @@
   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);
@@ -301,5 +267,5 @@
 void setMgrid (Catalog *catalog) {
 
-  int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ncal, Nrel, Nsys;
+  int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ncal, Nrel, Nsys, Nfit;
   double *list, *dlist;
   float Msys, Mrel, Mcal, Mmos;
@@ -315,5 +281,5 @@
   ALLOCATE (dlist, double, Nmax);
 
-  Nbad = Ncal = Nmos = Nrel = Nsys = 0;
+  Nbad = Ncal = Nmos = Nrel = Nsys = Nfit = 0;
 
   for (i = 0; i < Ngrid; i++) {
@@ -356,11 +322,27 @@
     }
 
+    // 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, N, &stats);
     gridM[i] = stats.mean;
     gridS[i] = stats.sigma;
     gridN[i] = N;
-  }
-
-  fprintf (stderr, "grid cells having too few measurements (Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Nbad, Nmos, Ncal, Nrel, Nsys);
+    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);
@@ -531,44 +513,4 @@
   gfits_free_matrix (&matrix);
 
-# ifdef GRID_V1
-  /* calculate pixel values for each CCD pixel, write out CCD images */
-  /* grid pixels are in RA,DEC coords, transform to image and interpolate */
-  for (i = 0; i < Nimage; i++) {
-    image = getimage (imlist[i]);
-    pname = GetPhotcodeNamebyCode (image[0].photcode);
-
-    /* this is kind of bogus... */
-    /* pname is CAMERA.FILTER.CCD, grab the CCD */
-    p = strrchr (pname, '.');
-    if (p == (char *) NULL) {
-      fprintf (stderr, "error parsing photcode %s\n", pname);
-      exit (2);
-    }
-    p++;
-    sprintf (ccdname, "ccd%s", p);
-
-    gfits_modify (&theader, "EXTNAME", "%s", 1, ccdname);
-    gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
-    gfits_modify (&theader, "PHOTCODE", "%s", 1, pname);
-    gfits_modify (&theader, "NX", "%d", 1, image[i].NX);
-    gfits_modify (&theader, "NY", "%d", 1, image[i].NY);
-    write_coords (&theader, &image[0].coords);
-
-    Nx = 2 * image[i].NX / RELPHOT_GRID_BINNING;
-    Ny = 2 * image[i].NY / RELPHOT_GRID_BINNING;
-    theader.Naxis[0] = Nx;
-    theader.Naxis[1] = Ny;
-    gfits_modify (&theader, "NAXIS1", "%d", 1, Nx);
-    gfits_modify (&theader, "NAXIS2", "%d", 1, Ny);
-    gfits_create_matrix  (&theader, &matrix);
-
-    InterpolateGrid ((float *)matrix.buffer, Nx, Ny, &image[0].coords, &refmosaic[0].coords);
-    gfits_fwrite_header (f, &theader);
-    gfits_fwrite_matrix (f, &matrix);
-    gfits_free_matrix (&matrix);
-  }
-# endif
-
-# ifdef GRID_V2
   /* calculate value for each CCD pixel, write out CCD images */
   /* grid pixels are tied to detector pixels, but are flipped to match focal plane */
@@ -610,5 +552,4 @@
     gfits_free_matrix (&matrix);
   }
-# endif
 
   free (filename);
Index: /trunk/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- /trunk/Ohana/src/relphot/src/ImageOps.c	(revision 20322)
+++ /trunk/Ohana/src/relphot/src/ImageOps.c	(revision 20323)
@@ -178,5 +178,4 @@
     if (measure[0].t > stop[i]) continue;
     
-# ifdef GRID_V2 /* this section is added to support GridOps.v2.c */
     if (USE_GRID) {
 
@@ -199,13 +198,9 @@
 	 constructed, there will be null values for undefined ccdnums */
 
-# if (0)
-      /* add this measurement to the grid cell for this chip */
-      ave = measure[0].averef;
-      ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
-      dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
-       
-      /* X,Y always positive-definite in range 0,0 - dX, dY */
-      RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
-# endif
+      // old code to add this measurement to the grid cell for this chip 
+      // ave = measure[0].averef;
+      // ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
+      // dec = catalog[cat].average[ave].D - measure[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)
@@ -214,5 +209,4 @@
       setGridMeasure (meas, cat, X, Y, ccdnum);
     }
-# endif
 
     bin[cat][meas] = i;
Index: /trunk/Ohana/src/relphot/src/MosaicOps.c
===================================================================
--- /trunk/Ohana/src/relphot/src/MosaicOps.c	(revision 20322)
+++ /trunk/Ohana/src/relphot/src/MosaicOps.c	(revision 20323)
@@ -320,4 +320,5 @@
 
   if (!MOSAIC_ZEROPT) return (FALSE);
+  if (FREEZE_MOSAICS) return (FALSE);
 
   image = getimages (&N);
@@ -453,4 +454,5 @@
   bzero (&stats, sizeof (StatType));
   if (!MOSAIC_ZEROPT) return (stats);
+  if (FREEZE_MOSAICS) return (stats); 
 
   ALLOCATE (list, double, Nmosaic);
@@ -481,4 +483,5 @@
   bzero (&stats, sizeof (StatType));
   if (!MOSAIC_ZEROPT) return (stats);
+  if (FREEZE_MOSAICS) return (stats); 
 
   ALLOCATE (list, double, Nmosaic);
@@ -523,4 +526,5 @@
   bzero (&stats, sizeof (StatType));
   if (!MOSAIC_ZEROPT) return (stats);
+  if (FREEZE_MOSAICS) return (stats); 
 
   ALLOCATE (list, double, Nmosaic);
@@ -551,4 +555,5 @@
 
   if (!MOSAIC_ZEROPT) return;
+  if (FREEZE_MOSAICS) return;
 
   if (VERBOSE) fprintf (stderr, "marking poor mosaics\n");
@@ -645,4 +650,5 @@
 
   if (!MOSAIC_ZEROPT) return;
+  if (FREEZE_MOSAICS) return;
 
   ALLOCATE (xlist, double, Nmosaic);
Index: /trunk/Ohana/src/relphot/src/StarOps.c
===================================================================
--- /trunk/Ohana/src/relphot/src/StarOps.c	(revision 20322)
+++ /trunk/Ohana/src/relphot/src/StarOps.c	(revision 20323)
@@ -85,5 +85,8 @@
 	N++;
       }
-      if (N < STAR_TOOFEW) { /* too few measurements */
+
+      // when performing the grid analysis, STAR_TOOFEW will be set to 1;
+
+      if (N <= STAR_TOOFEW) { /* too few measurements */
 	catalog[i].average[j].code |= ID_STAR_FEW;
 	Nfew ++;
@@ -367,5 +370,5 @@
 	N++;
       }
-      if (N < TOOFEW) continue;
+      if (N <= TOOFEW) continue;
 
       /* 3-sigma clip based on stats of inner 50% */
Index: /trunk/Ohana/src/relphot/src/args.c
===================================================================
--- /trunk/Ohana/src/relphot/src/args.c	(revision 20322)
+++ /trunk/Ohana/src/relphot/src/args.c	(revision 20323)
@@ -93,4 +93,16 @@
     remove_argument (N, &argc, argv);
   }
+  if ((N = get_argument (argc, argv, "-nloop"))) {
+    remove_argument (N, &argc, argv);
+    NLOOP = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  NGRID = 8;
+  if ((N = get_argument (argc, argv, "-ngrid"))) {
+    remove_argument (N, &argc, argv);
+    NGRID = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
 
   RESET = FALSE;
@@ -139,4 +151,10 @@
     remove_argument (N, &argc, argv);
     FREEZE_IMAGES = TRUE;
+  }
+
+  FREEZE_MOSAICS = FALSE;
+  if ((N = get_argument (argc, argv, "-mosfreeze"))) {
+    remove_argument (N, &argc, argv);
+    FREEZE_MOSAICS = TRUE;
   }
 
Index: /trunk/Ohana/src/relphot/src/bcatalog.c
===================================================================
--- /trunk/Ohana/src/relphot/src/bcatalog.c	(revision 20322)
+++ /trunk/Ohana/src/relphot/src/bcatalog.c	(revision 20323)
@@ -90,5 +90,5 @@
 
     // XXXX test : what checks do I need to make elsewhere to avoid problems here?
-    if (Nm < STAR_TOOFEW) { /* enough measurements in band? */
+    if (Nm <= STAR_TOOFEW) { /* enough measurements in band? */
       Nmeasure -= Nm;
       continue; 
Index: /trunk/Ohana/src/relphot/src/relphot.c
===================================================================
--- /trunk/Ohana/src/relphot/src/relphot.c	(revision 20322)
+++ /trunk/Ohana/src/relphot/src/relphot.c	(revision 20323)
@@ -70,4 +70,21 @@
   }
 
+  // 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_STAR_FEW
+      // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel
+      star_toofew = STAR_TOOFEW;
+      for (i = 0; i < NGRID; i++) {
+	  STAR_BAD = ID_STAR_POOR;
+	  setMrel  (catalog, Ncatalog);
+	  STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
+	  setMgrid (catalog);
+      }
+      STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
+      STAR_TOOFEW = star_toofew;
+  }
+
   /* determine fit values */
   for (i = 0; i < NLOOP; i++) {