# include "relphot.h"

Image *getimages (int *N);

static int    Nmosaic;
static Mosaic *mosaic;

static int    *Nimlist;
static int   **imlist; /* mosaic -> image[] */
static int   **bin;    /* catalog, measure -> mosaic */

static int   **clist;  /* mosaic -> catalog[] */
static int   **mlist;  /* mosiac -> measure[] */
static int    *Nlist;
static int    *NLIST;

/* find mosaic frames (unique time periods & photcode name matches mosaic) */
void initMosaics (Image *image, int Nimage) {

  int i, j, status, found, NMOSAIC, *NIMLIST;
  unsigned int start, stop;
  char *pname;

  if (!MOSAICNAME[0]) return;

  Nmosaic = 0;
  NMOSAIC = 10;
  ALLOCATE (mosaic, Mosaic, NMOSAIC);

  ALLOCATE (imlist, int *, NMOSAIC);
  ALLOCATE (Nimlist, int, NMOSAIC);
  ALLOCATE (NIMLIST, int, NMOSAIC);

  /* generate list of unique mosaics */
  for (i = 0; i < Nimage; i++) {

    /* select valid mosaic images by photcode */
    pname = GetPhotcodeNamebyCode (image[i].source);
    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
    if (status) continue;

    /* set image time range */
    start = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1);
    stop  = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1);

    /* find existing mosaic with this time range */
    found = FALSE;
    for (j = 0; !found && (j < Nmosaic); j++) { 
      if (stop  < mosaic[j].start) continue;
      if (start > mosaic[j].stop)  continue;
      found = TRUE;

      /* add image to mosaic image list */
      imlist[j][Nimlist[j]] = i;
      Nimlist[j] ++;
      if (Nimlist[j] == NIMLIST[j]) {
	NIMLIST[j] += 10;
	REALLOCATE (imlist[j], int, NIMLIST[j]);
      }

    }
    if (found) continue;
    
    /* a new mosaic, define ranges */
    mosaic[Nmosaic].start = start;
    mosaic[Nmosaic].stop  = stop;
    mosaic[Nmosaic].Mcal  = 0.0;
    mosaic[Nmosaic].dMcal = 0.0;
    mosaic[Nmosaic].Xm    = 0.0;
    mosaic[Nmosaic].code  = image[i].code;
    mosaic[Nmosaic].secz  = image[i].secz;

    /* add image to mosaic image list */
    NIMLIST[Nmosaic] = 10;
    Nimlist[Nmosaic] = 1;
    ALLOCATE (imlist[Nmosaic], int, NIMLIST[Nmosaic]);
    imlist[Nmosaic][0] = i;

    Nmosaic ++;
    if (Nmosaic == NMOSAIC) {
      NMOSAIC += 10;
      REALLOCATE (mosaic, Mosaic, NMOSAIC);
      REALLOCATE (imlist, int *, NMOSAIC);
      REALLOCATE (Nimlist, int, NMOSAIC);
      REALLOCATE (NIMLIST, int, NMOSAIC);
    }
  }

  initMosaicGrid (image, Nimage);
  return;
}

void initMosaicGrid (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 dR, dD, 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 = image[m].Xm = NO_MAG;
    }
    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;
}

void setMcalFinal (Image *image, int Nimage) {

  int i, j, im;

  if (!MOSAICNAME[0]) return;

  for (i = 0; i < Nmosaic; i++) {
    for (j = 0; j < Nimlist[i]; j++) {
      im = imlist[i][j];
      image[im].Mcal = mosaic[i].Mcal;
      image[im].dMcal = mosaic[i].dMcal;
      image[im].Xm = mosaic[i].Xm;
      image[im].code |= (mosaic[i].code & ID_IMAGE_FEW);
      image[im].code |= (mosaic[i].code & ID_IMAGE_POOR);
    }
  }      
}

void initMosaicBins (Catalog *catalog, int Ncatalog) {

  int i, j;

  /* measure -> mosaic */
  if (!MOSAICNAME[0]) return;
  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;
  }

  /* mosaic -> measure */
  ALLOCATE (Nlist, int, Nmosaic);
  ALLOCATE (NLIST, int, Nmosaic);
  ALLOCATE (clist, int *, Nmosaic);
  ALLOCATE (mlist, int *, Nmosaic);

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

void freeMosaicBins (Ncatalog) {

  int i;

  /* measure -> mosaic */
  if (!MOSAICNAME[0]) return;
  for (i = 0; i < Ncatalog; i++) {
    free (bin[i]);
  }
  free (bin);

  /* mosaic -> measure */
  for (i = 0; i < Nmosaic; i++) {
    free (clist[i]);
    free (mlist[i]);
  }
  free (Nlist);
  free (NLIST);
  free (clist);
  free (mlist);
}

int findMosaics (Catalog *catalog, int Ncatalog) {
  
  int i, j, ecode;

  if (!MOSAICNAME[0]) return (FALSE);
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Nmeasure; j++) {
      if (TimeSelect) {
	if (catalog[i].measure[j].t < TSTART) continue;
	if (catalog[i].measure[j].t > TSTOP) continue;
      }
      ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].source);
      if (photcode[0].code != ecode) continue;
      matchMosaics (catalog, j, i);
    }
  }
  
}

void matchMosaics (Catalog *catalog, int meas, int cat) {

  int i, ave;
  double ra, dec, X, Y;

  for (i = 0; i < Nmosaic; i++) {
    if (catalog[cat].measure[meas].t < mosaic[i].start) continue;
    if (catalog[cat].measure[meas].t > mosaic[i].stop) continue;
    
# ifdef GRID_V1
    if (USE_GRID) {
      ave = catalog[cat].measure[meas].averef;
      ra  = catalog[cat].average[ave].R - catalog[cat].measure[meas].dR / 360000.0;
      dec = catalog[cat].average[ave].D - catalog[cat].measure[meas].dD / 360000.0;

      /* X,Y always positive-definite in range 0,0 - dX, dY */
      RD_to_XY (&X, &Y, ra, dec, &mosaic[i].coords);
      setGridMeasure (meas, cat, X, Y);
    }
# 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, "missed measurement\n");
  return;
}

short getMmos (int meas, int cat) {

  int i;
  short value;

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

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

int setMmos (Catalog *catalog, int PoorImages) {

  int i, j, m, c, n, N, Nmax, mark, bad;
  short Msys, Mrel, Mcal, Mgrid;
  double *list, *dlist, *Mlist, *dMlist;
  StatType stats;
  Image *image;
  Graphdata graphdata;

  if (!MOSAICNAME[0]) return (FALSE);

  image = getimages (&N);

  if (PoorImages) {
    IMAGE_BAD = 0;
  }

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

  for (i = 0; i < Nmosaic; i++) {
    
    /* on PoorImages run, skip good images */
    if (PoorImages) {
      bad = mosaic[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].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;
      
      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 - Mgrid;
      dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
      Mlist[N] = 0.001*Msys;
      dMlist[N] = 0.001*list[N];
      N++;
    }
    /* Nlist[i] is all measurements, N is good measurements */

    /* too few good measurements or too many bad measurements (skip in PoorImages run) */

    if (!PoorImages) {
      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
      if (mark) {
	fprintf (stderr, "marked image %s (%d), (%d < %d) || (%d < %f*%d)\n", image[imlist[i][0]].name, i, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, Nlist[i]);
	mosaic[i].code |= ID_IMAGE_FEW;
      } else {
	mosaic[i].code &= ~ID_IMAGE_FEW;
      }
    }
    liststats (list, dlist, N, &stats);
    if (PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
    mosaic[i].Mcal  = stats.mean;
    mosaic[i].dMcal = stats.sigma;
    mosaic[i].Xm    = 100.0*log10(stats.chisq);
  }
  free (list);
  free (dlist);
  free (Mlist);
  free (dMlist);

  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 | ID_MEAS_SKIP | ID_MEAS_AREA;
  }
  return (TRUE);
}
  
StatType statsMosaicM (Catalog *catalog) {

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

  bzero (&stats, sizeof (StatType));
  if (!MOSAICNAME[0]) return (stats);

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

  n = 0;
  for (i = 0; i < Nmosaic; i++) {
    if (mosaic[i].code & IMAGE_BAD) continue;
    list[n] = 0.001*mosaic[i].Mcal;
    dlist[n] = 1;
    n++;
  }

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

StatType statsMosaicdM (Catalog *catalog) {

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

  bzero (&stats, sizeof (StatType));
  if (!MOSAICNAME[0]) return (stats);

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

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

    if (mosaic[i].code & IMAGE_BAD) continue;
    list[n] = 0.001*mosaic[i].dMcal;
    dlist[n] = 1;
    n++;
  }

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

StatType statsMosaicN (Catalog *catalog) {

  int i, j, k, m, c, n, N, Nmax;
  double *list, *dlist;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (!MOSAICNAME[0]) return (stats);

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

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

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

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

      if (getMcal  (m, c) == NO_MAG) continue;
      if (getMgrid (m, c) == NO_MAG) continue;
      if (getMrel  (catalog, m, c) == NO_MAG) continue;
      N++;
    }
    list[n] = N;
    dlist[n] = 1;
    n++;
  }
  fprintf (stderr, "Nmosaic: %d, n: %d\n", Nmosaic, n);

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

StatType statsMosaicX (Catalog *catalog) {

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

  bzero (&stats, sizeof (StatType));
  if (!MOSAICNAME[0]) return (stats);

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

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

    if (mosaic[i].code & IMAGE_BAD) continue;
    list[n] = pow(10.0, 0.01*mosaic[i].Xm);
    dlist[n] = 1;
    n++;
  }

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

/* mark mosaic if: abs(Mcal - <Mcal>) too large, dMcal too large */
void clean_mosaics () {

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

  if (!MOSAICNAME[0]) return;

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

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

  for (i = N = 0; i < Nmosaic; i++) {
    if (mosaic[i].code & IMAGE_BAD) continue;
    mlist[N] = mosaic[i].Mcal;
    slist[N] = mosaic[i].dMcal;
    dlist[N] = 1;
    N++;
  }
  initstats ("MEAN");
  liststats (mlist, dlist, N, &stats);
  MaxOffset = MAX (IMAGE_OFFSET, 2*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", 0.001*MedOffset, 0.001*stats.median, 0.001*MaxScatter, 0.001*MaxOffset);
  
  Nmark = 0;
  for (i = 0; i < Nmosaic; i++) {
    mark = FALSE;
    mark = (mosaic[i].dMcal > MaxScatter) || (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset);
    if (mark) { 
      Nmark ++;
      mosaic[i].code |= ID_IMAGE_POOR;
    } else {
      mosaic[i].code &= ~ID_IMAGE_POOR;
    }
  }

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

void plot_mosaic_fields (Catalog *catalog) {

  int i, j, m, c, N, ave, Nimage;
  short Mrel, Mcal, Mmos;
  double *xlist, *ylist;
  double Xmin, Xmax, Ymin, Ymax;
  char string[64];
  Image *image;
  Graphdata graphdata;

  if (!MOSAICNAME[0]) return;

  image = getimages (&Nimage);

  N = 0;
  for (i = 0; i < Nmosaic; i++) 
    N = MAX (N, Nlist[i]);

  ALLOCATE (xlist, double, N);
  ALLOCATE (ylist, double, N);

  for (i = 0; i < Nmosaic; i++) {
    N = 0;
    Xmin = Ymin = +360.0;
    Xmax = Ymax = -360.0;
    for (j = 0; j < Nlist[i]; j++) {
      
      m = mlist[i][j];
      c = clist[i][j];
      
      ave = catalog[c].measure[m].averef;
      xlist[N] = catalog[c].average[ave].R - catalog[c].measure[m].dR / 360000.0;
      ylist[N] = catalog[c].average[ave].D - catalog[c].measure[m].dD / 360000.0;
      N++;
    }
  
    sprintf (string, "Mosaic %d", i);
    plot_defaults (&graphdata);
    plot_list (&graphdata, xlist, ylist, N, string, NULL);
  }

  free (ylist);
  free (xlist);
}

void plot_mosaics () {

  int i, j, m, c, N, Nmax, bin;
  short Mrel, Mcal, Mmos;
  double *xlist, *Mlist, *dlist;
  Graphdata graphdata;

  if (!MOSAICNAME[0]) return;

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

  for (i = 0; i < Nmosaic; i++) {
    Mlist[i] = 0.001*mosaic[i].Mcal;
    dlist[i] = 0.001*mosaic[i].dMcal;
    xlist[i] = 0.001*mosaic[i].secz;
  }

  plot_defaults (&graphdata);
  graphdata.xmin = 0.95;
  graphdata.xmax = 2.50;
  graphdata.ymin = PlotdMmin;
  graphdata.ymax = PlotdMmax;
  plot_list (&graphdata, xlist, Mlist, Nmosaic, "airmass vs Mcal", "airmass.png");
  plot_defaults (&graphdata);
  plot_list (&graphdata, Mlist, dlist, Nmosaic, "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.0025*i;
  bzero (Mlist, NBIN*sizeof(double));
  for (i = 0; i < Nmosaic; i++) {
    bin = 0.001*mosaic[i].dMcal / 0.0025;
    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);
}

int *SelectRefMosaic (Mosaic **refmosaic, int *Nimage) {

  int i, Imax, Nmax;

  Imax = 0;
  Nmax = Nimlist[0];
  for (i = 0; i < Nmosaic; i++) {
    if (Nimlist[i] > Nmax) {
      Imax = i;
      Nmax = Nimlist[i];
    }
  }

  *refmosaic = &mosaic[Imax];
  *Nimage = Nmax;
  return (imlist[Imax]);
}
