Index: unk/Ohana/src/relastro/src/GridOps.v1.c
===================================================================
--- /trunk/Ohana/src/relastro/src/GridOps.v1.c	(revision 12066)
+++ 	(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/relastro/src/GridOps.v2.c
===================================================================
--- /trunk/Ohana/src/relastro/src/GridOps.v2.c	(revision 12066)
+++ 	(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;
-    }
-  }
-}