# include "relphot.h"

static unsigned int *start;
static unsigned int *stop;
static int         **bin;

static int         **clist;
static int         **mlist;
static int          *Nlist;
static int          *NLIST;

static Image        *image;
static int          Nimage;

void initImages (Image *input, int N) {

  int i;

  image = input;
  Nimage = N;

  ALLOCATE (start,   unsigned, Nimage);
  ALLOCATE (stop,    unsigned, Nimage);

  for (i = 0; i < Nimage; i++) {
    start[i] = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
    stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
  }
}

void initImageBins (Catalog *catalog, int Ncatalog) {

  int i, j;

  ALLOCATE (bin, int *, Ncatalog);
  for (i = 0; i < Ncatalog; i++) {
    ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1));
    for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
  }

  ALLOCATE (Nlist, int, Nimage);
  ALLOCATE (NLIST, int, Nimage);
  ALLOCATE (clist, int *, Nimage);
  ALLOCATE (mlist, int *, Nimage);

  for (i = 0; i < Nimage; i++) {
    Nlist[i] = 0;
    NLIST[i] = 100;
    ALLOCATE (clist[i], int, NLIST[i]);
    ALLOCATE (mlist[i], int, NLIST[i]);
  }
}

void freeImageBins (int Ncatalog) {

  int i;

  for (i = 0; i < Ncatalog; i++) {
    free (bin[i]);
  }
  free (bin);
  for (i = 0; i < Nimage; i++) {
    free (clist[i]);
    free (mlist[i]);
  }
  free (clist);
  free (mlist);
}

/* **** use this to mimic the Mosaic version ***/
# if (0)
void initImageGrid (Image *image, int Nimage) {

  /* find max dR, dD range for all mosaics */
  /* define mosaic.coords to cover dR, dD */
  /* send results to initGridBins */

  int i, j, m, NX, NY;
  int dXmax, dYmax;
  double dS, dX, dY;
  double R, D, Rmin, Rmax, Dmin, Dmax;
  double Mcal, dMcal, Xm;

  dXmax = dYmax = 0.0;
  for (i = 0; i < Nmosaic; i++) {
    Dmin = Rmin =  360.0;
    Dmax = Rmax = -360.0;
    dS = 0.0;
    Mcal = dMcal = Xm = 0;
    for (j = 0; j < Nimlist[i]; j++) {
      m = imlist[i][j];
      NX = image[m].NX;
      NY = image[m].NY;
      dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
      XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
      Rmin = MIN (Rmin, R);
      Rmax = MAX (Rmax, R);
      Dmin = MIN (Dmin, D);
      Dmax = MAX (Dmax, D);
      XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
      Rmin = MIN (Rmin, R);
      Rmax = MAX (Rmax, R);
      Dmin = MIN (Dmin, D);
      Dmax = MAX (Dmax, D);
      XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
      Rmin = MIN (Rmin, R);
      Rmax = MAX (Rmax, R);
      Dmin = MIN (Dmin, D);
      Dmax = MAX (Dmax, D);
      XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
      Rmin = MIN (Rmin, R);
      Rmax = MAX (Rmax, R);
      Dmin = MIN (Dmin, D);
      Dmax = MAX (Dmax, D);
      Mcal  += image[m].Mcal;
      dMcal += image[m].dMcal;
      Xm    += image[m].Xm;
      /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
      image[m].Mcal  = 0.0;
      image[m].dMcal = NAN;
      image[m].Xm    = NAN_S_SHORT;
    }
    dS /= Nimlist[i];
    strcpy (mosaic[i].coords.ctype, "RA---TAN");
    mosaic[i].coords.crval1 = Rmin;
    mosaic[i].coords.crval2 = Dmin;
    mosaic[i].coords.crpix1 = 0.0;
    mosaic[i].coords.crpix2 = 0.0;
    mosaic[i].coords.cdelt1 = dS;
    mosaic[i].coords.cdelt2 = dS;
    mosaic[i].coords.pc1_1  = 1.0;
    mosaic[i].coords.pc2_2  = 1.0;
    mosaic[i].coords.pc1_2  = 0.0;
    mosaic[i].coords.pc2_1  = 0.0;
    RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);

    mosaic[i].Mcal  = Mcal / Nimlist[i];
    mosaic[i].dMcal = dMcal / Nimlist[i];
    mosaic[i].Xm    = Xm / Nimlist[i];
  }
  if (!USE_GRID) return;

  dXmax = MAX (dXmax, dX);
  dYmax = MAX (dYmax, dY);
  initGrid (dXmax, dYmax);
  return;
}
# endif

/* select all image equivalent to the current photcode */
void findImages (Catalog *catalog, int Ncatalog) {

  int i, j, ecode;

  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Nmeasure; j++) {
      ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode);
      if (photcode[0].code != ecode) continue;
      matchImage (catalog, j, i);
    }
  }
}

/* XXX convert this to use the values of measure->Xccd,Yccd */
/* XXX need to apply the conversion for ELIXIR and LONEOS formats */
void matchImage (Catalog *catalog, int meas, int cat) {

  int i, ave, ccdnum;
  double ra, dec, X, Y;
  char *pname, *filter, *p, base[256];
  Measure *measure;
  
  measure = &catalog[cat].measure[meas];
  for (i = 0; i < Nimage; i++) {
    if (image[0].photcode == -1) continue;
    if (measure[0].photcode != image[i].photcode) continue;
    if (measure[0].t < start[i]) continue;
    if (measure[0].t > stop[i]) continue;
    
# ifdef GRID_V2 /* this section is added to support GridOps.v2.c */
    if (USE_GRID) {

      /* identify the ccd on the basis of the photcode name */
      pname = GetPhotcodeNamebyCode (image[i].photcode);
      filter = photcode[0].name;
      sprintf (base, "%s.%s.", MOSAICNAME, filter);
      if (strncmp (pname, base, strlen (base))) continue;
      p = pname + strlen(base);
      if (*p == 0) continue;
      ccdnum = atoi (p);
      /* some ambiguity here between seq number and id number */

      /* add this measurement to the grid cell for this chip */
      ave = measure[0].averef;
      ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
      dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
       
      /* X,Y always positive-definite in range 0,0 - dX, dY */
      RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
      setGridMeasure (meas, cat, X, Y, ccdnum);
    }
# endif

    bin[cat][meas] = i;

    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;
  }
  /*  fprintf (stderr, "can't find source image for this measurement: %d (%d)\n", measure[0].t, measure[0].photcode); */
}

int getImageEntry (int meas, int cat) {

  int i;

  i = bin[cat][meas];
  return (i);
}

float getMcal (int meas, int cat) {

  int i;
  float value;

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

  if (image[i].code & IMAGE_BAD) return (NAN);  
  value = image[i].Mcal;
  return (value);
}

Coords *getCoords (int meas, int cat) {

  int i;

  i = bin[cat][meas];
  if (i == -1) return (NULL);
  return (&image[i].coords);
}

/* determine Mcal values for all images */
void setMcal (Catalog *catalog, int PoorImages) {

  int i, j, m, c, n, N, Nmax, mark, bad;
  float Msys, Mrel, Mmos, Mgrid;
  double *list, *dlist;
  StatType stats;

  if (FREEZE_IMAGES) return;

  if (PoorImages) {
    IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
  }

  Nmax = 0;
  for (i = 0; i < Nimage; i++) {
    Nmax = MAX (Nmax, Nlist[i]);
  }
  ALLOCATE (list, double, Nmax);
  ALLOCATE (dlist, double, Nmax);

  for (i = 0; i < Nimage; i++) {
    
    /* on PoorImages run, skip good images */
    if (PoorImages) {
      bad = image[i].code & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP);
      if (!bad) continue;
    }      

    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) continue;
      Mmos  = getMmos  (m, c);
      if (isnan(Mmos)) continue;
      Mgrid = getMgrid (m, c);
      if (isnan(Mgrid)) 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]);
      list[N] = Msys - Mrel - Mmos - Mgrid;
      dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
      N++;
    }
    /* Nlist[i] is all measurements, N is good measurements */

    /* too few good measurements or too many bad measurements */
    if (!PoorImages) {
      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
      if (mark) {
	image[i].code |= ID_IMAGE_FEW;
      } else {
	image[i].code &= ~ID_IMAGE_FEW;
      }      
    }
    
    liststats (list, dlist, N, &stats);
    image[i].Mcal  = stats.mean;
    image[i].dMcal = stats.sigma;
    image[i].Xm    = 100.0*log10(stats.chisq);
  }
  free (list);
  free (dlist);
  if (PoorImages) {
    IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
    STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
  }
  return;
}

/* mark image if: abs(Mcal) too large, dMcal too large */
void clean_images () {

  int i, N, mark, Nmark;
  double *mlist, *slist, *dlist;
  double MaxOffset, MaxScatter, MedOffset;
  StatType stats;

  if (FREEZE_IMAGES) return;

  if (VERBOSE) fprintf (stderr, "marking poor images\n");

  ALLOCATE (mlist, double, Nimage);
  ALLOCATE (slist, double, Nimage);
  ALLOCATE (dlist, double, Nimage);

  for (i = N = 0; i < Nimage; i++) {
    if (image[i].code & IMAGE_BAD) continue;
    mlist[N] = fabs (image[i].Mcal);
    slist[N] = image[i].dMcal;
    dlist[N] = 1;
    N++;
  }
  initstats ("MEAN");
  liststats (mlist, dlist, N, &stats);
  MaxOffset = MAX (IMAGE_OFFSET, 3*stats.sigma);
  MedOffset = stats.median;
  liststats (slist, dlist, N, &stats);
  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
  fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
  
  Nmark = 0;
  for (i = 0; i < Nimage; i++) {
    mark = FALSE;
    image[i].code &= ~ID_IMAGE_POOR;
    mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].Mcal - MedOffset) > MaxOffset);
    if (mark) { 
      Nmark ++;
      image[i].code |= ID_IMAGE_POOR;
    } else {
      image[i].code &= ~ID_IMAGE_POOR;
    }
  }

  fprintf (stderr, "%d images marked poor\n", Nmark);
  initstats (STATMODE);
  free (mlist);
  free (slist);
  free (dlist);
}

void plot_images () {

  int i, bin;
  double *xlist, *Mlist, *dlist;
  Graphdata graphdata;

  if (FREEZE_IMAGES) return;

  ALLOCATE (xlist, double, Nimage);
  ALLOCATE (dlist, double, Nimage);
  ALLOCATE (Mlist, double, Nimage);

  /**** dMcal vs airmass ****/
  for (i = 0; i < Nimage; i++) {
    Mlist[i] = image[i].Mcal;
    dlist[i] = image[i].dMcal;
    xlist[i] = image[i].secz;
  }

  plot_defaults (&graphdata);
  graphdata.ymin = PlotdMmin; 
  graphdata.ymax = PlotdMmax;
  plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "airmass.png");
  plot_defaults (&graphdata);
  plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", NULL);

# define NBIN 200
  REALLOCATE (xlist, double, NBIN);
  REALLOCATE (Mlist, double, NBIN);

  /**** dMcal histgram ****/
  for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i;
  bzero (Mlist, NBIN*sizeof(double));
  for (i = 0; i < Nimage; i++) {
    bin = image[i].dMcal / 0.00025;
    bin = MAX (0, MIN (NBIN - 1, bin));
    Mlist[bin] += 1.0;
  }

  plot_defaults (&graphdata);
  graphdata.style = 1;
  plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "dMcalhist.png");

  free (dlist);
  free (xlist);
  free (Mlist);
}

StatType statsImageN (Catalog *catalog) {

  int i, j, m, c, n, N;
  double *list, *dlist;
  float Mcal, Mmos, Mgrid;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (FREEZE_IMAGES) return (stats);

  ALLOCATE (list, double, Nimage);
  ALLOCATE (dlist, double, Nimage);

  n = 0;
  for (i = 0; i < Nimage; i++) {
    if (image[i].code & IMAGE_BAD)  continue;

    N = 0;
    for (j = 0; j < Nlist[i]; j++) {

      m = mlist[i][j];
      c = clist[i][j];

      Mcal  = getMcal  (m, c);
      if (isnan(Mcal)) continue;
      Mmos  = getMmos  (m, c);
      if (isnan(Mmos)) continue;
      Mgrid = getMgrid (m, c);
      if (isnan(Mgrid)) continue;
      N++;
    }
    list[n] = N;
    dlist[n] = 1;
    n++;
  }

  liststats (list, dlist, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

StatType statsImageX (Catalog *catalog) {

  int i, n;
  double *list, *dlist;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (FREEZE_IMAGES) return (stats);

  ALLOCATE (list, double, Nimage);
  ALLOCATE (dlist, double, Nimage);

  n = 0;
  for (i = 0; i < Nimage; i++) {

    if (image[i].code & IMAGE_BAD)  continue;

    list[n] = pow (10.0, 0.01*image[i].Xm);
    dlist[n] = 1;
    n++;
  }

  liststats (list, dlist, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

StatType statsImageM (Catalog *catalog) {

  int i, n;
  double *list, *dlist;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (FREEZE_IMAGES) return (stats);

  ALLOCATE (list, double, Nimage);
  ALLOCATE (dlist, double, Nimage);

  n = 0;
  for (i = 0; i < Nimage; i++) {

    if (image[i].code & IMAGE_BAD)  continue;

    list[n] = image[i].Mcal;
    dlist[n] = 1;
    n++;
  }

  liststats (list, dlist, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

StatType statsImagedM (Catalog *catalog) {

  int i, n;
  double *list, *dlist;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (FREEZE_IMAGES) return (stats);

  ALLOCATE (list, double, Nimage);
  ALLOCATE (dlist, double, Nimage);

  n = 0;
  for (i = 0; i < Nimage; i++) {

    if (image[i].code & IMAGE_BAD)  continue;

    list[n] = image[i].dMcal;
    dlist[n] = 1;
    n++;
  }

  liststats (list, dlist, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

Image *getimages (int *N) {

  *N = Nimage;
  return (image);
}

Image *getimage (int N) {
  return (&image[N]);
}

