# include "relphot.h"

enum {
    GRID_FITTED,
    GRID_FROZEN,
    GRID_REFERENCE,
};

static int     Ngrid;
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;

static int    **bin;
static int    **Xmeas;
static int    **Ymeas;

static int    **clist; // link from measurement on a cell to catalog containing measurement
static int    **mlist; // link from measurement on a cell to measurement in a catalog
static int     *Nlist; // list of measurements for each grid cell
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 */
  int *valid;
  char **ccdname;
} camera;    

static int *ccdnum;
static char *config;

void initGrid (int dX, int dY) {

  int i, N, ccdnum_max, refX, refY, refBin;
  char *p, field[64], line[256];
  int *Fx, *Fy;  /* chip flip */
  int *Ox, *Oy;  /* chip offset */
  char **ccdname;

  /* load camera config file */
  config = LoadConfigFile (CameraConfig);
  if (config == (char *) NULL) {
    fprintf (stderr, "ERROR: can't find camera config file %s\n", CameraConfig);
    exit (1);
  }

  /* load basic mosaic parameters */
  ScanConfig (config, "NCCD", "%d", 1, &camera.Nchip);
  ScanConfig (config, "MOSAIC_X", "%d", 1, &camera.Mx);
  ScanConfig (config, "MOSAIC_Y", "%d", 1, &camera.My);
  ScanConfig (config, "NAXIS1", "%d", 1, &camera.Nx);
  ScanConfig (config, "NAXIS2", "%d", 1, &camera.Ny);

  ScanConfig (config, "REFCELL.X", "%d", 1, &refX);
  ScanConfig (config, "REFCELL.Y", "%d", 1, &refY);

  // temporary storage
  ALLOCATE (Fx,      int,    camera.Nchip);
  ALLOCATE (Fy,      int,    camera.Nchip);
  ALLOCATE (Ox,      int,    camera.Nchip);
  ALLOCATE (Oy,      int,    camera.Nchip);
  ALLOCATE (ccdname, char *, camera.Nchip);
  ALLOCATE (ccdnum,  int,    camera.Nchip);

  /* load per-chip parameters */
  for (i = 0; i < camera.Nchip; i++) {
    ALLOCATE (ccdname[i], char, 256);
    sprintf (field, "CCD.%d", i);
    ScanConfig (config, field, "%s", 1, line);
    // XXX get error status!

    sscanf (line, "%s %d %d %d %d", ccdname[i], &Ox[i], &Oy[i], &Fx[i], &Fy[i]);
    
    p = ccdname[i];
    while (!isdigit(*p) && *p) p++;
    if (*p == 0) continue;
    ccdnum[i] = atoi (p);
  }

  /* we now have the parameters loaded into a minimal length list; reshuffle to the a list of length 0 - MAX(ccdnum) */
  ccdnum_max = 0;
  for (i = 0; i < camera.Nchip; i++) {
    ccdnum_max = MAX(ccdnum_max, ccdnum[i]);
  }
  ccdnum_max ++;

  if (ccdnum_max < camera.Nchip) {
    fprintf (stderr, "problem with camera config: duplicate ccd IDs\n");
    exit (1);
  }
  
  if (ccdnum_max > 0x1000) {
    fprintf (stderr, "problem with camera config: absurd max ccd ID number %d\n", ccdnum_max);
    exit (1);
  }

  ALLOCATE (camera.Fx, 	    int,    ccdnum_max);
  ALLOCATE (camera.Fy, 	    int,    ccdnum_max);
  ALLOCATE (camera.Ox, 	    int,    ccdnum_max);
  ALLOCATE (camera.Oy, 	    int,    ccdnum_max);
  ALLOCATE (camera.valid,   int,    ccdnum_max);
  ALLOCATE (camera.ccdname, char *, ccdnum_max);

  for (i = 0; i < ccdnum_max; i++) {
    camera.valid[i] = FALSE;
  }

  for (i = 0; i < camera.Nchip; i++) {
    N = ccdnum[i];
    camera.Fx[N] = Fx[i];
    camera.Fy[N] = Fy[i];
    camera.Ox[N] = Ox[i];
    camera.Oy[N] = Oy[i];
    camera.ccdname[N] = ccdname[i];
    camera.valid[N] = TRUE;
  }

  /* define mosaic 2d correction grid: 
   * GRID_X is the number of grid pixels per chip in the X direction */
  gridX = RELPHOT_GRID_X * camera.Mx;
  gridY = RELPHOT_GRID_Y * camera.My;
  Ngrid = gridX * gridY;

  ALLOCATE (gridM, float, Ngrid);
  ALLOCATE (gridS, float, Ngrid);
  ALLOCATE (gridN, int,   Ngrid);
  ALLOCATE (gridV, int,   Ngrid);

  // the grid bins may have one of three possible states: fitted, frozen, reference
  // set the initial values to indicate that the bins are frozen 
  for (i = 0; i < Ngrid; i++) {
    gridM[i] = 0.0;
    gridS[i] = 0.0;
    gridN[i] = 0;
    gridV[i] = GRID_FROZEN;
  }

  // refBin is the index of the grid cell which is kept at 0.0 (all others are relative to this)
  refBin = refX + refY*gridX;
  gridV[refBin] = GRID_REFERENCE;
}

void initGridBins (Catalog *catalog, int Ncatalog) {

  int i, j;

  if (!USE_GRID) return;

  /* define cat,meas -> grid pointers */
  ALLOCATE (bin, int *, Ncatalog);
  ALLOCATE (Xmeas, int *, Ncatalog);
  ALLOCATE (Ymeas, int *, Ncatalog);
  for (i = 0; i < Ncatalog; i++) {
    ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1));
    ALLOCATE (Xmeas[i], int, MAX (catalog[i].Nmeasure, 1));
    ALLOCATE (Ymeas[i], int, MAX (catalog[i].Nmeasure, 1));
    for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
  }

  /* define grid -> cat,meas pointers */
  ALLOCATE (Nlist, 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 ccdnum) {

  int ix, iy, Cx, Cy, i;
  double x, y;

  /* X, Y are chip coords on chip ccdnum */

  /* normalize X & Y */
  x = X;
  if (camera.Fx[ccdnum]) x = camera.Nx - X;
  y = Y;
  if (camera.Fy[ccdnum]) y = camera.Ny - Y;

  /* grid coords on the chip */
  Cx = MIN (MAX ((x / camera.Nx) * RELPHOT_GRID_X, 0), RELPHOT_GRID_X - 1);
  Cy = MIN (MAX ((y / camera.Ny) * RELPHOT_GRID_Y, 0), RELPHOT_GRID_Y - 1);
  
  /* coordinates in the grid */
  ix = Cx + camera.Ox[ccdnum]*RELPHOT_GRID_X;
  iy = Cy + camera.Oy[ccdnum]*RELPHOT_GRID_Y;

  i = ix + iy*gridX;

  bin[cat][meas] = i;
  Xmeas[cat][meas] = x + camera.Ox[ccdnum]*camera.Nx;
  Ymeas[cat][meas] = y + camera.Oy[ccdnum]*camera.Ny;
  clist[i][Nlist[i]] = cat;
  mlist[i][Nlist[i]] = meas;

  // for the moment, add up the total grid count
  gridN[i]++;

  Nlist[i] ++;
  if (Nlist[i] == NLIST[i]) {
    NLIST[i] += 100;
    REALLOCATE (clist[i], int, NLIST[i]);
    REALLOCATE (mlist[i], int, NLIST[i]);
  }	
  return (TRUE);

  fprintf (stderr, "error: star out of grid\n");
  exit (1);
}

int showGridCount() {

    int ix, iy, i;

    for (iy = 0; iy < gridY; iy++) {
	for (ix = 0; ix < gridX; ix++) {
	    i = ix + iy*gridX;
	    fprintf (stderr, "%3d  ", gridN[i]);
	}
	fprintf (stderr, "\n");
    }
    return (TRUE);
}

float getMgrid (int meas, int cat) {

  int i;
  float value;

  if (!USE_GRID) return (0);
  i = bin[cat][meas];
  if (i == -1) return (NAN);

  // during the grid annealing process, we skip over grid cells until they have enough
  // valid stars to be fitted
  if (gridV[i] == GRID_FROZEN) {
    return (NAN);
  }

  value = gridM[i];
  return (value);
}

/* direct (non-iterative) solution for Mgrid values for all grid bins */
void setMgridDirect (Catalog *catalog, int Ncatalog) {

  int **gotstar, **gridmeas;
  int i, j, k, N, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys;
  double **A, **B, *Mjx, *Wjx;
  float Msys, Mcal, Mmos, Merr, Wsys;
  double Mj, Wj;
  
  if (!USE_GRID) return;

  ALLOCATE (A, double *, Ngrid);
  ALLOCATE (B, double *, Ngrid);
  for (i = 0; i < Ngrid; i++) {
    ALLOCATE (A[i], double, Ngrid);
    ALLOCATE (B[i], double, 1);
    memset (A[i], 0, Ngrid*sizeof(double));
    memset (B[i], 0, sizeof(double));
  }

  Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = 0;

  ALLOCATE (gotstar, int *, Ncatalog);
  for (i = 0; i < Ncatalog; i++) {
    ALLOCATE (gotstar[i], int, catalog[i].Naverage);
  }  

  // set up gridmeas table : grid index for each measurement
  ALLOCATE (gridmeas, int *, Ncatalog);
  for (i = 0; i < Ncatalog; i++) {
    ALLOCATE (gridmeas[i], int, catalog[i].Nmeasure);
    for (j = 0; j < catalog[i].Nmeasure; j++) {
      gridmeas[i][j] = -1;
    }
  }  
  for (i = 0; i < Ngrid; i++) {
    for (j = 0; j < Nlist[i]; j++) {
      int m, c;
      m = mlist[i][j];
      c = clist[i][j];
      gridmeas[c][m] = i;
    }
  }

  // as we loop over the grid cells, we need to accumulate the values of Wjx for each star
  ALLOCATE (Wjx, double, Ngrid);
  ALLOCATE (Mjx, double, Ngrid);

  // accumulate the elements of the matrix equation.  We have an equation of the form: Ax = B
  // where x is the vector of grid cell values G_x (x = 0 - Ngrid), A is an Ngrid x Ngrid matrix,
  // and B is an Ngrid vector.  For a cell A(x,y), we need the following elements from each
  // star which touches grid cell (x):
  // 
  // Mj  : sum over all measurements of Msys / dMsys^2 
  // Wj  : sum over all measurements of 1.0 / dMsys^2 
  // Mjx : sum over all measurements which touch cell (x) of Msys / dMsys^2
  // Wjx : sum over all measurements which touch cell (x) of 1.0 / dMsys^2

  // Mjx and Wjx can be calculated by summing over all measurements which touch the cell 
  // Mj requires looping over stars which touch (x)

  // this is tricky because we need to know both the measurements which touch a cell
  // and the stars for which any measurement touches a cell.  we need to accumulate
  // sums for each star which touches as cell on both bases.

  for (i = 0; i < Ngrid; i++) {
    
    for (j = 0; j < Ncatalog; j++) {
      memset (gotstar[j], 0, catalog[j].Naverage*sizeof(int));
    }  

    // we are looping over the stars, but doing so by looping over the set of measurements:
    // every star which touches this grid cell has a measurement in Nlist[i]
    for (j = 0; j < Nlist[i]; j++) {
      
      int mx, c, n, m0, Npts;

      mx = mlist[i][j];
      c  = clist[i][j];
      n  = catalog[c].measure[mx].averef;
      
      // if we have already visited this star, skip the stuff below
      if (gotstar[c][n]) continue;
      gotstar[c][n] = TRUE;

      // skip stars marked as BAD
      if (catalog[c].average[n].code & STAR_BAD) {
	Nrel ++;
	continue;
      }

      m0 = catalog[c].average[n].measureOffset;

      // we accumuate an entry for each cell
      memset (Wjx, 0, Ngrid*sizeof(double));
      memset (Mjx, 0, Ngrid*sizeof(double));
      Npts = Mj = Wj = 0.0;

      // if we have not yet visited this star, accumulate the Mj, Wj entries
      for (k = 0; k < catalog[c].average[n].Nmeasure; k++) {

	int m, Ng;

	m = m0 + k;

	// skip measurements marked as BAD
	if (catalog[c].measure[m].dbFlags & MEAS_BAD) {
	  Nbad ++;
	  continue;
	}

	// skip images marked as BAD
	Mcal = getMcal  (m, c);
	if (isnan(Mcal)) {
	  Ncal ++;
	  continue;
	}

	// skip mosaics marked as BAD
	Mmos = getMmos  (m, c);
	if (isnan(Mmos)) {
	  Nmos ++;
	  continue;
	}

	// select the color- and airmass-corrected observed magnitude for this star
	// XXX need to be able to turn off the color-correction until initial average mags are found
	// Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
	Msys = PhotCat (&catalog[c].measure[m]);
	if (isnan(Msys)) {
	  Nsys++;
	  continue;
	}

	// mag-error for this measurement
	Merr =  MAX (catalog[c].measure[m].dM, MIN_ERROR);

	// Wsys = 1.0 / SQ(Merr);
	Wsys = 1.0;

	Ng = gridmeas[c][m];
	if (Ng == -1) continue;  // skip measurements which do not touch any cell

	Mj += Msys * Wsys;  // we are only including measurements touching this cell
	Wj += Wsys;  // we are only including measurements touching this cell
	Npts ++;
	Ngood ++;

	Mjx[Ng] += Msys * Wsys;  // we are only including measurements touching cell (x)
	Wjx[Ng] += Wsys;  // we are only including measurements touching cell (x)
      }

      // some stars will not have any valid measurements, skip these
      if (Npts == 0) continue;

      B[i][0] += Mj*Wjx[i]/Wj - Mjx[i];
      A[i][i] -= Wjx[i];
      for (k = 0; k < Ngrid; k++) {
	A[i][k] += Wjx[i]*Wjx[k]/Wj;
	// fprintf (stderr, "%3.0f ", Wjx[k]);
      }
    }
  }

  if (1) {

    FILE *f;
    Header theader;
    Matrix matrix;

    /* we are writing to this file */
    f = fopen ("matrix.fits", "w");
    if (f == (FILE *) NULL) { 
      fprintf (stderr, "cannot open matrix.fits for output\n");
      return;
    }

    /* save grid mag values */
    gfits_init_header (&theader);
    theader.Naxes = 2;
    theader.Naxis[0] = Ngrid;
    theader.Naxis[1] = Ngrid;
    theader.bitpix   = -32;
    gfits_create_Theader (&theader, "IMAGE");
    gfits_modify (&theader, "EXTNAME", "%s", 1, "MATRIX");
    gfits_create_matrix  (&theader, &matrix);
    for (i = 0; i < Ngrid; i++) {
      for (j = 0; j < Ngrid; j++) {
	gfits_set_matrix_value (&matrix, i, j, (double) A[i][j]);
      }
    }
    gfits_fwrite_header (f, &theader);
    gfits_fwrite_matrix (f, &matrix);
    gfits_free_matrix (&matrix);

    gfits_modify (&theader, "EXTNAME", "%s", 1, "TRPOSE");
    gfits_create_matrix  (&theader, &matrix);
    for (i = 0; i < Ngrid; i++) {
      for (j = 0; j < Ngrid; j++) {
	gfits_set_matrix_value (&matrix, i, j, (double) A[j][i]);
      }
    }
    gfits_fwrite_header (f, &theader);
    gfits_fwrite_matrix (f, &matrix);
    gfits_free_matrix (&matrix);
    fclose (f);
  }

  dgaussjordan (A, B, Ngrid, 1);

  fprintf (stderr, "grid cells fitted (Ngood: %d, Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Ngood, Nbad, Nmos, Ncal, Nrel, Nsys);

  for (i = 0; i < Ngrid; i++) {
    gridM[i] = B[i][0];
    gridS[i] = sqrt(A[i][i]);
    gridN[i] = N;
  }

  free (Wjx);
  free (Mjx);

  for (i = 0; i < Ngrid; i++) {
    free (A[i]);
  }
  free (A);
  free (B);

  for (i = 0; i < Ncatalog; i++) {
    free (gotstar[i]);
    free (gridmeas[i]);
  }  
  free (gotstar);
  free (gridmeas);
}

/* determine Mgrid values for all grid bins */
void setMgrid (Catalog *catalog) {

  int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ncal, Nrel, Nsys, Nfit;
  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);

  Nbad = Ncal = Nmos = Nrel = Nsys = Nfit = 0;

  for (i = 0; i < Ngrid; i++) {

    N = 0;
    for (j = 0; j < Nlist[i]; j++) {
      
      m = mlist[i][j];
      c = clist[i][j];
      
      if (catalog[c].measure[m].dbFlags & MEAS_BAD) {
	Nbad ++;
	continue;
      }
      Mcal = getMcal  (m, c);
      if (isnan(Mcal)) {
	Ncal ++;
	continue;
      }
      Mmos = getMmos  (m, c);
      if (isnan(Mmos)) {
	Nmos ++;
	continue;
      }
      Mrel  = getMrel  (catalog, m, c);
      if (isnan(Mrel)) {
	Nrel ++;
	continue;
      }
      
      n = catalog[c].measure[m].averef;
      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
      if (isnan(Msys)) {
	Nsys++;
	continue;
      }
      list[N] = Msys - Mrel - Mcal - Mmos;
      dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
      N++;
    }

    // the reference Cell is forced to have a value of 0.0, and is never changed to GRID_FITTED
    if (gridV[i] == GRID_REFERENCE) {
      gridM[i] = 0.0;
      gridS[i] = 0.0;
      gridN[i] = N;
      continue;
    }

    // until we have enough valid measurements on this grid cell, skip it
    if (N < GRID_TOOFEW) {
      gridV[i] = GRID_FROZEN;
      continue;
    }

    liststats (list, dlist, N, &stats);
    gridM[i] = stats.mean;
    gridS[i] = stats.sigma;
    gridN[i] = N;
    gridV[i] = GRID_FITTED;
    Nfit++;
  }

  fprintf (stderr, "%d of %d grid cells fitted (+ reference cell) (Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Nfit, Ngrid, Nbad, Nmos, Ncal, Nrel, Nsys);

  free (list);
  free (dlist);
}

void plot_grid (Catalog *catalog) {

  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].dbFlags & MEAS_BAD) {
	Narea ++;
	continue;
      }
      Mcal  = getMcal  (m, c);
      if (isnan(Mcal)) continue;
      Mmos  = getMmos  (m, c);
      if (isnan(Mmos)) continue;
      Mrel  = getMrel  (catalog, m, c);
      if (isnan(Mrel)) 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", "%s.XdM.png", OUTROOT);
  plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "%s.XdMf.png", OUTROOT);
  plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "%s.YdMf.png", OUTROOT);

  plot_defaults (&graphdata);
  plot_list (&graphdata, xlist, ylist, N, "X vs Y", "%s.XY.png", OUTROOT);

  free (ylist);
  free (xlist);
  free (Mlist);
  free (dlist);

}

void dump_grid () { 

  int i, j, Nimage, Nbytes, Nformat;
  int *imlist;
  FILE *f;
  Header header, theader;
  Matrix matrix;
  Mosaic *refmosaic;
  char *filename;
  char formatline[32], key[32], value[64];
    
  Nbytes = strlen (OUTROOT) + 6;
  ALLOCATE (filename, char, Nbytes);
  snprintf (filename, Nbytes, "%s.fits", OUTROOT);

  /* select reference mosaic image */
  imlist = SelectRefMosaic (&refmosaic, &Nimage);

  /* we are writing to this file */
  f = fopen (filename, "w");
  if (f == (FILE *) NULL) { 
    fprintf (stderr, "cannot open %s for output\n", filename);
    free (filename);
    return;
  }

  /* create empty phu */
  gfits_init_header (&header);
  header.extend = TRUE;
  gfits_create_header (&header);
  gfits_create_matrix (&header, &matrix);
  gfits_modify (&header, "NEXTEND", "%d", 1, Nimage + 3);
  gfits_modify (&header, "FILTER", "%s", 1, photcode[0].name);
  gfits_modify (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis");

  // we need to add lines to the PHU to identify the camera and format; these are used by the ipp config system
  // Note that config must have been loaded (and not freed) above.
  ScanConfig (config, "NFORMAT", "%d", 1, &Nformat);
  for (i = 1; i <= Nformat; i++) {
      ScanConfig (config, "FORMAT", "%s", i, formatline);
      sscanf (formatline, "%s %s", key, value);
      gfits_modify (&header, key, "%s", 1, value);
  }

  gfits_fwrite_header (f, &header);
  gfits_fwrite_matrix (f, &matrix);
  gfits_free_matrix (&matrix);

  /* save grid mag values */
  gfits_init_header (&theader);
  theader.Naxes = 2;
  theader.Naxis[0] = gridX;
  theader.Naxis[1] = gridY;
  theader.bitpix   = -32;
  gfits_create_Theader (&theader, "IMAGE");
  gfits_modify (&theader, "FILTER", "%s", 1, 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 */
  /* grid pixels are tied to detector pixels, but are flipped to match focal plane */
  for (i = 0; i < camera.Nchip; i++) {
    int N, ix, iy, x, y, X, Y, bin;

    N = ccdnum[i];

    gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[N]);
    gfits_modify (&theader, "FILTER", "%s", 1, 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[N]) x = RELPHOT_GRID_X - X - 1;
	y = Y;
	if (camera.Fy[N]) y = RELPHOT_GRID_Y - Y - 1;
	      
	/* coordinates in the grid */
	ix = x + camera.Ox[N]*RELPHOT_GRID_X;
	iy = y + camera.Oy[N]*RELPHOT_GRID_Y;
	      
	bin = ix + iy*gridX;
	gfits_set_matrix_value (&matrix, X, Y, (double) gridM[bin]);
      }
    }
    gfits_fwrite_header (f, &theader);
    gfits_fwrite_matrix (f, &matrix);
    gfits_free_matrix (&matrix);
  }

  free (filename);
}

void InterpolateGrid (float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords) {

  int i, j;
  double x, y, r, d, X, Y, dx, dy;
  double V00, V01, V10, V11;
  double wV00, wV01, wV10, wV11;
  double dV00, dV01, dV10, dV11;
  double v1, v2, value;
  int ix, iy, N;

  for (i = 0; i < Nx; i++) {
    for (j = 0; j < Ny; j++) {
      x = i * RELPHOT_GRID_BINNING / 2;
      y = j * RELPHOT_GRID_BINNING / 2;
      XY_to_RD (&r, &d, x, y, ccd);
      RD_to_XY (&X, &Y, r, d, gcoords);

      X = X / RELPHOT_GRID_BINNING;
      Y = Y / RELPHOT_GRID_BINNING;

      ix = (int) X;
      dx = X - ix;
      iy = (int) Y;
      dy = Y - iy;

      if (ix < 0) continue;
      if (iy < 0) continue;
      if (ix >= gridX) continue;
      if (iy >= gridY) continue;

      N = ix + iy*gridX;
      V00 = gridM[N];
      V10 = gridM[N + 1];
      V01 = gridM[N + gridX];
      V11 = gridM[N + gridX + 1];

      dV00 = gridS[N];
      dV10 = gridS[N + 1];
      dV01 = gridS[N + gridX];
      dV11 = gridS[N + gridX + 1];

      wV00 = (dV00 == 0) ? 0.0 : 1 / SQ(dV00);
      wV01 = (dV01 == 0) ? 0.0 : 1 / SQ(dV01);
      wV10 = (dV10 == 0) ? 0.0 : 1 / SQ(dV10);
      wV11 = (dV11 == 0) ? 0.0 : 1 / SQ(dV11);

      v1 = wV00*V00*(1 + dx*dy - dx - dy) +
	wV10*V10*(dx - dx*dy) +
	wV01*V01*(dy - dx*dy) +
	wV11*V11*(dx*dy);

      v2 = wV00*(1 + dx*dy - dx - dy) +
	wV10*(dx - dx*dy) +
	wV01*(dy - dx*dy) +
	wV11*(dx*dy);

      value = v1 / v2;
      buffer[j*Nx + i] = value;
    }
  }
}
