# include "relastro.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;

Image *getimages (int *N) {
  *N = Nimage;
  return (image);
}

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

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.05*image[i].trate*image[i].NY, 1);
    stop[i]  = image[i].tzero + MAX(1.05*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);
}

/* match measurements to images */
void findImages (Catalog *catalog, int Ncatalog) {

  int i, j;

  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Nmeasure; j++) {
      matchImage (catalog, j, i);
    }
  }
}

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

  int i;
  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;
    
    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;
  }
}

Coords *getCoords (int meas, int cat) {

  int i;

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

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

/* XXX re-write to use image->nstars if that is possible */
StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) {

  int i, m, c;

  Image *mosaic;
  Coords *moscoords;
  StarData *raw;
  
  ALLOCATE (raw, StarData, Nlist[im]);

  mosaic = NULL;
  moscoords = NULL;
  if (isMosaic) {
    mosaic = getMosaicForImage (im);
    moscoords = &mosaic[0].coords;
  }

  for (i = 0; i < Nlist[im]; i++) {
    m = mlist[im][i];
    c = clist[im][i];

    /* apply the current image transformation or use the current value of R+dR, D+dD? */
    raw[i].X = catalog[c].measure[m].Xccd;
    raw[i].Y = catalog[c].measure[m].Yccd;
    
    raw[i].Mag  = catalog[c].measure[m].M;
    raw[i].dMag = catalog[c].measure[m].dM;

    /* note that for a Simple image, L,M = P,Q */
    XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
    if (isMosaic) {
      XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords);
      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
    } else {
      raw[i].P = raw[i].L;
      raw[i].Q = raw[i].M;
      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords);
    }
  }
  
  *Nstars = Nlist[im];
  return (raw);
}

StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) {

  int i, m, c, n;

  Image *mosaic;
  Coords *moscoords;
  StarData *ref;
  
  ALLOCATE (ref, StarData, Nlist[im]);

  mosaic = NULL;
  moscoords = NULL;
  if (isMosaic) {
    mosaic = getMosaicForImage (im);
    moscoords = &mosaic[0].coords;
  }

  for (i = 0; i < Nlist[im]; i++) {
    m = mlist[im][i];
    c = clist[im][i];
    n = catalog[c].measure[m].averef;

    /* apply the current image transformation or use the current value of R+dR, D+dD? */
    ref[i].R = catalog[c].average[n].R;
    ref[i].D = catalog[c].average[n].D;

    /* note that for a Simple image, L,M = P,Q */
    RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
    if (isMosaic) {
      LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, moscoords);
      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, moscoords);
    } else {
      ref[i].L = ref[i].P;
      ref[i].M = ref[i].Q;
      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
    }
  }
  
  *Nstars = Nlist[im];
  return (ref);
}
