Index: /trunk/Ohana/src/opihi/dvo/Makefile
===================================================================
--- /trunk/Ohana/src/opihi/dvo/Makefile	(revision 40407)
+++ /trunk/Ohana/src/opihi/dvo/Makefile	(revision 40408)
@@ -27,4 +27,5 @@
 $(SRC)/mySequence.$(ARCH).o		\
 $(SRC)/photcode_ops.$(ARCH).o	  	\
+$(SRC)/PeriodogramOps.$(ARCH).o	  	\
 $(SRC)/find_matches.$(ARCH).o           
 
@@ -46,4 +47,5 @@
 $(SRC)/avmatch.$(ARCH).o	  	\
 $(SRC)/avperiodogram.$(ARCH).o	  	\
+$(SRC)/avperiodomatch.$(ARCH).o	  	\
 $(SRC)/badimages.$(ARCH).o	  	\
 $(SRC)/catdir.$(ARCH).o             	\
Index: /trunk/Ohana/src/opihi/dvo/avperiodogram.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/avperiodogram.c	(revision 40407)
+++ /trunk/Ohana/src/opihi/dvo/avperiodogram.c	(revision 40408)
@@ -1,28 +1,3 @@
 # include "dvoshell.h"
-# define MIN_VAR 1e-8
-# define NPERIODS 50
-
-/* basic outline of the code:
-
- * parse optional arguments
- * -parallel : launch parallel jobs [below is called for non-parallel or for dvo_client
-
- * parse optional -where-measure clause
- * 
-
- */
-
-typedef struct {
-  int Nperiods;
-  opihi_flt *period;
-  opihi_flt *power;
-} PeriodogramResult;
-
-void PeriodogramResultFree (PeriodogramResult *result);
-void PeriodogramResultSave (PeriodogramResult *result, char *extname, FILE *foutput, Average *average);
-PeriodogramResult *periodogram_fm_raw (opihi_flt *time, opihi_flt *flux, opihi_flt *dflux, int Npts);
-
-float minP = 0.2;
-float maxP = 200.0;
 
 int avperiodogram (int argc, char **argv) {
@@ -57,4 +32,8 @@
     PARALLEL = TRUE;
   }
+
+  /* load photcode information (this needs to come before the SelectRegions command below to define the catdir */
+  if (!InitPhotcodes ()) goto escape;
+  int Nsecfilt = GetPhotcodeNsecfilt ();
 
   // parse skyregion options.  NOTE: this is stripped off in parallel operation and always
@@ -70,7 +49,7 @@
 
   // **** launch parallel / remote jobs ****
-  // this does all the work of re-packaging the command, calling it on the remote machines, then loading in the results
+  // the result files are left behind for the user to access (no automatic re-merge)
   if (PARALLEL && !HOST_ID) {
-    int status = HostTableParallelOps (skylist, argc, argv, NULL, TRUE, 0, VERBOSE);
+    int status = HostTableParallelOps (skylist, argc, argv, NULL, FALSE, 0, VERBOSE);
     SkyListFree (skylist);
     FreeSkyRegionSelection (selection);
@@ -82,5 +61,5 @@
   if ((N = get_argument (argc, argv, "-min-period"))) {
     remove_argument (N, &argc, argv);
-    minP = atof (argv[N]);
+    PeriodogramSetOptions (atof (argv[N]), NAN);
     remove_argument (N, &argc, argv);
   }
@@ -88,5 +67,5 @@
   if ((N = get_argument (argc, argv, "-max-period"))) {
     remove_argument (N, &argc, argv);
-    maxP = atof (argv[N]);
+    PeriodogramSetOptions (NAN, atof (argv[N]));
     remove_argument (N, &argc, argv);
   }
@@ -179,10 +158,7 @@
   dvo_catalog_init (&catalog, TRUE);
 
-  /* load photcode information */
-  if (!InitPhotcodes ()) goto escape;
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
   // init locally static variables (time refs)
   dbExtractAveragesInit (); 
+  dbExtractMeasuresInit(HOST_ID);
 
   // examine line for 'where' or 'match to'.  'match to' is forbidden
@@ -207,4 +183,6 @@
 
   // parse stack elements into fields and scalars as needed
+  Nfields = 0;
+  ALLOCATE (fields, dbField, 1);
   if (!dbCheckStack (stack, Nstack, DVO_TABLE_AVERAGE, &fields, &Nfields)) goto escape;
 
@@ -288,7 +266,12 @@
 	// The first 3 values[] / fields[] will be time, mag, dmag
 	for (off_t n = 0; n < Nmfields; n++) {
-	  mvalues[n] = dbExtractMeasures (average, secfilt, measure, NULL, NULL, &mfields[n]);
+	  mvalues[n] = dbExtractMeasures (average, secfilt, &measure[k], NULL, NULL, &mfields[n]);
 	}
 	if (!dbBooleanCond (mstack, Nmstack, mvalues)) continue;
+
+	// do not allow any NAN or Inf values to polute the result
+	if (!isfinite(mvalues[0].Flt)) continue;
+	if (!isfinite(mvalues[1].Flt)) continue;
+	if (!isfinite(mvalues[2].Flt)) continue;
 
 	time[Npts] = mvalues[0].Flt;
@@ -299,5 +282,12 @@
       }
 
-      PeriodogramResult *result = periodogram_fm_raw (time, mag, dmag, Npts);
+      PeriodogramResult *result = PeriodogramRawFloatingMean (time, mag, dmag, Npts);
+      free (time);
+      free (mag);
+      free (dmag);
+
+      if (!result) continue;
+
+      if (VERBOSE2) gprint (GP_ERR, "periodogram for %d, %d\n", (int) average->catID, (int) average->objID);
 
       // XXX if we have 1e6 objects, we will have output files with sizes of ~80GB
@@ -308,18 +298,26 @@
       Nobject ++;
     }
-
-    fclose (foutput);
-    fflush (foutput);
-
     dvo_catalog_free (&catalog);
   }
   ClearInterrupt (old_sigaction);
 
-  if (values) free (values);
+  fclose (foutput);
+  fflush (foutput);
+
+  // free measure stack stuff
+  FREE (mvalues);
+  dbFreeFields(mfields, Nmfields);
+  dbFreeStack(mstack, Nmstack);
+  FREE (mstack);
+
+  FREE (values);
   dbFreeFields (fields, Nfields);
   dbFreeStack (stack, Nstack);
-  if (stack) free (stack);
+  FREE (stack);
+
   SkyListFree (skylist);
   FreeSkyRegionSelection (selection);
+
+  free (output);
   return (TRUE);
 
@@ -341,166 +339,2 @@
   return (FALSE);
 }
-
-// low-level periodogram-fm function
-PeriodogramResult *periodogram_fm_raw (opihi_flt *time, opihi_flt *flux, opihi_flt *dflux, int Npts) {
-
-  /* find the max baseline, sum the inverse variances */
-  double minT = time[0];
-  double maxT = time[0];
-  opihi_flt Weight = 0;
-
-  opihi_flt *tv =  time;
-  opihi_flt *df = dflux;
-  opihi_flt *fv =  flux;
-
-  for (int i = 0; i < Npts; i++, tv++, df++) {
-    minT = MIN (minT, *tv);
-    maxT = MAX (maxT, *tv);
-    // skip points with 0.0 error? or add minimum variance?
-    Weight += 1.0 / (SQ(*df) + MIN_VAR); // MIN_VAR : added in quadrature to avoid Inf (XXX make this a user parameter?)
-  }
-
-  double WeightInv = 1.0 / Weight;
-  double dTime = maxT - minT;
-  if (dTime == 0) {
-    gprint (GP_ERR, "ERROR: time range is zero\n");
-    return NULL;
-  }
-
-  int Np = 0;
-  int NP = NPERIODS;
-  ALLOCATE_PTR (period, double, NP);
-  ALLOCATE_PTR (power, double, NP);
-
-  // for testing, we are going to write out the terms explicitly 
-  // for optimization, some of the trig functions can be calculated more efficiently 
-  // by storing cos,sin(w t) and using some recurrence rules
-
-  // maxP = minP * dP^n [n = 0 -- Nperiods-1]
-  // log(maxP) = log(minP) + (Nperiods - 1) * log (dP)
-  // log(dP) = (log(maxP) - log(minP)) / (Nperiods - 1)
-  // double dP = LINEAR ? (maxP - minP)/(NPERIODS - 1) : pow(10.0, ((log10(maxP) - log10(minP))/(NPERIODS - 1)));
-
-  double P = minP;
-  for (Np = 0; (Np < NPERIODS) && (P <= maxP); Np++) {
-    double w = 2*M_PI/P;
-    
-    /* find the period offset tau  */
-    tv = time;
-    df = dflux;
-
-    double cs = 0.0, sn = 0.0, sn2 = 0.0, cs2 = 0.0;
-
-    for (int i = 0; i < Npts; i++, tv++, df++) {
-      opihi_flt wt = WeightInv / (SQ(*df) + MIN_VAR);
-      cs  += wt*cos (*tv*w);
-      sn  += wt*sin (*tv*w);
-      cs2 += wt*cos (*tv*w*2);
-      sn2 += wt*sin (*tv*w*2);
-    }
-    double ytan = sn2 - 2*cs*sn;
-    double xtan = cs2 - (SQ(cs) - SQ(sn));
-    double tau = 0.5*atan2 (ytan, xtan) / w;
-      
-    /* find the power at this period */
-    tv = time;
-    df = dflux;
-    fv =  flux;
-
-    // YY = YY_s - Y_s*Y_s
-    // CC = CC_s - C_s*C_s
-    // SS = SS_s - S_s*S_s
-    // YC = YC_s - Y_s*C_s
-    // YS = YS_s - Y_s*S_s
-    
-    double YY_s = 0.0, CC_s = 0.0, SS_s = 0.0, YC_s = 0.0, YS_s = 0.0;
-    double Y_s = 0.0, C_s = 0.0, S_s = 0.0;
-    for (int i = 0; i < Npts; i++, tv++, fv++, df++) {
-      opihi_flt wt = WeightInv / (SQ(*df) + MIN_VAR);
-      cs = cos (w*(*tv-tau));
-      sn = sin (w*(*tv-tau));
-
-      double wtf = *fv*wt;
-      Y_s += wtf;
-      C_s += wt*cs;
-      S_s += wt*sn;
-
-      YY_s += wtf*(*fv);
-      YC_s += wtf*cs;
-      YS_s += wtf*sn;
-
-      CC_s += wt*cs*cs;
-      SS_s += wt*sn*sn;
-    }
-    double YY = YY_s - SQ(Y_s);
-    double YC = YC_s - Y_s*C_s;
-    double YS = YS_s - Y_s*S_s;
-    double CC = CC_s - SQ(C_s);
-    double SS = SS_s - SQ(S_s);
-
-    double Pa = SQ(YC) / CC;
-    double Pb = SQ(YS) / SS;
-    double Po = (Pa + Pb) / YY;
-
-    power[Np] = Po;
-    period[Np] = P;
-    if (Np >= NP) {
-      NP += 100;
-      REALLOCATE (power, opihi_flt, NP);
-      REALLOCATE (period, opihi_flt, NP);
-    }
-
-    float dP = 0.4 * SQ(P) / dTime;
-    P += dP;
-  }
-
-  ALLOCATE_PTR (result, PeriodogramResult, 1);
-  result->Nperiods = Np;
-  result->period = period;
-  result->power = power;
-  
-  return result;
-}
-
-void PeriodogramResultFree (PeriodogramResult *result) {
-  if (!result) return;
-
-  FREE (result->period);
-  FREE (result->power);
-  free (result);
-  return;
-}
-
-// write the period and power vectors to an extension in this output file
-void PeriodogramResultSave (PeriodogramResult *result, char *extname, FILE *foutput, Average *average) {
-
-  // generate a binary table 
-  Header outheader;
-  FTable outtable;
-  gfits_init_header (&outheader);
-  gfits_init_table  (&outtable);
-  outtable.header = &outheader;
-
-  gfits_create_table_header (&outheader, "BINTABLE", extname);
-
-  // add some metadata derived from average:
-  gfits_modify (&outheader, "RA_OBJ",  "%lf", 1, average->R);
-  gfits_modify (&outheader, "DEC_OBJ", "%lf", 1, average->D);
-  gfits_modify (&outheader, "OBJ_ID",  "%d", 1, average->objID);
-  gfits_modify (&outheader, "CAT_ID",  "%d", 1, average->catID);
-
-  gfits_define_bintable_column (&outheader, "D", "PERIOD", NULL, NULL, 1.0, 0.0);
-  gfits_define_bintable_column (&outheader, "D", "POWER", NULL, NULL, 1.0, 0.0);
-
-  gfits_create_table (&outheader, &outtable);
-
-  gfits_set_bintable_column (&outheader, &outtable, "PERIOD", result->period, result->Nperiods);
-  gfits_set_bintable_column (&outheader, &outtable, "POWER",  result->power,  result->Nperiods);
-
-  // write the actual table data
-  gfits_fwrite_Theader (foutput, &outheader);
-  gfits_fwrite_table   (foutput, &outtable);
-    
-  gfits_free_header (&outheader);
-  gfits_free_table  (&outtable);
-}
Index: /trunk/Ohana/src/opihi/dvo/avperiodomatch.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/avperiodomatch.c	(revision 40408)
+++ /trunk/Ohana/src/opihi/dvo/avperiodomatch.c	(revision 40408)
@@ -0,0 +1,370 @@
+# include "dvoshell.h"
+
+int avperiodomatch (int argc, char **argv) {
+  
+  int N;
+
+  SkyList *skylist = NULL;
+  SkyRegionSelection *selection = NULL;
+
+  int Ninvec = 0;
+
+  Vector **invec = NULL;
+  Vector *RAvec = NULL;
+  Vector *DECvec = NULL;
+
+  // **** parse the optional arguments ****
+  if ((N = get_argument (argc, argv, "-h"))) goto help;
+  if ((N = get_argument (argc, argv, "--help"))) goto help;
+
+  int VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = TRUE;
+  }
+  int VERBOSE2 = FALSE;
+  if ((N = get_argument (argc, argv, "-vv"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = TRUE;
+    VERBOSE2 = TRUE;
+  }
+
+  int PARALLEL = FALSE;
+  if ((N = get_argument (argc, argv, "-parallel"))) {
+    remove_argument (N, &argc, argv);
+    PARALLEL = TRUE;
+  }
+
+  /* load photcode information (this needs to come before the SelectRegions command below to define the catdir */
+  if (!InitPhotcodes ()) goto escape;
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  // parse skyregion options.  NOTE: this is stripped off in parallel operation and always
+  // defined for the client via the -skyregion option.  The dvo_client parses this
+  // argument in the main program, before it is passed to the command (like mextract)
+  if ((selection = SetRegionSelection (&argc, argv)) == NULL) {
+    print_error(); 
+    goto escape; 
+  }
+
+  /* load region corresponding to selection above */
+  if ((skylist = SelectRegions (selection)) == NULL) goto escape;
+
+  // dump results directly to fits file (esp for parallel dvo)
+  char *CoordsFile = NULL;
+  if ((N = get_argument (argc, argv, "-coords"))) {
+    remove_argument (N, &argc, argv);
+    CoordsFile = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  // XXX arg-parsing is a bit confusing with the -where-measure option
+  if (!CoordsFile && (argc < 5)) goto help;
+  if ( CoordsFile && (argc < 3)) goto help;
+
+  RAvec  = NULL;
+  DECvec = NULL;
+  // get vectors corresponding to coordinates of interest
+  if (CoordsFile) {
+    // read RAvec, DECvec from coords file (1st 2 fields?)
+    Ninvec = 0;
+    invec = ReadVectorTableFITS (CoordsFile, "COORDS", &Ninvec);
+    RAvec = invec[0];
+    DECvec = invec[1];
+  } else {
+    if ((RAvec  = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto help;
+    if ((DECvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) goto help;
+    remove_argument (2, &argc, argv); // strip off RA
+    remove_argument (2, &argc, argv); // strip off DEC
+  }
+
+  // **** launch parallel / remote jobs ****
+  // the result files are left behind for the user to access (no automatic re-merge)
+  if (PARALLEL && !HOST_ID) {
+    // allocate the temp array and copy all but (RA) (DEC)
+    int targc = 0;
+    char **targv = NULL;
+    ALLOCATE (targv, char *, argc + 2);
+    for (int i = 0; i < argc; i++) {
+      targv[targc] = strcreate (argv[i]);
+      targc ++;
+    }
+
+    // if not specified, create the coords.fits input file
+    // NOTE: RAvec, DECvec were set above
+    if (!CoordsFile) {
+      ALLOCATE_PTR (vec, Vector *, 2);
+      vec[0] = RAvec;
+      vec[1] = DECvec;
+
+      CoordsFile = abspath("coords.fits", 1024);
+      int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, FALSE, NULL, 0);
+      if (!status) goto escape;
+      free (vec);
+    }
+
+    // add the coords file to the args list
+    targv[targc+0] = strcreate ("-coords");
+    targv[targc+1] = CoordsFile; // this gets freed with targv
+    targc += 2;
+
+    int status = HostTableParallelOps (skylist, targc, targv, NULL, FALSE, RAvec->Nelements, VERBOSE);
+
+    // free up targv
+    for (int i = 0; i < targc; i++) {
+      free (targv[i]);
+    }
+    free (targv);
+
+    SkyListFree (skylist);
+    FreeSkyRegionSelection (selection);
+    return status;
+  }
+
+  // NOTE: optional arguments below are parsed after the possible remote call so they are also passed to the dvo_client above
+
+  if ((N = get_argument (argc, argv, "-min-period"))) {
+    remove_argument (N, &argc, argv);
+    PeriodogramSetOptions (atof (argv[N]), NAN);
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-max-period"))) {
+    remove_argument (N, &argc, argv);
+    PeriodogramSetOptions (NAN, atof (argv[N]));
+    remove_argument (N, &argc, argv);
+  }
+
+  // the first three Measure fields must be time, mag, dmag (mag:err)
+  int Nmfields = 3;
+  ALLOCATE_PTR (mfields, dbField, Nmfields);
+  dbInitField (&mfields[0]); ParseMeasureField (&mfields[0], "time");
+  dbInitField (&mfields[1]); ParseMeasureField (&mfields[1], "mag");
+  dbInitField (&mfields[2]); ParseMeasureField (&mfields[2], "mag:err");
+  
+  int Nmstack = 0;
+  dbStack *mstack = NULL;
+  dbValue *mvalues = NULL;
+
+  // we can restrict these with a -where-measure clause
+  if ((N = get_argument (argc, argv, "-where-measure"))) {
+    remove_argument (N, &argc, argv);
+
+    // parse the remainder of the line as a boolean math expression
+    unsigned int Nmcstack;
+    char **mcstack = isolate_elements (argc-N, &argv[N], &Nmcstack);
+  
+    // construct the db Boolean math stack (frees cmstack)
+    mstack = dbRPN (Nmcstack, mcstack, &Nmstack);
+    if (Nmcstack && !Nmstack) {
+      print_error(); 
+      goto escape; 
+    }
+    
+    // parse stack elements into fields and scalars as needed (supplement mfields)
+    if (!dbCheckStack (mstack, Nmstack, DVO_TABLE_MEASURE, &mfields, &Nmfields)) goto escape;
+
+    argc = N; // hide remaining entries from the rest of the code below
+  }    
+  // output values
+  ALLOCATE (mvalues, dbValue, Nmfields);
+
+  // use remote tables, but not dvo_client..
+  int PARALLEL_LOCAL = FALSE;
+  HostTable *table = NULL;
+  if ((N = get_argument (argc, argv, "-parallel-local"))) {
+    remove_argument (N, &argc, argv);
+    PARALLEL_LOCAL = TRUE;
+
+    char *CATDIR = GetCATDIR();
+    if (!CATDIR) {
+      gprint (GP_ERR, "CATDIR is not set\n");
+      return FALSE;
+    }
+    SkyTable *sky = GetSkyTable();
+    if (!sky) {
+      gprint (GP_ERR, "failed to load sky table for database\n");
+      return FALSE;
+    }
+    table = HostTableLoad (CATDIR, sky->hosts);
+    if (!table) {
+      gprint (GP_ERR, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR);
+      return FALSE;
+    }    
+  }
+
+  // command-line is of the form: avperiodogram (output) where (field op value)...
+  if (argc != 3) goto help;
+
+  char *output = strcreate (argv[1]);
+  float RADIUS = atof (argv[2]);
+
+  // **** generate output file ****
+  // we save the results in a single FITS file, one extension per object
+  FILE *foutput = NULL;
+  if (RESULT_FILE) {
+    foutput = fopen (RESULT_FILE, "w");
+  } else {
+    foutput = fopen (output, "w");
+  }
+  // generate the PHU and write to disk
+  // XXX add some metadata here?
+  Header header;
+  Matrix matrix;
+  gfits_init_header (&header);
+  gfits_init_matrix (&matrix);
+  header.extend = TRUE;
+  gfits_create_header (&header);
+  gfits_create_matrix (&header, &matrix);
+  gfits_fwrite_header (foutput, &header);
+  gfits_fwrite_matrix (foutput, &matrix);
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+  
+  Catalog catalog;
+  dvo_catalog_init (&catalog, TRUE);
+
+  // init locally static variables (time refs)
+  dbExtractMeasuresInit(HOST_ID);
+
+  int Nobject = 0;
+  char extname[80];
+
+  int NPTS = RAvec->Nelements;
+  ALLOCATE_PTR (idxValue, off_t, NPTS);
+
+  // grab data from all selected sky regions
+  struct sigaction *old_sigaction = SetInterrupt();
+  for (off_t i = 0; (i < skylist[0].Nregions) && !interrupt; i++) {
+
+    // does this host ID match the desired location for the table?
+    if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue;
+
+    /* lock, load, unlock catalog */
+    char hostfile[1024];
+    if (PARALLEL_LOCAL) {
+      int hostID = (skylist[0].regions[i]->hostFlags & DATA_USE_BCK) ? skylist[0].regions[i]->backupID : skylist[0].regions[i]->hostID;
+      int seq = table->index[hostID];
+      HOSTDIR = table->hosts[seq].pathname;
+    }
+    snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name);
+
+    dvo_catalog_init (&catalog, TRUE);
+    catalog.filename = (HOST_ID) ? hostfile : skylist[0].filename[i];
+    catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_SECFILT | DVO_LOAD_MEASURE;
+    catalog.Nsecfilt = 0;
+
+    if (VERBOSE) gprint (GP_ERR, "trying %s ("OFF_T_FMT" of "OFF_T_FMT")\n", catalog.filename,  i,  skylist[0].Nregions);
+      
+    // an error exit status here is a significant error
+    if (!dvo_catalog_open (&catalog, NULL, VERBOSE2, "r")) {
+      gprint (GP_ERR, "ERROR: failure to open catalog file %s\n", catalog.filename);
+      exit (2);
+    }
+    dvo_catalog_unlock (&catalog);
+
+    find_matches_by_vectors (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, idxValue);
+
+    // Scan the catalog for objects which match the WHERE clause
+    for (off_t j = 0; (j < NPTS) && !interrupt; j++) {
+      off_t Ncat = idxValue[j];
+
+      if (Ncat == -1) continue; // this point is not in this catalog file
+      if (Ncat == -2) continue; // no matches to this point
+
+      dbExtractAveragesInitAve ();  // reset counters for saved fields (costs very little)
+
+      Average *average = &catalog.average[Ncat];
+
+      Measure *measure = &catalog.measure[average->measureOffset];
+
+      off_t m = Ncat*Nsecfilt;
+      SecFilt *secfilt = &catalog.secfilt[m];
+
+      // generate time (mjd), mag, dmag vectors
+      ALLOCATE_PTR (time, opihi_flt, average->Nmeasure);
+      ALLOCATE_PTR (mag,  opihi_flt, average->Nmeasure);
+      ALLOCATE_PTR (dmag, opihi_flt, average->Nmeasure);
+
+      int Npts = 0;
+
+      // this object passes the conditions above: now grab its measures which match our conditions?
+      for (off_t k = 0; (k < average->Nmeasure); k++) {
+	if (measure[k].averef != Ncat) {
+	  gprint (GP_ERR, "ERROR: inconsistent measure->average link.  Unsorted database?\n");
+	  goto escape;
+	}
+
+	// extract the relevant values for this measurement
+	dbExtractMeasuresInitMeas (); // reset counters for saved fields  (costs very little
+
+	// XXX I need to have a different boolean expression here to restrict the measurements
+	// The first 3 values[] / fields[] will be time, mag, dmag
+	for (off_t n = 0; n < Nmfields; n++) {
+	  mvalues[n] = dbExtractMeasures (average, secfilt, &measure[k], NULL, NULL, &mfields[n]);
+	}
+	if (!dbBooleanCond (mstack, Nmstack, mvalues)) continue;
+
+	// do not allow any NAN or Inf values to polute the result
+	if (!isfinite(mvalues[0].Flt)) continue;
+	if (!isfinite(mvalues[1].Flt)) continue;
+	if (!isfinite(mvalues[2].Flt)) continue;
+
+	time[Npts] = mvalues[0].Flt;
+	mag[Npts]  = mvalues[1].Flt;
+	dmag[Npts] = mvalues[2].Flt;
+
+	Npts++;
+      }
+
+      PeriodogramResult *result = PeriodogramRawFloatingMean (time, mag, dmag, Npts);
+      free (time);
+      free (mag);
+      free (dmag);
+
+      if (!result) continue;
+
+      if (VERBOSE2) gprint (GP_ERR, "periodogram for %d, %d\n", (int) average->catID, (int) average->objID);
+
+      // XXX if we have 1e6 objects, we will have output files with sizes of ~80GB
+      // XXX warn or exit if Nboject > 1e6?
+      snprintf (extname, 80, "OBJ_%06d", Nobject);
+      PeriodogramResultSave (result, extname, foutput, average);
+      PeriodogramResultFree (result);
+      Nobject ++;
+    }
+    dvo_catalog_free (&catalog);
+  }
+  ClearInterrupt (old_sigaction);
+
+  fclose (foutput);
+  fflush (foutput);
+
+  FREE (idxValue);
+
+  // free measure stack stuff
+  FREE (mvalues);
+  dbFreeFields(mfields, Nmfields);
+  dbFreeStack(mstack, Nmstack);
+  FREE (mstack);
+
+  SkyListFree (skylist);
+  FreeSkyRegionSelection (selection);
+
+  free (output);
+  return (TRUE);
+
+ escape:
+  SkyListFree (skylist);
+  FreeSkyRegionSelection (selection);
+  dvo_catalog_free (&catalog);
+
+  if (invec) FreeVectorArray (invec, Ninvec);
+  return (FALSE);
+
+
+ help:
+  gprint (GP_ERR, "USAGE: avperiodomatch (output) (RA) (DEC) (RADIUS) -where-measure (expression)\n");
+  gprint (GP_ERR, "  NOTE: the optional -where-measure (expression) MUST come after the required arguments\n");
+  return (FALSE);
+}
Index: /trunk/Ohana/src/opihi/dvo/dvo_host_utils.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/dvo_host_utils.c	(revision 40407)
+++ /trunk/Ohana/src/opihi/dvo/dvo_host_utils.c	(revision 40408)
@@ -212,6 +212,6 @@
 	continue;
       }
-      free (table->hosts[i].results);
-      table->hosts[i].results = NULL;
+      // free (table->hosts[i].results);
+      // table->hosts[i].results = NULL;
       set_int_variable (name, 1); // result file has been read
 
@@ -225,4 +225,7 @@
 	FreeVectorArray (invec, Ninvec);
       }
+    } else {
+      // free (table->hosts[i].results);
+      // table->hosts[i].results = NULL;
     }
   }
@@ -242,5 +245,5 @@
   free (vec);
 
-  free (table);
+  FreeHostTable (table);
   return TRUE;
 }
Index: /trunk/Ohana/src/opihi/dvo/init.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/init.c	(revision 40407)
+++ /trunk/Ohana/src/opihi/dvo/init.c	(revision 40408)
@@ -4,4 +4,5 @@
 int avmatch         PROTO((int, char **));
 int avperiodogram   PROTO((int, char **));
+int avperiodomatch  PROTO((int, char **));
 int badimages       PROTO((int, char **));
 int calextract      PROTO((int, char **));
@@ -68,5 +69,6 @@
   {1, "avextract",   avextract,    "extract average data values"},
   {1, "avmatch",     avmatch,      "extract average data values matched to RA,DEC points"},
-  {1, "avperiodogram", avperiodogram, "perform periodogram on objects based on restrictions"},
+  {1, "avperiodogram",  avperiodogram, "perform periodogram on objects based on restrictions"},
+  {1, "avperiodomatch", avperiodomatch, "perform periodogram on objects based on ra,dec list"},
   {1, "badimages",   badimages,    "look for images with anomalous astrometry"},
 //  {1, "calextract",  calextract,   "extract photometry calibration"},
Index: /trunk/Ohana/src/opihi/dvo/region_list.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/region_list.c	(revision 40407)
+++ /trunk/Ohana/src/opihi/dvo/region_list.c	(revision 40408)
@@ -181,4 +181,16 @@
   SkyList *skylist;
 
+  // the list of regions comes directly from a file
+  if (selection->list != NULL) {
+    skylist = SkyListLoadFile (selection->list);
+    return (skylist);
+  }
+
+  // all other options require sky to be set
+  if (!sky) {
+    gprint (GP_ERR, "CATDIR not set\n");
+    return NULL;
+  }
+
   /* determine region-file names */
   if (selection->name != NULL) {
@@ -186,9 +198,4 @@
     return (skylist);
   } 
-
-  if (selection->list != NULL) {
-    skylist = SkyListLoadFile (selection->list);
-    return (skylist);
-  }
 
   if (selection->useDisplay) {
Index: /trunk/Ohana/src/opihi/include/dvoshell.h
===================================================================
--- /trunk/Ohana/src/opihi/include/dvoshell.h	(revision 40407)
+++ /trunk/Ohana/src/opihi/include/dvoshell.h	(revision 40408)
@@ -41,4 +41,20 @@
 char *HOSTDIR;
 char *RESULT_FILE;
+
+typedef struct {
+  int Nmeasure;
+  int Nperiods;
+  opihi_flt *period;
+  opihi_flt *power;
+  double P50;
+  double R50;
+  double R90;
+} PeriodogramResult;
+
+/*** Periodogram functions (avperiodogram and avperiodomatch) ***/
+void PeriodogramResultFree (PeriodogramResult *result);
+void PeriodogramResultSave (PeriodogramResult *result, char *extname, FILE *foutput, Average *average);
+void PeriodogramSetOptions (float minimumPeriod, float maximumPeriod);
+PeriodogramResult *PeriodogramRawFloatingMean (opihi_flt *time, opihi_flt *flux, opihi_flt *dflux, int Npts);
 
 /*** dvo prototypes ***/