# 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;
    }
  }
}
