# include "dvoshell.h"

int mextract (int argc, char **argv) {
  
  int i, j, k, m, n, N, Npts, NPTS, last, Nfields, Nreturn, Ncstack, Nstack;
  int Nsecfilt, VERBOSE, loadImages, mosaicMode;
  char **cstack, name[1024];
  float *values;
  void *Signal;

  Catalog catalog;
  SkyList *skylist;
  PhotCode *code;
  Vector **vec;
  dbField *fields;
  dbStack *stack;
  SkyRegionSelection *selection;

  /* defaults */
  skylist = NULL;
  code = NULL;
  fields = NULL;
  stack = NULL;

  if ((N = get_argument (argc, argv, "-h"))) goto help;
  if ((N = get_argument (argc, argv, "--help"))) goto help;

  VERBOSE = FALSE;
  if ((N = get_argument (argc, argv, "-v"))) {
    remove_argument (N, &argc, argv);
    VERBOSE = TRUE;
  }

  dvo_catalog_init (&catalog, TRUE);

  /* load photcode information */
  if (!InitPhotcodes ()) goto escape;
  Nsecfilt = GetPhotcodeNsecfilt ();
  
  // init locally static variables (time refs)
  dbExtractMeasuresInit();

  // parse skyregion options
  if ((selection = SetRegionSelection (&argc, argv)) == NULL) goto escape;

  // command-line is of the form: avextract field,field, field [where (field op value)...]

  // parse the fields to be extracted and returned
  fields = dbCmdlineFields (argc, argv, DVO_TABLE_MEASURE, &last, &Nfields);
  if (fields == NULL) return (FALSE);

  // XXX add the skyregion limits as if it were a where statement

  // parse the remainder of the line as a boolean math expression
  cstack = isolate_elements (argc-last, &argv[last], &Ncstack);
  
  // construct the db Boolean math stack (frees cstack)
  stack = dbRPN (Ncstack, cstack, &Nstack);
  if ((Ncstack > 0) && (Nstack < 1)) {
    print_error ();
    goto escape;
  }

  // add the skyregion limits to the where statement (or create)
  dbAstroRegionLimits (&stack, &Nstack, selection, DVO_TABLE_MEASURE);

  // parse stack elements into fields and scalars as needed
  Nreturn = Nfields; 
  dbCheckStack (stack, Nstack, DVO_TABLE_MEASURE, &fields, &Nfields);
  // XXX handle errors

  /* load region corresponding to selection above */
  if ((skylist = SelectRegions (selection)) == NULL) goto escape;

  // load image data if needed (for fields listed below)
  loadImages = FALSE;
  mosaicMode = FALSE;
  for (i = 0; !loadImages && (i < Nfields); i++) {
    if (fields[i].ID == MEAS_XCCD) loadImages = TRUE;
    if (fields[i].ID == MEAS_YCCD) loadImages = TRUE;
    if (fields[i].ID == MEAS_XMOSAIC) loadImages = mosaicMode = TRUE;
    if (fields[i].ID == MEAS_YMOSAIC) loadImages = mosaicMode = TRUE;
  }
  if (loadImages && !SetImageSelection (mosaicMode, selection)) goto escape;

  /* create storage vector */
  ALLOCATE (values, float, Nfields);
  ALLOCATE (vec, Vector *, Nreturn);
  for (i = 0; i < Nreturn; i++) {
    if (ISNUM(fields[i].name[0])) {
      sprintf (name, "v_%s", fields[i].name);
    } else {
      sprintf (name, "%s", fields[i].name);
    }
    if ((vec[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) goto escape;
  }

  Npts = 0;
  NPTS = 1;

  // grab data from all selected sky regions
  Signal = signal (SIGINT, handle_interrupt);
  interrupt = FALSE;
  for (i = 0; (i < skylist[0].Nregions) && !interrupt; i++) {
    /* lock, load, unlock catalog */
    catalog.filename = skylist[0].filename[i];
    catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF;
    catalog.Nsecfilt = Nsecfilt;

    if (VERBOSE) gprint (GP_ERR, "trying %s (%d of %d)\n", catalog.filename, i, skylist[0].Nregions);
      
    // an error exit status here is a significant error
    if (!dvo_catalog_open (&catalog, NULL, FALSE, "r")) {
      gprint (GP_ERR, "ERROR: failure to open catalog file %s\n", catalog.filename);
      exit (2);
    }
    dvo_catalog_unlock (&catalog);

    /* XXX need to call dvo_catalog_chipcoords here passing the loaded images */

    for (j = 0; (j < catalog.Naverage) && !interrupt; j++) {
      m = catalog.average[j].measureOffset;
      dbExtractMeasuresInitAve (); // reset counters for saved fields 

      for (k = 0; (k < catalog.average[j].Nmeasure); k++, m++) {

	// extract the relevant values for this measurement
	dbExtractMeasuresInitMeas (); // reset counters for saved fields 
	for (n = 0; n < Nfields; n++) {
	  values[n] = dbExtractMeasures (&catalog.average[j], &catalog.secfilt[j*Nsecfilt], &catalog.measure[m], &fields[n]);
	}
	// fprintf (stderr, "object: ave: %f, cat: %f, averef %d\n", fields[n].name, values[2], values[3], catalog.measure[m].averef);

	// test the conditional statement
	if (!dbBooleanCond (stack, Nstack, values)) continue;
	for (n = 0; n < Nreturn; n++) {
	  vec[n][0].elements[Npts] = values[n];
	  // fprintf (stderr, "keep : field: %s, value: %f\n", fields[n].name, values[n]);
	}
	Npts++;
	if (Npts >= NPTS) {
	  NPTS += 2000;
	  for (n = 0; n < Nreturn; n++) {
	    REALLOCATE (vec[n][0].elements, float, NPTS);
	  }
	}
      }
    }
    dvo_catalog_free (&catalog);
    // dbStackAllocPrint ();
    // dbStackAllocReset ();
    // dbStackFreePrint ();
    // dbStackFreeReset ();
  }
  signal (SIGINT, Signal);
  interrupt = FALSE;

  for (n = 0; n < Nreturn; n++) {
    vec[n][0].Nelements = Npts;
    REALLOCATE (vec[n][0].elements, float, MAX(1,Npts));
  }

  dbFreeFields (fields, Nfields);
  dbFreeStack (stack, Nstack);
  free (stack);
  FreeImageSelection ();
  SkyListFree (skylist);
  FreeSkyRegionSelection (selection);
  return (TRUE);

escape:
  dbFreeFields (fields, Nfields);
  dbFreeStack (stack, Nstack);
  free (stack);
  FreeImageSelection ();
  SkyListFree (skylist);
  FreeSkyRegionSelection (selection);
  dvo_catalog_free (&catalog);
  return (FALSE);

 help:
  gprint (GP_ERR, "USAGE: mextract field[,field,field...] where (expression)\n");

  if ((argc > N + 1) && !strcasecmp (argv[N+1], "fields")) {
    gprint (GP_ERR, " USAGE: avextract field[,field,field...] where (expression)\n");
    gprint (GP_ERR, "  RA : right ascension (J2000) for detection\n");
    gprint (GP_ERR, "  DEC : declination for detection\n");
    gprint (GP_ERR, "  RA:ave : average right ascension (J2000) for object\n");
    gprint (GP_ERR, "  DEC:ave : average declination for object\n");
    gprint (GP_ERR, "  RA:err : ra scatter \n");
    gprint (GP_ERR, "  DEC:err : dec scatter\n");
    gprint (GP_ERR, "  uRA : proper motion in ra\n");
    gprint (GP_ERR, "  uDEC : proper motion in dec\n");
    gprint (GP_ERR, "  duRA : proper motion error in ra\n");
    gprint (GP_ERR, "  duDEC : proper motion error in dec\n");
    gprint (GP_ERR, "  PAR : parallax\n");
    gprint (GP_ERR, "  dPAR : parallax error \n");
    gprint (GP_ERR, "  nmeas : number of measurements\n");
    gprint (GP_ERR, "  nmiss : number of non-detections\n");
    gprint (GP_ERR, "  xp : positional chi-square\n");
    gprint (GP_ERR, "  objflag : object flags\n");
    gprint (GP_ERR, "  photcode:ave : average magnitude for photcode (or equivalent)\n");
    gprint (GP_ERR, "  photcode:ref : reference magnitude system for photcode (or equivalent)\n");
    gprint (GP_ERR, "  photcode:inst : instrumental magnitude for photcode\n");
    gprint (GP_ERR, "  photcode:cat :  catalog magnitude for photcode\n");
    gprint (GP_ERR, "  photcode:sys :  system magnitude for photcode\n");
    gprint (GP_ERR, "  photcode:rel :  relative magnitude for photcode\n");
    gprint (GP_ERR, "  photcode:cal :  calibrated magnitude for photcode \n");
    gprint (GP_ERR, "  photcode:err : magnitude error for photcode\n");
    gprint (GP_ERR, "  photcode:chisq : chi-square of magnitude fit\n");
    gprint (GP_ERR, "  photcode:ncode : number of measurements in photcode\n");
    gprint (GP_ERR, "  photcode:nphot : number of measurements used for average magnitude\n");
    gprint (GP_ERR, "  airmass : airmass of detection\n");
    gprint (GP_ERR, "  exptime : exposure time\n");
    gprint (GP_ERR, "  photcode : photcode \n");
    gprint (GP_ERR, "  time : time of exposure\n");
    gprint (GP_ERR, "  dR : ra offset\n");
    gprint (GP_ERR, "  dD : dec offset\n");
    gprint (GP_ERR, "  fwhm : fwhm (average)\n");
    gprint (GP_ERR, "  fwhm_maj : fwhm (major axis)\n");
    gprint (GP_ERR, "  fwhm_min : fwhm (minor axis)\n");
    gprint (GP_ERR, "  theta : position angle\n");
    gprint (GP_ERR, "  flags : detection flags\n");
    gprint (GP_ERR, "  xccd : ccd x position\n");
    gprint (GP_ERR, "  yccd : ccd y position\n");
    gprint (GP_ERR, "  xmosaic : mosaic x position\n");
    gprint (GP_ERR, "  ymosaic : mosaic y position\n");
    gprint (GP_ERR, "  xchip : chip x position\n");
    gprint (GP_ERR, "  ychip : chip y position\n");
    gprint (GP_ERR, "  xfpa : fpa x position\n");
    gprint (GP_ERR, "  yfpa : fpa y position\n");
    return (FALSE);
  }
  gprint (GP_ERR, " mextract --help fields : for a complete listing of allowed fields\n");
  return (FALSE);
}
