Index: /trunk/Ohana/src/relphot/Makefile
===================================================================
--- /trunk/Ohana/src/relphot/Makefile	(revision 9632)
+++ /trunk/Ohana/src/relphot/Makefile	(revision 9633)
@@ -21,5 +21,5 @@
 RELPHOT = \
 $(SRC)/ConfigInit.$(ARCH).o	 \
-$(SRC)/GridOps.v2.$(ARCH).o	 \
+$(SRC)/GridOps.$(ARCH).o	 \
 $(SRC)/ImageOps.$(ARCH).o	 \
 $(SRC)/MosaicOps.$(ARCH).o	 \
Index: /trunk/Ohana/src/relphot/doc/notes.txt
===================================================================
--- /trunk/Ohana/src/relphot/doc/notes.txt	(revision 9632)
+++ /trunk/Ohana/src/relphot/doc/notes.txt	(revision 9633)
@@ -1,2 +1,41 @@
+
+2006.10.17
+
+  I am working on a partial upgrade of relphot to handle more general
+  problems than I have treated in the past.  In particular, I would
+  like to be able to have it assign the average magnitudes, regardless
+  of whether any or all images and/or objects are well-treated for
+  relative photometry.  In fact, I would like it to be able to
+  determine average magnitudes even for sources which have only
+  externally supplied data, and thus have no matching images.  I would
+  eventually like to have more flexibility about the filtering which
+  is performed on the measurements to determine the average
+  magnitudes.  Some of this work will require a bit of a more careful
+  treatment of the process.  
+
+  At the moment, we can consider the relphot process to consist of the
+  following steps:
+
+  1) load the complete catalog data
+  2) select a subset of the sources:
+     - appropriate photcode (equiv == selected photcode)
+     - in time range, if specified
+     - appropriate dophot type, if specified
+     - bright enough:
+       - mag < MAG_LIM	
+       - dMag > SIGMA_LIM
+       - ImagMin < iMag < ImagMax
+     - other restrictions
+       - within valid image region for chip
+  3) iterate to a solution for Mcal (image offsets)
+     - mark specific images as bad
+     - mark specific sources as bad
+     - mark specific measurements as bad
+     - include desired restrictions on Mcal (eg, all chips of mosaic
+       matched)
+     - only keep the image offsets
+
+  4) if the results are to be kept, recalculate the average magnitudes
+     using the image exclusions and Mcal values determined above 
 
 2006.05.03
Index: /trunk/Ohana/src/relphot/include/relphot.h
===================================================================
--- /trunk/Ohana/src/relphot/include/relphot.h	(revision 9632)
+++ /trunk/Ohana/src/relphot/include/relphot.h	(revision 9633)
@@ -6,4 +6,5 @@
 /* # define GRID_V1 */
 # define GRID_V2
+# define NO_IMAGE -100
 
 # if (0)
@@ -91,4 +92,6 @@
 
 int ImagSelect, ImagMin, ImagMax;
+
+int DophotSelect, DophotValue;
 
 double  PlotMmin, PlotMmax, PlotdMmin, PlotdMmax;
Index: /trunk/Ohana/src/relphot/src/GridOps.c
===================================================================
--- /trunk/Ohana/src/relphot/src/GridOps.c	(revision 9633)
+++ /trunk/Ohana/src/relphot/src/GridOps.c	(revision 9633)
@@ -0,0 +1,585 @@
+# include "relphot.h"
+
+static int     Ngrid;
+static float   *gridM;
+static float   *gridS;
+static int     *gridN;
+static int      gridX;
+static int      gridY;
+
+static int    **bin;
+static int    **Xmeas;
+static int    **Ymeas;
+
+static int    **clist;
+static int    **mlist;
+static int     *Nlist;
+static int     *NLIST;
+
+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 */
+  char **ccdname;
+} camera;    
+
+# 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;
+  char field[64], *config, line[256];
+
+  /* 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);
+
+  ALLOCATE (camera.Fx, int, camera.Nchip);
+  ALLOCATE (camera.Fy, int, camera.Nchip);
+  ALLOCATE (camera.Ox, int, camera.Nchip);
+  ALLOCATE (camera.Oy, int, camera.Nchip);
+  ALLOCATE (camera.ccdname, char *, camera.Nchip);
+
+  /* load per-chip parameters */
+  for (i = 0; i < camera.Nchip; i++) {
+    ALLOCATE (camera.ccdname[i], char, 256);
+    sprintf (field, "CCD.%d", i);
+    ScanConfig (config, field, "%s", 1, line);
+    sscanf (line, "%s %d %d %d %d", camera.ccdname[i], &camera.Ox[i], &camera.Oy[i], &camera.Fx[i], &camera.Fy[i]);
+  }
+  free (config);
+
+  /* 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);
+  bzero (gridM, Ngrid*sizeof(float));
+  bzero (gridS, Ngrid*sizeof(float));
+  bzero (gridN, Ngrid*sizeof(int));
+
+}
+# endif
+
+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, int, Ngrid);
+  ALLOCATE (NLIST, int, Ngrid);
+  ALLOCATE (clist, int *, Ngrid);
+  ALLOCATE (mlist, int *, Ngrid);
+
+  for (i = 0; i < Ngrid; i++) {
+    Nlist[i] = 0;
+    NLIST[i] = 100;
+    ALLOCATE (clist[i], int, NLIST[i]);
+    ALLOCATE (mlist[i], int, 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);
+}
+
+# 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;
+
+  i = ix + iy*gridX;
+
+  bin[cat][meas] = i;
+  Xmeas[cat][meas] = X;
+  Ymeas[cat][meas] = Y;
+  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);
+
+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) {
+
+  int i;
+  float value;
+
+  if (!USE_GRID) return (0);
+  i = bin[cat][meas];
+  if (i == -1) return (NO_MAG);
+
+  value = gridM[i];
+  return (value);
+}
+
+/* determine Mgrid values for all grid bins */
+void setMgrid (Catalog *catalog) {
+
+  int i, j, m, c, n, N, Nmax;
+  double *list, *dlist;
+  float Msys, Mrel, Mcal, Mmos;
+  StatType stats;
+  
+  if (!USE_GRID) return;
+
+  Nmax = Nlist[0];
+  for (i = 0; i < Ngrid; i++) {
+    Nmax = MAX (Nmax, Nlist[i]);
+  }
+  ALLOCATE (list, double, Nmax);
+  ALLOCATE (dlist, double, Nmax);
+
+  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].measure[m].flags & MEAS_BAD) continue;
+      Mcal = getMcal  (m, c);
+      if (Mcal == NO_MAG) continue;
+      if (Mcal == NO_IMAGE) continue;
+      Mmos = getMmos  (m, c);
+      if (Mmos == NO_MAG) continue;
+      if (Mmos == NO_IMAGE) continue;
+      if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
+      
+      n = catalog[c].measure[m].averef;
+      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+      list[N] = Msys - Mrel - Mcal - Mmos;
+      dlist[N] = MAX (catalog[c].measure[m].dM_PS, MIN_ERROR);
+      N++;
+    }
+
+    liststats (list, dlist, N, &stats);
+    gridM[i] = stats.mean;
+    gridS[i] = stats.sigma;
+    gridN[i] = N;
+  }
+  free (list);
+  free (dlist);
+}
+
+void plot_grid (Catalog *catalog) {
+
+  int i, j, m, c, n, N, Narea;
+  float Msys, Mrel, Mcal, Mmos;
+  double *xlist, *Mlist, *dlist, *ylist;
+  Graphdata graphdata;
+
+  if (!USE_GRID) return;
+
+  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].measure[m].flags & MEAS_BAD) {
+	Narea ++;
+	continue;
+      }
+      Mcal  = getMcal  (m, c);
+      if (Mcal == NO_MAG) continue;
+      if (Mcal == NO_IMAGE) continue;
+      Mmos  = getMmos  (m, c);
+      if (Mmos == NO_MAG) continue;
+      if (Mmos == NO_IMAGE) continue;
+      if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
+
+      n = catalog[c].measure[m].averef;
+      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+
+      xlist[N] = Xmeas[c][m];
+      ylist[N] = Ymeas[c][m];
+      Mlist[N] = Msys - Mrel - Mcal - Mmos;
+      dlist[N] = Msys - Mrel - Mcal - Mmos - 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", "XdM.png");
+  plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "XdMf.png");
+  plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "YdMf.png");
+
+  plot_defaults (&graphdata);
+  plot_list (&graphdata, xlist, ylist, N, "X vs Y", "XY.png");
+
+  free (ylist);
+  free (xlist);
+  free (Mlist);
+  free (dlist);
+
+}
+
+void dump_grid () { 
+
+  int i, j, Nimage;
+  int *imlist;
+  FILE *f;
+  Header header, theader;
+  Matrix matrix;
+  Mosaic *refmosaic;
+    
+  /* select reference mosaic image */
+  imlist = SelectRefMosaic (&refmosaic, &Nimage);
+
+  /* we are writing to this file */
+  f = fopen ("mosaic.fits", "w");
+  if (f == (FILE *) NULL) { 
+    fprintf (stderr, "cannot open %s for output\n", "mosaic.fits");
+    return;
+  }
+
+  /* create empty phu */
+  gfits_init_header (&header);
+  header.extend = TRUE;
+  gfits_create_header (&header);
+  gfits_create_matrix (&header, &matrix);
+  gfits_modify (&header, "NEXTEND", "%d", 1, Nimage + 3);
+  gfits_modify (&header, "FILTER", "%s", 1, photcode[0].name);
+  gfits_modify (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis");
+  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, photcode[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, photcode[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, photcode[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);
+
+# 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].source);
+
+    /* 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 */
+  for (i = 0; i < camera.Nchip; i++) {
+    int ix, iy, x, y, X, Y, bin;
+
+    gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[i]);
+    gfits_modify (&theader, "FILTER", "%s", 1, photcode[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[i]) x = RELPHOT_GRID_X - X - 1;
+	y = Y;
+	if (camera.Fy[i]) y = RELPHOT_GRID_Y - Y - 1;
+	      
+	/* coordinates in the grid */
+	ix = x + camera.Ox[i]*RELPHOT_GRID_X;
+	iy = y + camera.Oy[i]*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);
+  }
+# endif
+
+}
+
+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;
+    }
+  }
+}
Index: unk/Ohana/src/relphot/src/GridOps.v1.c
===================================================================
--- /trunk/Ohana/src/relphot/src/GridOps.v1.c	(revision 9632)
+++ 	(revision )
@@ -1,427 +1,0 @@
-# include "relphot.h"
-
-static int     Ngrid;
-static short   *gridM;
-static short   *gridS;
-static int     *gridN;
-static int      gridX;
-static int      gridY;
-
-static int    **bin;
-static int    **Xmeas;
-static int    **Ymeas;
-
-static int    **clist;
-static int    **mlist;
-static int     *Nlist;
-static int     *NLIST;
-
-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, short, Ngrid);
-  ALLOCATE (gridS, short, Ngrid);
-  ALLOCATE (gridN, int,   Ngrid);
-  bzero (gridM, Ngrid*sizeof(short));
-  bzero (gridS, Ngrid*sizeof(short));
-  bzero (gridN, Ngrid*sizeof(int));
-
-}
-
-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, int, Ngrid);
-  ALLOCATE (NLIST, int, Ngrid);
-  ALLOCATE (clist, int *, Ngrid);
-  ALLOCATE (mlist, int *, Ngrid);
-
-  for (i = 0; i < Ngrid; i++) {
-    Nlist[i] = 0;
-    NLIST[i] = 100;
-    ALLOCATE (clist[i], int, NLIST[i]);
-    ALLOCATE (mlist[i], int, 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 (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;
-
-  i = ix + iy*gridX;
-
-  bin[cat][meas] = i;
-  Xmeas[cat][meas] = X;
-  Ymeas[cat][meas] = Y;
-  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);
-
- escape:
-  fprintf (stderr, "error: star out of grid\n");
-  exit (1);
-}
-
-short getMgrid (int meas, int cat) {
-
-  int i;
-  short value;
-
-  if (!USE_GRID) return (0);
-  i = bin[cat][meas];
-  if (i == -1) return (NO_MAG);
-
-  value = gridM[i];
-  return (value);
-}
-
-/* determine Mgrid values for all grid bins */
-void setMgrid (Catalog *catalog) {
-
-  int i, j, m, c, n, N, Nmax;
-  double *list, *dlist;
-  short Msys, Mrel, Mcal, Mmos;
-  StatType stats;
-  
-  if (!USE_GRID) return;
-
-  Nmax = Nlist[0];
-  for (i = 0; i < Ngrid; i++) {
-    Nmax = MAX (Nmax, Nlist[i]);
-  }
-  ALLOCATE (list, double, Nmax);
-  ALLOCATE (dlist, double, Nmax);
-
-  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].measure[m].flags & MEAS_BAD) continue;
-      if ((Mcal  = getMcal  (m, c)) == NO_MAG) continue;
-      if ((Mmos  = getMmos  (m, c)) == NO_MAG) continue;
-      if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
-      
-      n = catalog[c].measure[m].averef;
-      Msys = iPhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
-      list[N] = Msys - Mrel - Mcal - Mmos;
-      dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
-      N++;
-    }
-
-    liststats (list, dlist, N, &stats);
-    gridM[i] = stats.mean;
-    gridS[i] = stats.sigma;
-    gridN[i] = N;
-  }
-  free (list);
-  free (dlist);
-}
-
-void plot_grid (Catalog *catalog) {
-
-  int i, j, m, c, n, N;
-  short Msys, Mrel, Mcal, Mmos;
-  double *xlist, *Mlist, *dlist, *ylist;
-  Graphdata graphdata;
-
-  if (!USE_GRID) return;
-
-  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);
-
-  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].measure[m].flags & MEAS_BAD) continue;
-      if ((Mcal  = getMcal  (m, c)) == NO_MAG) continue;
-      if ((Mmos  = getMmos  (m, c)) == NO_MAG) continue;
-      if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
-
-      n = catalog[c].measure[m].averef;
-      Msys = iPhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
-
-      xlist[N] = Xmeas[c][m];
-      ylist[N] = Ymeas[c][m];
-      Mlist[N] = Msys - Mrel - Mcal - Mmos;
-      dlist[N] = Msys - Mrel - Mcal - Mmos - gridM[i];
-      Mlist[N] *= 0.001;
-      dlist[N] *= 0.001;
-      N++;
-    }
-  }
-
-  plot_defaults (&graphdata);
-  graphdata.ymin = PlotdMmin;
-  graphdata.ymax = PlotdMmax;
-  plot_list (&graphdata, xlist, Mlist, N, "X vs dM raw", "XdM.png");
-  plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "XdMf.png");
-  plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "YdMf.png");
-
-  plot_defaults (&graphdata);
-  plot_list (&graphdata, xlist, ylist, N, "X vs Y", "XY.png");
-
-  free (ylist);
-  free (xlist);
-  free (Mlist);
-  free (dlist);
-
-}
-
-int *SelectRefMosaic (Mosaic **refmosaic, int *Nimage);
-Image *getimage (int N);
-
-void dump_grid () { 
-
-  int i, j, Nx, Ny, Nimage;
-  int *imlist;
-  char *p, *pname, ccdname[80];
-  FILE *f;
-  Header header, theader;
-  Matrix matrix;
-  Mosaic *refmosaic;
-  Image *image;
-    
-  /* select reference mosaic image */
-  imlist = SelectRefMosaic (&refmosaic, &Nimage);
-
-  /* we are writing to this file */
-  f = fopen ("mosaic.fits", "w");
-  if (f == (FILE *) NULL) { 
-    fprintf (stderr, "cannot open %s for output\n", "mosaic.fits");
-    return;
-  }
-
-  /* create empty phu */
-  gfits_init_header (&header);
-  header.extend = TRUE;
-  gfits_create_header (&header);
-  gfits_create_matrix (&header, &matrix);
-  gfits_modify (&header, "NEXTEND", "%d", 1, Nimage + 3);
-  gfits_modify (&header, "FILTER", "%s", 1, photcode[0].name);
-  gfits_modify (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis");
-  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, photcode[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, photcode[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, photcode[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 */
-  for (i = 0; i < Nimage; i++) {
-    image = getimage (imlist[i]);
-    pname = GetPhotcodeNamebyCode (image[0].source);
-
-    /* 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);
-  }
-}
-
-InterpolateGrid (float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords) {
-
-  int i, j;
-  double x, y, r, d, X, Y, dx, dy, dVx, dVy, dVxy;
-  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;
-    }
-  }
-}
Index: unk/Ohana/src/relphot/src/GridOps.v2.c
===================================================================
--- /trunk/Ohana/src/relphot/src/GridOps.v2.c	(revision 9632)
+++ 	(revision )
@@ -1,577 +1,0 @@
-# include "relphot.h"
-
-static int     Ngrid;
-static float   *gridM;
-static float   *gridS;
-static int     *gridN;
-static int      gridX;
-static int      gridY;
-
-static int    **bin;
-static int    **Xmeas;
-static int    **Ymeas;
-
-static int    **clist;
-static int    **mlist;
-static int     *Nlist;
-static int     *NLIST;
-
-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 */
-  char **ccdname;
-} camera;    
-
-# 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;
-  char field[64], *config, line[256];
-
-  /* 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);
-
-  ALLOCATE (camera.Fx, int, camera.Nchip);
-  ALLOCATE (camera.Fy, int, camera.Nchip);
-  ALLOCATE (camera.Ox, int, camera.Nchip);
-  ALLOCATE (camera.Oy, int, camera.Nchip);
-  ALLOCATE (camera.ccdname, char *, camera.Nchip);
-
-  /* load per-chip parameters */
-  for (i = 0; i < camera.Nchip; i++) {
-    ALLOCATE (camera.ccdname[i], char, 256);
-    sprintf (field, "CCD.%d", i);
-    ScanConfig (config, field, "%s", 1, line);
-    sscanf (line, "%s %d %d %d %d", camera.ccdname[i], &camera.Ox[i], &camera.Oy[i], &camera.Fx[i], &camera.Fy[i]);
-  }
-  free (config);
-
-  /* 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);
-  bzero (gridM, Ngrid*sizeof(float));
-  bzero (gridS, Ngrid*sizeof(float));
-  bzero (gridN, Ngrid*sizeof(int));
-
-}
-# endif
-
-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, int, Ngrid);
-  ALLOCATE (NLIST, int, Ngrid);
-  ALLOCATE (clist, int *, Ngrid);
-  ALLOCATE (mlist, int *, Ngrid);
-
-  for (i = 0; i < Ngrid; i++) {
-    Nlist[i] = 0;
-    NLIST[i] = 100;
-    ALLOCATE (clist[i], int, NLIST[i]);
-    ALLOCATE (mlist[i], int, 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);
-}
-
-# 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;
-
-  i = ix + iy*gridX;
-
-  bin[cat][meas] = i;
-  Xmeas[cat][meas] = X;
-  Ymeas[cat][meas] = Y;
-  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);
-
-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) {
-
-  int i;
-  float value;
-
-  if (!USE_GRID) return (0);
-  i = bin[cat][meas];
-  if (i == -1) return (NO_MAG);
-
-  value = gridM[i];
-  return (value);
-}
-
-/* determine Mgrid values for all grid bins */
-void setMgrid (Catalog *catalog) {
-
-  int i, j, m, c, n, N, Nmax;
-  double *list, *dlist;
-  float Msys, Mrel, Mcal, Mmos;
-  StatType stats;
-  
-  if (!USE_GRID) return;
-
-  Nmax = Nlist[0];
-  for (i = 0; i < Ngrid; i++) {
-    Nmax = MAX (Nmax, Nlist[i]);
-  }
-  ALLOCATE (list, double, Nmax);
-  ALLOCATE (dlist, double, Nmax);
-
-  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].measure[m].flags & MEAS_BAD) continue;
-      if ((Mcal  = getMcal  (m, c)) == NO_MAG) continue;
-      if ((Mmos  = getMmos  (m, c)) == NO_MAG) continue;
-      if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
-      
-      n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
-      list[N] = Msys - Mrel - Mcal - Mmos;
-      dlist[N] = MAX (catalog[c].measure[m].dM_PS, MIN_ERROR);
-      N++;
-    }
-
-    liststats (list, dlist, N, &stats);
-    gridM[i] = stats.mean;
-    gridS[i] = stats.sigma;
-    gridN[i] = N;
-  }
-  free (list);
-  free (dlist);
-}
-
-void plot_grid (Catalog *catalog) {
-
-  int i, j, m, c, n, N, Narea;
-  float Msys, Mrel, Mcal, Mmos;
-  double *xlist, *Mlist, *dlist, *ylist;
-  Graphdata graphdata;
-
-  if (!USE_GRID) return;
-
-  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].measure[m].flags & MEAS_BAD) {
-	Narea ++;
-	continue;
-      }
-      if ((Mcal  = getMcal  (m, c)) == NO_MAG) continue;
-      if ((Mmos  = getMmos  (m, c)) == NO_MAG) continue;
-      if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
-
-      n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
-
-      xlist[N] = Xmeas[c][m];
-      ylist[N] = Ymeas[c][m];
-      Mlist[N] = Msys - Mrel - Mcal - Mmos;
-      dlist[N] = Msys - Mrel - Mcal - Mmos - 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", "XdM.png");
-  plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "XdMf.png");
-  plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "YdMf.png");
-
-  plot_defaults (&graphdata);
-  plot_list (&graphdata, xlist, ylist, N, "X vs Y", "XY.png");
-
-  free (ylist);
-  free (xlist);
-  free (Mlist);
-  free (dlist);
-
-}
-
-void dump_grid () { 
-
-  int i, j, Nimage;
-  int *imlist;
-  FILE *f;
-  Header header, theader;
-  Matrix matrix;
-  Mosaic *refmosaic;
-    
-  /* select reference mosaic image */
-  imlist = SelectRefMosaic (&refmosaic, &Nimage);
-
-  /* we are writing to this file */
-  f = fopen ("mosaic.fits", "w");
-  if (f == (FILE *) NULL) { 
-    fprintf (stderr, "cannot open %s for output\n", "mosaic.fits");
-    return;
-  }
-
-  /* create empty phu */
-  gfits_init_header (&header);
-  header.extend = TRUE;
-  gfits_create_header (&header);
-  gfits_create_matrix (&header, &matrix);
-  gfits_modify (&header, "NEXTEND", "%d", 1, Nimage + 3);
-  gfits_modify (&header, "FILTER", "%s", 1, photcode[0].name);
-  gfits_modify (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis");
-  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, photcode[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, photcode[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, photcode[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);
-
-# 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].source);
-
-    /* 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 */
-  for (i = 0; i < camera.Nchip; i++) {
-    int ix, iy, x, y, X, Y, bin;
-
-    gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[i]);
-    gfits_modify (&theader, "FILTER", "%s", 1, photcode[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[i]) x = RELPHOT_GRID_X - X - 1;
-	y = Y;
-	if (camera.Fy[i]) y = RELPHOT_GRID_Y - Y - 1;
-	      
-	/* coordinates in the grid */
-	ix = x + camera.Ox[i]*RELPHOT_GRID_X;
-	iy = y + camera.Oy[i]*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);
-  }
-# endif
-
-}
-
-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;
-    }
-  }
-}
Index: /trunk/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- /trunk/Ohana/src/relphot/src/ImageOps.c	(revision 9632)
+++ /trunk/Ohana/src/relphot/src/ImageOps.c	(revision 9633)
@@ -142,5 +142,5 @@
 
   i = bin[cat][meas];
-  if (i == -1) return (NO_MAG);
+  if (i == -1) return (NO_IMAGE);
 
   if (image[i].code & IMAGE_BAD)  return (NO_MAG);  
@@ -194,5 +194,7 @@
       
       if (catalog[c].measure[m].flags & MEAS_BAD) continue;
-      if ((Mmos  = getMmos  (m, c)) == NO_MAG) continue;
+      Mmos  = getMmos  (m, c);
+      if (Mmos == NO_MAG) continue;
+      if (Mmos == NO_IMAGE) continue;
       if ((Mgrid = getMgrid (m, c)) == NO_MAG) continue;
       if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
@@ -334,4 +336,5 @@
   int i, j, m, c, n, N;
   double *list, *dlist;
+  float Mcal, Mmos;
   StatType stats;
 
@@ -352,6 +355,10 @@
       c = clist[i][j];
 
-      if (getMcal  (m, c) == NO_MAG) continue;
-      if (getMmos  (m, c) == NO_MAG) continue;
+      Mcal = getMcal  (m, c);
+      if (Mcal == NO_MAG) continue;
+      if (Mcal == NO_IMAGE) continue;
+      Mmos = getMmos  (m, c);
+      if (Mmos == NO_MAG) continue;
+      if (Mmos == NO_IMAGE) continue;
       if (getMgrid (m, c) == NO_MAG) continue;
       N++;
Index: /trunk/Ohana/src/relphot/src/MosaicOps.c
===================================================================
--- /trunk/Ohana/src/relphot/src/MosaicOps.c	(revision 9632)
+++ /trunk/Ohana/src/relphot/src/MosaicOps.c	(revision 9633)
@@ -299,5 +299,5 @@
   if (!MOSAICNAME[0]) return (0);
   i = bin[cat][meas];
-  if (i == -1) return (NO_MAG);
+  if (i == -1) return (NO_IMAGE);
 
   if (mosaic[i].code & IMAGE_BAD) return (NO_MAG);  
@@ -345,6 +345,8 @@
       c = clist[i][j];
       
+      Mcal  = getMcal  (m, c);
+      if (Mcal == NO_MAG) continue;
+      if (Mcal == NO_IMAGE) continue;
       if (catalog[c].measure[m].flags & MEAS_BAD) continue;
-      if ((Mcal  = getMcal  (m, c)) == NO_MAG) continue;
       if ((Mgrid = getMgrid (m, c)) == NO_MAG) continue;
       if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
@@ -447,4 +449,5 @@
   int i, j, m, c, n, N;
   double *list, *dlist;
+  float Mcal;
   StatType stats;
 
@@ -465,5 +468,7 @@
       c = clist[i][j];
 
-      if (getMcal  (m, c) == NO_MAG) continue;
+      Mcal = getMcal  (m, c);
+      if (Mcal == NO_MAG) continue;
+      if (Mcal == NO_IMAGE) continue;
       if (getMgrid (m, c) == NO_MAG) continue;
       if (getMrel  (catalog, m, c) == NO_MAG) continue;
Index: /trunk/Ohana/src/relphot/src/StarOps.c
===================================================================
--- /trunk/Ohana/src/relphot/src/StarOps.c	(revision 9632)
+++ /trunk/Ohana/src/relphot/src/StarOps.c	(revision 9633)
@@ -45,10 +45,21 @@
       m = catalog[i].average[j].offset;
 
+      // XXX allow REF stars to be included in the calculation
+      // this should be optionally set, and should allow for 
+      // REF stars to be downweighted by more than their reported
+      // errors.  how such info is carried is unclear...
       N = 0;
       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
 	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
-	if ((Mcal  = getMcal  (m, i)) == NO_MAG) continue;
-	if ((Mmos  = getMmos  (m, i)) == NO_MAG) continue;
-	if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
+	Mcal  = getMcal  (m, i);
+	if (Mcal == NO_MAG) continue;
+	if (Mcal == NO_IMAGE) {
+	  Mcal = Mmos = Mgrid = 0;
+	} else {
+	  Mmos  = getMmos  (m, i);
+	  if (Mmos == NO_MAG) continue;
+	  if (Mmos == NO_IMAGE) continue;
+	  if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
+	}
 
 	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
@@ -106,7 +117,14 @@
       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
 	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
-	if ((Mcal  = getMcal  (m, i)) == NO_MAG) continue;
-	if ((Mmos  = getMmos  (m, i)) == NO_MAG) continue;
-	if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
+	Mcal  = getMcal  (m, i);
+	if (Mcal == NO_MAG) continue;
+	if (Mcal == NO_IMAGE) {
+	  Mcal = Mmos = Mgrid = 0;
+	} else {
+	  Mmos  = getMmos  (m, i);
+	  if (Mmos == NO_MAG) continue;
+	  if (Mmos == NO_IMAGE) continue;
+	  if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
+	}
 
 	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
@@ -153,6 +171,10 @@
       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
 	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
-	if ((Mcal  = getMcal  (m, i)) == NO_MAG) continue;
-	if ((Mmos  = getMmos  (m, i)) == NO_MAG) continue;
+	Mcal  = getMcal  (m, i);
+	if (Mcal == NO_MAG) continue;
+	if (Mcal == NO_IMAGE) continue;
+	Mmos  = getMmos  (m, i);
+	if (Mmos == NO_MAG) continue;
+	if (Mmos == NO_IMAGE) continue;
 	if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
 	catalog[i].measure[m].Mcal_PS = Mcal + Mmos + Mgrid;
@@ -263,6 +285,10 @@
       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
 	/* if (catalog[i].measure[m].flags & MEAS_BAD) continue; */
-	if ((Mcal  = getMcal  (m, i)) == NO_MAG) continue;
-	if ((Mmos  = getMmos  (m, i)) == NO_MAG) continue;
+	Mcal  = getMcal  (m, i);
+	if (Mcal == NO_MAG) continue;
+	if (Mcal == NO_IMAGE) continue;
+	Mmos  = getMmos  (m, i);
+	if (Mmos == NO_MAG) continue;
+	if (Mmos == NO_IMAGE) continue;
 	if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
 
@@ -295,6 +321,10 @@
       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
 	/* if (catalog[i].measure[m].flags & MEAS_BAD) continue; */
-	if ((Mcal  = getMcal  (m, i)) == NO_MAG) continue;
-	if ((Mmos  = getMmos  (m, i)) == NO_MAG) continue;
+	Mcal  = getMcal  (m, i);
+	if (Mcal == NO_MAG) continue;
+	if (Mcal == NO_IMAGE) continue;
+	Mmos  = getMmos  (m, i);
+	if (Mmos == NO_MAG) continue;
+	if (Mmos == NO_IMAGE) continue;
 	if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
 
@@ -328,4 +358,5 @@
   int i, j, k, m, n, N, Ntot;
   double *list, *dlist;
+  float Mcal, Mmos;
   StatType stats;
 
@@ -348,6 +379,10 @@
       N = 0;
       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
-	if (getMcal  (m, i) == NO_MAG) continue;
-	if (getMmos  (m, i) == NO_MAG) continue;
+	Mcal = getMcal  (m, i);
+	if (Mcal == NO_MAG) continue;
+	if (Mcal == NO_IMAGE) continue;
+	Mmos = getMmos  (m, i);
+	if (Mmos == NO_MAG) continue;
+	if (Mmos == NO_IMAGE) continue;
 	if (getMgrid (m, i) == NO_MAG) continue;
 	N++;
Index: /trunk/Ohana/src/relphot/src/args.c
===================================================================
--- /trunk/Ohana/src/relphot/src/args.c	(revision 9632)
+++ /trunk/Ohana/src/relphot/src/args.c	(revision 9633)
@@ -171,4 +171,12 @@
     remove_argument (N, &argc, argv);
     ImagSelect = TRUE;
+  }
+
+  DophotSelect = FALSE;
+  if ((N = get_argument (argc, argv, "-dophot"))) {
+    remove_argument (N, &argc, argv);
+    DophotValue = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    DophotSelect = TRUE;
   }
 
Index: /trunk/Ohana/src/relphot/src/bcatalog.c
===================================================================
--- /trunk/Ohana/src/relphot/src/bcatalog.c	(revision 9632)
+++ /trunk/Ohana/src/relphot/src/bcatalog.c	(revision 9633)
@@ -18,11 +18,4 @@
   for (i = 0; i < catalog[0].Naverage; i++) {
     if (catalog[0].average[i].Nm < 2) continue; 
-
-    /* XXX this limitation is absurd
-    if (catalog[0].average[i].R < fullregion[0].RA[0]) continue;
-    if (catalog[0].average[i].R > fullregion[0].RA[1]) continue;
-    if (catalog[0].average[i].D < fullregion[0].DEC[0]) continue;
-    if (catalog[0].average[i].D > fullregion[0].DEC[1]) continue;
-    */
 
     /* start with all stars good */
@@ -59,5 +52,6 @@
 
       /* select measurements by quality */
-      if (catalog[0].measure[offset].dophot != 1) continue;
+      // XXX ignore this criterion for REF measurements?
+      if (DophotSelect && (catalog[0].measure[offset].dophot != DophotValue)) continue;
 
       /* select measurements by mag limit */
@@ -66,5 +60,5 @@
 
       /* select measurements by measurement error */
-      if (catalog[0].measure[offset].dM_PS > SIGMA_LIM) continue;
+      if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM_PS > SIGMA_LIM)) continue;
 
       /* select measurements by mag limit */
@@ -92,8 +86,12 @@
       }
     }
+
+    // XXXX test : what checks do I need to make elsewhere to avoid problems here?
+    # if 1
     if (Nm < 2) { /* enough measurements in band? */
       Nmeasure -= Nm;
       continue; 
     }
+    # endif
     subcatalog[0].average[Naverage].Nm = Nm;
     Naverage ++;
Index: /trunk/Ohana/src/relphot/src/load_catalogs.c
===================================================================
--- /trunk/Ohana/src/relphot/src/load_catalogs.c	(revision 9632)
+++ /trunk/Ohana/src/relphot/src/load_catalogs.c	(revision 9633)
@@ -17,5 +17,6 @@
     tcatalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
     tcatalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
-    tcatalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point
+    tcatalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF;    // don't need to load all data at this point
+    tcatalog.Nsecfilt  = GetPhotcodeNsecfilt ();               // set the desired number in case we need to create the catalog
 
     if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE, "r")) {
@@ -37,6 +38,5 @@
   }
   if (Nstar < 2) { 
-    fprintf (stderr, "insufficient stars %d\n", Nstar);
-    exit (0);
+    fprintf (stderr, "warning: insufficient stars %d\n", Nstar);
   }
 
Index: /trunk/Ohana/src/relphot/src/plot_scatter.c
===================================================================
--- /trunk/Ohana/src/relphot/src/plot_scatter.c	(revision 9632)
+++ /trunk/Ohana/src/relphot/src/plot_scatter.c	(revision 9633)
@@ -28,6 +28,10 @@
       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
 	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
-	if ((Mcal  = getMcal  (m, i)) == NO_MAG) continue;
-	if ((Mmos  = getMmos  (m, i)) == NO_MAG) continue;
+	Mcal  = getMcal  (m, i);
+	if (Mcal == NO_MAG) continue;
+	if (Mcal == NO_IMAGE) continue;
+	Mmos  = getMmos  (m, i);
+	if (Mmos == NO_MAG) continue;
+	if (Mmos == NO_IMAGE) continue;
 	if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
 
Index: /trunk/Ohana/src/relphot/src/reload_catalogs.c
===================================================================
--- /trunk/Ohana/src/relphot/src/reload_catalogs.c	(revision 9632)
+++ /trunk/Ohana/src/relphot/src/reload_catalogs.c	(revision 9633)
@@ -24,4 +24,5 @@
     catalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
     catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
+    catalog.Nsecfilt  = GetPhotcodeNsecfilt ();               // set the desired number in case we need to create the catalog
 
     if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
Index: /trunk/Ohana/src/relphot/src/relphot.c
===================================================================
--- /trunk/Ohana/src/relphot/src/relphot.c	(revision 9632)
+++ /trunk/Ohana/src/relphot/src/relphot.c	(revision 9633)
@@ -14,10 +14,22 @@
   /* register database handle with shutdown procedure */
   set_db (&db);
+  db.mode   = dvo_catalog_catmode (CATMODE);
+  db.format = dvo_catalog_catformat (CATFORMAT);
 
   /* lock and load the image db table */
   status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT));
-  if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
-  if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
-  if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
+  if (!status && UPDATE) {
+    fprintf (stderr, "error\n");
+    Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
+  }
+  // if the file is missing, db.dbstate will have a value of either:
+  // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE)
+  if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) {
+    // XXX get ZERO_POINT from config
+    dvo_image_create (&db, 25.0);
+    // Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
+  } else {
+    if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
+  }
 
   /* load regions and images based on specified sky patch */
Index: /trunk/Ohana/src/relphot/src/setMrelFinal.c
===================================================================
--- /trunk/Ohana/src/relphot/src/setMrelFinal.c	(revision 9632)
+++ /trunk/Ohana/src/relphot/src/setMrelFinal.c	(revision 9633)
@@ -110,4 +110,5 @@
 
       /* skip measurements from BAD images and mosaics */
+      /* do NOT skip measurements without a matching image */
       if ((getMcal  (m, 0)) == NO_MAG) goto skip;
       if ((getMmos  (m, 0)) == NO_MAG) goto skip;