Index: trunk/Ohana/src/opihi/dvo/mmextract.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/mmextract.c	(revision 19570)
+++ trunk/Ohana/src/opihi/dvo/mmextract.c	(revision 19570)
@@ -0,0 +1,310 @@
+# include "dvoshell.h"
+
+int mmextract (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, name1[1024], name2[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)...]
+  // mmextract [-noauto] field,field field where conditions matched to conditions
+  // mmextract ra,dec,mag where photcode equiv g matched to photcode == 2MASS_J
+  // mmextract -noauto ra,dec,mag where photcode equiv g matched to photcode == 2MASS_J
+  // objectID matched is implied
+  // -noauto means imageID is not matched
+
+  // parse the fields to be extracted and returned
+  fields = dbCmdlineFields (argc, argv, DVO_TABLE_MEASURE, &last, &Nfields);
+  if (fields == NULL) return (FALSE);
+
+  // XXX require a 'where' and a 'matched to'? if not, all are cross-matched (that's ok)
+  // XXX add the skyregion limits as if it were a where statement
+
+  // need to find the location of 'where' and 'matched to' expressions
+  for (i = last; i < argc; i++) {
+    if (!strcasecmp (argv[i], "matched")) {
+    }
+  }
+
+  whereS = N; whereE = N;
+  matchS = N; matchE = N;
+
+  // parse the 'where' and 'matched to' segments of the line as boolean math expressions
+  cstack1 = isolate_elements (Nwhere, &argv[whereS], &Ncstack1);
+  cstack2 = isolate_elements (Nmatch, &argv[matchS], &Ncstack2);
+  
+  // construct the db Boolean math stack (frees cstack)
+  stack1 = dbRPN (Ncstack1, cstack1, &Nstack1);
+  if ((Ncstack1 > 0) && (Nstack1 < 1)) {
+    print_error ();
+    goto escape;
+  }
+
+  // construct the db Boolean math stack (frees cstack)
+  stack2 = dbRPN (Ncstack2, cstack2, &Nstack2);
+  if ((Ncstack2 > 0) && (Nstack2 < 1)) {
+    print_error ();
+    goto escape;
+  }
+
+  // XXX disallow skyregion limits in the matched to expression?
+
+  // 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_base = Nfields;
+  Nreturn = 2*Nfields; // we are returning fieldi_1, fieldi_2 for the selected fields
+
+  dbCheckStack (stack1, Nstack1, DVO_TABLE_MEASURE, &fields, &Nfields);
+  dbCheckStack (stack2, Nstack2, 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_base; i++) {
+    if (ISNUM(fields[i].name[0])) {
+      sprintf (name1, "v_%s_1", fields[i].name);
+      sprintf (name2, "v_%s_2", fields[i].name);
+    } else {
+      sprintf (name1, "%s_1", fields[i].name);
+      sprintf (name2, "%s_2", fields[i].name);
+    }
+    if ((vec[2*i+0] = SelectVector (name1, ANYVECTOR, TRUE)) == NULL) goto escape;
+    if ((vec[2*i+1] = SelectVector (name2, ANYVECTOR, TRUE)) == NULL) goto escape;
+  }
+
+  Npts = 0;
+  NPTS = 1;
+
+  // we save the selected measures for each average to temporary tables 1 and 2
+  NTABLE = 100;
+  ALLOCATE (table1, double *, Nreturn_base);
+  ALLOCATE (table2, double *, Nreturn_base);
+  for (i = 0; i < Nreturn_base; i++) {
+    ALLOCATE (table1[i], double, NTABLE);
+    ALLOCATE (table2[i], double, NTABLE);
+  }
+
+  // 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 
+
+      // XXX check that we have space to keep all Nmeasure
+      if (NTABLE < catalog.average[j].Nmeasure) {
+	NTABLE = catalog.average[j].Nmeasure;
+	for (n = 0; n < Nreturn_base; n++) {
+	  REALLOCATE (table1[n], float, NTABLE);
+	  REALLOCATE (table2[n], float, NTABLE);
+	}
+      }
+
+      // extract the matching measures for this object into the temp tables
+      Nt1 = Nt2 = 0;
+      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 first conditional statement
+	if (dbBooleanCond (stack1, Nstack1, values)) {
+	  for (n = 0; n < Nreturn_base; n++) {
+	    table1[n][Nt1] = values[n];
+	    Nt1 ++;
+	    // fprintf (stderr, "keep : field: %s, value: %f\n", fields[n].name, values[n]);
+	  }
+	}
+	// test the second conditional statement
+	if (dbBooleanCond (stack2, Nstack2, values)) {
+	  for (n = 0; n < Nreturn_base; n++) {
+	    table2[n][Nt2] = values[n];
+	    Nt2 ++;
+	    // fprintf (stderr, "keep : field: %s, value: %f\n", fields[n].name, values[n]);
+	  }
+	}
+
+	// XXX now do the join :: need to filter against automatch if -noauto is selected (record the index value k to test)
+	for (n1 = 0; n1 < Nt1; n1++) {
+	  for (n2 = 0; n2 < Nt2; n2++) {
+	    for (n = 0; n < Nreturn_base; n++) {
+	      vec[n][0].elements[Npts] = table1[n][n1];
+	    }
+	    for (n = 0; n < Nreturn_base; n++) {
+	      vec[n+Nreturn_base][0].elements[Npts] = table2[n][n2];
+	    }
+	    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);
+}