# include "imregister.h"
# include "imreg.h"

static int REFCCD;

void DumpCADCTable (char *filename, RegImage *image, int *match, int Nmatch) {
  
  int i, Obsid, Nobsid, Nsubset, ref;
  char *datestr, *line, hdrname[99];
  unsigned long tsecond;
  Header header, theader;
  Matrix matrix;
  FTable table;
  RegImage *row;
  int *index, *entry, *obsid;
  MosaicLayout *layout;
  int *subset;
  double left, right, center, outer, top, bottom;
  double iqx, iqy, iqr, iqf, iqo;

  /* set REFCCD for use in GetREFIQ */
  REFCCD = MatchCCDName (SeeingREFCCD);
  if (REFCCD == -1) { 
    fprintf (stderr, "ERROR: can't get reference ccd\n");
    exit (1);
  }

  /* assign relevant mosaic layout structure */
  layout = (MosaicLayout *) NULL;
  if (!strcasecmp (Camera, "CFH12K"))  layout = CreateCFH12K ();
  if (!strcasecmp (Camera, "MegaCam")) layout = CreateMegaCam ();
  if (!strcasecmp (Camera, "MegaNorth")) layout = CreateMegaCam ();
  if (layout == (MosaicLayout *) NULL) {
    fprintf (stderr, "ERROR: invalid camera for CADC Table\n");
    exit (1);
  }

  /* create primary header */
  gfits_init_header (&header);    
  header.extend = TRUE;
  gfits_create_header (&header);
  gfits_create_matrix (&header, &matrix);
  gfits_print (&header, "NEXTEND", "%d", 1, 1);
  
  /* create table header */
  gfits_create_table_header (&theader, "TABLE", "CADC_RAW_IMAGES");
      
  /* add current date/time to header */
  str_to_time ("now", &tsecond);
  datestr = sec_to_date (tsecond);
  gfits_modify (&header,  "DATE", "%s", 1, datestr);
  gfits_modify (&theader, "DATE", "%s", 1, datestr);

  /* define table layout */
  gfits_define_table_column (&theader, "A99",   "FILENAME",       "filename in db",                  "");
  gfits_define_table_column (&theader, "A99",   "HDR_FILENAME",   "image header filename",           "");
  gfits_define_table_column (&theader, "I10",   "OBSID",          "image ID number",                 "");
  gfits_define_table_column (&theader, "F5.2",  "OBS_IQ_REFCCD",  "image quality on reference chip", "arcsec");
  gfits_define_table_column (&theader, "F5.2",  "OBS_IQ_CENTER",  "image quality center region",     "arcsec");
  gfits_define_table_column (&theader, "F5.2",  "OBS_IQ_R_RATIO", "IQ ratio (outer / center)",       "");
  gfits_define_table_column (&theader, "F5.2",  "OBS_IQ_X_RATIO", "IQ ratio (left / right)",         "");
  gfits_define_table_column (&theader, "F5.2",  "OBS_IQ_Y_RATIO", "IQ ratio (top / bottom)",         "");
  gfits_define_table_column (&theader, "F9.3",  "OBS_BG_VAL",     "background level",                "counts / pixel");

  /* define TNULL, TNVAL values */
  gfits_modify (&theader, "TNULL1",  "%s", 1, "NULL");  /* FILENAME       */
  gfits_modify (&theader, "TNULL2",  "%s", 1, "NULL");  /* HDR_FILENAME   */
  gfits_modify (&theader, "TNULL3",  "%s", 1, "0");     /* OBSID          */
  gfits_modify (&theader, "TNULL4",  "%s", 1, "0.00");  /* OBS_IQ_REFCCD  */
  gfits_modify (&theader, "TNULL5",  "%s", 1, "0.00");  /* OBS_IQ_CENTER  */
  gfits_modify (&theader, "TNULL6",  "%s", 1, "0.00");  /* OBS_IQ_R_RATIO */
  gfits_modify (&theader, "TNULL7",  "%s", 1, "0.00");  /* OBS_IQ_X_RATIO */
  gfits_modify (&theader, "TNULL8",  "%s", 1, "0.00");  /* OBS_IQ_Y_RATIO */
  gfits_modify (&theader, "TNULL9",  "%s", 1, "0.00");  /* OBS_BG_VAL     */

  gfits_modify (&theader, "TNVAL1",  "%s", 1, "NA");    /* FILENAME     */
  gfits_modify (&theader, "TNVAL2",  "%s", 1, "NA");    /* HDR_FILENAME */
  gfits_modify (&theader, "TNVAL3",  "%s", 1, "-1");    /* OBSID        */
  gfits_modify (&theader, "TNVAL4",  "%s", 1, "-1.00"); /* OBS_IQ_REFCCD  */
  gfits_modify (&theader, "TNVAL5",  "%s", 1, "-1.00"); /* OBS_IQ_CENTER  */
  gfits_modify (&theader, "TNVAL6",  "%s", 1, "-1.00"); /* OBS_IQ_R_RATIO */
  gfits_modify (&theader, "TNVAL7",  "%s", 1, "-1.00"); /* OBS_IQ_X_RATIO */
  gfits_modify (&theader, "TNVAL8",  "%s", 1, "-1.00"); /* OBS_IQ_Y_RATIO */
  gfits_modify (&theader, "TNVAL9",  "%s", 1, "-1.00"); /* OBS_BG_VAL   */

  /* create table, add data values */
  gfits_create_table (&theader, &table);
  
  /* prepare indicies to handle table data */
  GetObsIDIndex (image, match, Nmatch, &index, &entry);
  obsid = GetUniqueObsID (image, index, entry, Nmatch, &Nobsid);

  for (i = 0; i < Nobsid; i++) {
    ref = GetREFCCD (image, index, entry, Nmatch, obsid[i]);
    row = &image[ref];
    Obsid = index[obsid[i]];
    sprintf (hdrname, "%s.hdr", row[0].filename);

    subset = GetObsIDSubset (image, obsid[i], index, entry, Nmatch, &Nsubset);
    center = MosaicIQStats (image, subset, Nsubset, &layout[0].center);
    outer  = MosaicIQStats (image, subset, Nsubset, &layout[0].outer );
    top    = MosaicIQStats (image, subset, Nsubset, &layout[0].top   );
    bottom = MosaicIQStats (image, subset, Nsubset, &layout[0].bottom);
    left   = MosaicIQStats (image, subset, Nsubset, &layout[0].left  );
    right  = MosaicIQStats (image, subset, Nsubset, &layout[0].right );
    free (subset);

    iqx = (left   == 0) ? 0 : (right / left);
    iqy = (bottom == 0) ? 0 : (top / bottom);
    iqr = (center == 0) ? 0 : (outer / center);
    iqf = center * ARCSEC_PIXEL;
    iqo = row[0].fwhm*ARCSEC_PIXEL;
    line = gfits_table_print (&table, row[0].filename, hdrname, Obsid, iqo, iqf, iqr, iqx, iqy, row[0].sky);

    gfits_add_rows (&table, line, 1, strlen(line));
    free (line);
  }
  free (obsid);
  free (index);
  free (entry);

  gfits_write_header  (filename, &header);
  gfits_write_matrix  (filename, &matrix);
  gfits_write_Theader (filename, &theader);
  gfits_write_table   (filename, &table);
  exit (0);
}

int *GetObsIDSubset (RegImage *image, int start, int *index, int *entry, int Nindex, int *Nsubset) {

  int i, N, NSUBSET;
  int *subset;

  /* create output index */
  N = 0;
  NSUBSET = 64;
  ALLOCATE (subset, int, NSUBSET);

  /* find unique sequences */
  for (i = start; (i < Nindex) && (index[i] == index[start]); i++) {
    subset[N] = entry[i];
    N++;
    if (N == NSUBSET) {
      NSUBSET += 64;
      REALLOCATE (subset, int, NSUBSET);
    }
  }

  *Nsubset = N;
  return (subset);
}

/* return an index list of unique obs id entries at start of sequence */
int *GetUniqueObsID (RegImage *image, int *index, int *entry, int Nindex, int *Nmatch) {
  
  int i, j, N, NMATCH;
  int *match;

  /* create output index */
  N = 0;
  NMATCH = 1000;
  ALLOCATE (match, int, NMATCH);

  /* find unique sequences */
  for (i = 0; i < Nindex; ) {
    for (j = i + 1; (j < Nindex) && (index[i] == index[j]); j++);

    /* add unique entry to output list */
    match[N] = i;
    N ++;
    if (N == NMATCH) {
      NMATCH += 1000;
      REALLOCATE (match, int, NMATCH);
    }

    /* j always points to the next entry */
    i = j;
  }
  *Nmatch = N;
  return (match);
}

/* return seeing for ccd == REFCCD */
int GetREFCCD (RegImage *image, int *index, int *entry, int Nindex, int start) {
  
  int i, N, NMATCH;
  int *match;

  /* create output index */
  N = 0;
  NMATCH = 1000;
  ALLOCATE (match, int, NMATCH);

  /* find unique sequences */
  for (i = start; (i < Nindex) && (index[i] == index[start]); i++) {
    if (image[entry[i]].ccd == REFCCD) return (entry[i]);
  }
  return (start);
}

void GetObsIDIndex (RegImage *image, int *match, int Nmatch, int **Index, int **Entry) {

  int i;
  int *index, *entry;

  /* index = OBSID */
  ALLOCATE (index, int, Nmatch);
  ALLOCATE (entry, int, Nmatch);
  for (i = 0; i < Nmatch; i++) {
    index[i] = atoi (image[match[i]].filename);
    if (index[i] < 400000) fprintf (stderr, "warning: derived obsid < 400000\n");
    entry[i] = match[i];
  }
  sortpair (index, entry, Nmatch);
  *Index = index;
  *Entry = entry;
}

/* match is a list of image entries with the same obsid */
double MosaicIQStats (RegImage *image, int *match, int Nmatch, MosaicRegion *region) {

  int i, j, N, Nccd;
  double *list, value;

  Nccd = region[0].Nccd;
  ALLOCATE (list, double, Nccd);

  N = 0;
  for (i = 0; i < Nccd; i++) {
    for (j = 0; j < Nmatch; j++) {
      if (image[match[j]].ccd == region[0].ccd[i]) {
	list[N] = image[match[j]].fwhm;
	N++;
	break;
      }
    }
  }

  value = SigmaClipList (list, N);
  free (list);
  return (value);
}

double SigmaClipList (double *list, int N) {

  int i, n;
  double median, sigma3, m1, m2;

  if (N == 0) return (0.0);
  if (N == 1) return (list[0]);
  
  dsort (list, N);
  median = list[(int)(0.5*N)];
  
  m1 = m2 = 0;
  for (i = 0; i < N; i++) { m1 += list[i]; m2 += SQ(list[i]); }
  sigma3 = 3*sqrt (m2/N - m1*m1/N/N);
  
  m1 = n = 0;
  for (i = 0; i < N; i++) { 
    if (abs(list[i] - median) > sigma3) continue;
    m1 += list[i];
    n ++;
  }

  m2 = m1 / n;
  return (m2);
}
  

/* the CADC table is special: we need to report specfic IQ stats 
   which represent the variations in focus across the mosaic. 
   this representation is explicitly dependent on the mosaic,
   and does not make sense for other camera types 
   CADC table option cannot be combined with CCD filtering options
   CADC table forces -unique, needed to make calculation
*/

/* derived CADC parameters:
   OBS_IQ_CENTER - sigma-clip mean value of center region
   OBS_IQ_X      - ratio of left to right regions
   OBS_IQ_Y      - ratio of top to bottom regions
   OBS_IQ_R      - ratio of center to edge regions
*/

/* cfh12k:

   00 01 02 03 04 05
   06 07 08 09 10 11

   center: 02,03,08,09
   outer:  00,01,04,05,06,07,10,11
   top;    00-05
   bottom: 06-11
   right:  00,01,02,06,07,08
   left:   03,04,05,09,10,11

   for each region, calculate median, sigma, reject >3sigma, calculate mean
*/

/* we are guaranteed a unique set of filename / ccd values */
/* index = OBSID = atoi (filename) */

