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

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

  int i, j;
  char *name;

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

  for (i = 0; VERBOSE && (i < Nimage); i++) {
    name = GetPhotcodeNamebyCode (image[i].photcode);
    fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i], 
	     ohana_sec_to_date(image[i].tzero), name);
  } 
}

/* 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];

  /* find the image that supplied this measurement */
  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;
    
    // index for (catalog, measure) -> image
    bin[cat][meas] = i;

    // index for image, Nentry -> catalog
    clist[i][Nlist[i]] = cat;

    // index for image, Nentry -> measure
    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;

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

// return StarData values for detections in the specified image, converting coordinates from the
// chip positions: X,Y -> L,M -> P,Q -> R,D
void fixImageRaw (Catalog *catalog, int Ncatalog, int im) {

  int i, m, c, n;
  double X, Y, L, M, P, Q, R, D, dR, dD;
  
  Mosaic *mosaic;
  Coords *moscoords, *imcoords;
  
  moscoords = NULL;
  mosaic = getMosaicForImage (im);
  if (mosaic != NULL) {
      moscoords = &mosaic[0].coords;
  }
  imcoords = &image[im].coords;

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

    X = catalog[c].measure[m].Xccd;
    Y = catalog[c].measure[m].Yccd;
    n = catalog[c].measure[m].averef;

    if (moscoords == NULL) {
      // this is a Simple image (not a mosaic)
      // note that for a Simple image, L,M = P,Q
      XY_to_LM (&L, &M, X, Y, imcoords);
      LM_to_RD (&R, &D, L, M, imcoords);
    } else {
      XY_to_LM (&L, &M, X, Y, imcoords);
      XY_to_LM (&P, &Q, L, M, moscoords);
      LM_to_RD (&R, &D, P, Q, moscoords);
    }

    // new dR, dD : test
    dR = 3600.0*(catalog[c].average[n].R - R);
    dD = 3600.0*(catalog[c].average[n].D - D);
    
    if (fabs(catalog[c].measure[m].dR - dR) > 10.0) {
      // XXXXX running into this still for last megacam exposure: wrong mosaic?
      // ???? inconsistently hitting this????
      fprintf (stderr, "!");
      // abort ();
    }
    if (fabs(catalog[c].measure[m].dD - dD) > 10.0) {
      fprintf (stderr, "*");
      // abort ();
    }

    catalog[c].measure[m].dR = dR;
    catalog[c].measure[m].dD = dD;

    if (catalog[c].measure[m].dR > +180.0*3600.0) {
      // average on high end of boundary, move star up
      R += 360.0;
      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
    }
    if (catalog[c].measure[m].dR < -180.0*3600.0) {
      // average on low end of boundary, move star down
      R -= 360.0;
      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
    }
  }  
  return;
}

// return StarData values for detections in the specified image, converting coordinates from the
// chip positions: X,Y -> L,M -> P,Q -> R,D.  This function is used by the image fitting steps, for
// which the detections have already been filtered when they were loaded (bcatalog)
StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) {

  int i, m, c, n;
  
  Mosaic *mosaic;
  Coords *moscoords;
  StarData *raw;
  
  ALLOCATE (raw, StarData, Nlist[im]);

  mosaic = NULL;
  moscoords = NULL;
  if (mode == MODE_MOSAIC) {
      mosaic = getMosaicForImage (im);
      if (mosaic == NULL) {
	fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
	exit (1);
      }
      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;
    raw[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);

    n = catalog[c].measure[m].averef;

    // an object with only one detection provides no information about the image calibration
    raw[i].mask = FALSE;
    if (catalog[c].average[n].Nmeasure < 2) {
      raw[i].mask = TRUE;
    }

    switch (mode) {
      case MODE_SIMPLE:
	/* 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);
	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);
	break;
      case MODE_MOSAIC:
	XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
	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);
	break;
    default:
      fprintf (stderr, "error: invalid mode in getImageRaw");
      abort ();
    }
  }  

  *Nstars = Nlist[im];
  return (raw);
}

// return StarData values for averages positions in the specified image, converting coordinates from
// the sky positions: R,D -> P,Q -> L,M -> X,Y

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

  int i, m, c, n;

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

  mosaic = NULL;
  moscoords = NULL;
  if (mode == MODE_MOSAIC) {
    mosaic = getMosaicForImage (im);
    if (mosaic == NULL) {
      fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
      exit (1);
    }
    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;
    
    ref[i].Mag  = catalog[c].measure[m].M;
    ref[i].dMag = catalog[c].measure[m].dM;
    ref[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);

    ref[i].mask = FALSE;

    /* note that for a Simple image, L,M = P,Q */
    switch (mode) {
      case MODE_SIMPLE:
      RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
      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);
      break;
      case MODE_MOSAIC:
      RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords);
      LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords);
      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
      break;
      default:
	fprintf (stderr, "invalid case");
	abort();
    }
  }
  
  *Nstars = Nlist[im];
  return (ref);
}
