# include "relastro.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) */
/* XXX what is a mosaic?  do we need to use the 'DIS/WRP' info to track? */
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].photcode);
    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);
    }
  }
  return;
}

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 (int 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;

  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;
      }
      matchMosaics (catalog, j, i);
    }
  }
  return (TRUE);
}

/* modify this function to use the measure->imageID field ? */
void matchMosaics (Catalog *catalog, int meas, int cat) {

  int i;

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

void plot_mosaic_fields (Catalog *catalog) {

  int i, j, m, c, N, ave, Nimage;
  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];
      
      if (catalog[c].measure[m].flags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;

      ave = catalog[c].measure[m].averef;
      xlist[N] = catalog[c].average[ave].R - catalog[c].measure[m].dR / 3600.0;
      ylist[N] = catalog[c].average[ave].D - catalog[c].measure[m].dD / 3600.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, bin;
  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] = mosaic[i].Mcal;
    dlist[i] = mosaic[i].dMcal;
    xlist[i] = 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);
  graphdata.size = 1.5;
  graphdata.ptype = 7;
  plot_list (&graphdata, Mlist, dlist, Nmosaic, "Mcal vs dMcal", "MdM.png");

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

  /**** dMcal histgram ****/
  for (i = 0; i < NBIN; i++) xlist[i] = 0.00005*i;
  bzero (Mlist, NBIN*sizeof(double));
  for (i = 0; i < Nmosaic; i++) {
    bin = mosaic[i].dMcal / 0.00005;
    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);
}

Image *getMosaicForImage (int Nim) {

    Image *mosaic;

    mosaic = NULL;
    return mosaic;
}

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]);
}
