Index: trunk/Ohana/src/opihi/cmd.data/impeaks.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/impeaks.c	(revision 40397)
+++ trunk/Ohana/src/opihi/cmd.data/impeaks.c	(revision 40398)
@@ -100,4 +100,5 @@
   int Npeaks = 0;
   for (ix = 1; ix < Nx - 1; ix++) {
+    if (!isfinite(row[ix])) continue; // ignore NAN values
     if (row[ix] <  threshold) continue;   // only accept pixels above threshold
     if (row[ix] <  row[ix - 1]) continue; // peak pixel must be at least preceeding pixel
Index: trunk/Ohana/src/opihi/cmd.data/sort.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/sort.c	(revision 40397)
+++ trunk/Ohana/src/opihi/cmd.data/sort.c	(revision 40398)
@@ -7,5 +7,7 @@
   itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \
 }
-# define COMPARE(A,B)(X[A] < X[B])
+
+// # define COMPARE(A,B)(X[A] < X[B])
+# define COMPARE(A,B)((!isfinite(X[A]) && isfinite(X[B])) || (X[A] < X[B]))
 
   OHANA_SORT (N, COMPARE, SWAPFUNC);
@@ -22,5 +24,7 @@
   itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \
 }
+
 # define COMPARE(A,B)(X[A] < X[B])
+// # define COMPARE(A,B)((!isfinite(X[A]) && isfinite(X[B])) || (X[A] < X[B]))
 
   OHANA_SORT (N, COMPARE, SWAPFUNC);
Index: trunk/Ohana/src/opihi/dvo/Makefile
===================================================================
--- trunk/Ohana/src/opihi/dvo/Makefile	(revision 40397)
+++ trunk/Ohana/src/opihi/dvo/Makefile	(revision 40398)
@@ -45,4 +45,5 @@
 $(SRC)/avextract.$(ARCH).o	  	\
 $(SRC)/avmatch.$(ARCH).o	  	\
+$(SRC)/avperiodogram.$(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 40398)
+++ trunk/Ohana/src/opihi/dvo/avperiodogram.c	(revision 40398)
@@ -0,0 +1,398 @@
+# include "dvoshell.h"
+# define MIN_VAR 1e-8
+# define NPERIODS 50
+
+typedef struct {
+  int Nperiods;
+  opihi_flt *period;
+  opihi_flt *power;
+} PeriodogramResult;
+
+void PeriodogramResultFree (PeriodogramResult *result);
+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) {
+  
+  int N, next, Nfields;
+
+  dbStack *stack = NULL;
+  dbField *fields = NULL;
+  dbValue *values = NULL;
+  SkyList *skylist = NULL;
+  SkyRegionSelection *selection = NULL;
+
+  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;
+  }
+
+  if ((N = get_argument (argc, argv, "-min-period"))) {
+    remove_argument (N, &argc, argv);
+    minP = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-max-period"))) {
+    remove_argument (N, &argc, argv);
+    maxP = 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);
+
+  // XXX : I've removed -parallel-local : see avextract for an example
+
+  Catalog catalog;
+  dvo_catalog_init (&catalog, TRUE);
+
+  /* load photcode information */
+  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; 
+  }
+
+  // init locally static variables (time refs)
+  dbExtractAveragesInit (); 
+
+  // command-line is of the form: avperiodogram where (field op value)...
+
+  // examine line for 'where' or 'match to'.  'match to' is forbidden
+  int state = dbCmdlineConditions (argc, argv, 1, &next);
+  if (state == DVO_DB_CMDLINE_ERROR) goto escape;
+  if (state == DVO_DB_CMDLINE_IS_MATCH) goto escape; // not allowed for avperiodogram
+
+  // parse the remainder of the line as a boolean math expression
+  unsigned int Ncstack;
+  char **cstack = isolate_elements (argc-next, &argv[next], &Ncstack);
+  
+  // construct the db Boolean math stack (frees cstack)
+  int Nstack;
+  stack = dbRPN (Ncstack, cstack, &Nstack);
+  if (Ncstack && !Nstack) {
+    print_error(); 
+    goto escape; 
+  }
+
+  // add the skyregion limits to the where statement (or create)
+  dbAstroRegionLimits (&stack, &Nstack, selection, DVO_TABLE_AVERAGE);
+
+  // parse stack elements into fields and scalars as needed
+  if (!dbCheckStack (stack, Nstack, DVO_TABLE_AVERAGE, &fields, &Nfields)) goto escape;
+
+  // output values
+  ALLOCATE (values, dbValue, Nfields);
+
+  /* load region corresponding to selection above */
+  if ((skylist = SelectRegions (selection)) == NULL) goto escape;
+
+  // this does all the work of re-packaging the command, calling it on the remote machines, then loading in the results
+  if (PARALLEL && !HOST_ID) {
+    int status = HostTableParallelOps (skylist, argc, argv, RESULT_FILE, TRUE, 0, VERBOSE);
+
+    dbFreeFields (fields, Nfields);
+    dbFreeStack (stack, Nstack);
+    free (stack);
+    FreeSkyRegionSelection (selection);
+    dvo_catalog_free (&catalog);
+
+    return status;
+  }
+
+  // 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];
+    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);
+
+    // Scan the catalog for objects which match the WHERE clause
+    for (off_t j = 0; (j < catalog.Naverage) && !interrupt; j++) {
+      dbExtractAveragesInitAve ();  // reset counters for saved fields (costs very little)
+
+      Average *average = &catalog.average[j];
+
+      off_t mOff = average->measureOffset;
+      Measure *measure = &catalog.measure[mOff];
+
+      off_t nSec = j*Nsecfilt;
+      SecFilt *secfilt = &catalog.secfilt[nSec];
+
+      // extract the values needed to evaluate the WHERE clause
+      for (off_t n = 0; n < Nfields; n++) {
+	// for Measure, we are passing in the *first* measure, but average->Nmeasure gives the count
+	// for Secfilt, we are passing in the first entry; photcode.equiv gives the entry
+	// for Lensobj, we are passing in the first entry or NULL; photcode.equiv gives the entry
+	values[n] = dbExtractAverages (average, secfilt, measure, NULL, NULL, NULL, &fields[n]);
+      }
+
+      // test the conditional statement
+      if (!dbBooleanCond (stack, Nstack, values)) continue;
+
+      // 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 != j) {
+	  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, NULL, NULL, &mfields[n]);
+	}
+	if (!dbBooleanCond (mstack, Nmstack, mvalues)) continue;
+
+	time[Npts] = mvalues[0].Flt;
+	mag[Npts]  = mvalues[1].Flt;
+	dmag[Npts] = mvalues[2].Flt;
+
+	Npts++;
+      }
+
+      PeriodogramResult *result = periodogram_fm_raw (time, mag, dmag, Npts);
+      // XXX write the periodogram result to a file
+      PeriodogramResultFree (result);
+    }
+    dvo_catalog_free (&catalog);
+  }
+  ClearInterrupt (old_sigaction);
+
+  if (values) free (values);
+  dbFreeFields (fields, Nfields);
+  dbFreeStack (stack, Nstack);
+  if (stack) free (stack);
+  SkyListFree (skylist);
+  FreeSkyRegionSelection (selection);
+  return (TRUE);
+
+ escape:
+  if (values) free (values);
+  dbFreeFields (fields, Nfields);
+  dbFreeStack (stack, Nstack);
+  if (stack) free (stack);
+  SkyListFree (skylist);
+  FreeSkyRegionSelection (selection);
+  dvo_catalog_free (&catalog);
+  return (FALSE);
+
+ help:
+  gprint (GP_ERR, "USAGE: avperiodogram where (expression) -where-measure (expression)\n");
+  gprint (GP_ERR, "  the where (expression) is a boolean to limit the objects analysed based on average properties\n");
+  gprint (GP_ERR, "  NOTE: the optional -where-measure (expression) MUST come after the average where expression\n");
+  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;
+}
Index: trunk/Ohana/src/opihi/dvo/init.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/init.c	(revision 40397)
+++ trunk/Ohana/src/opihi/dvo/init.c	(revision 40398)
@@ -3,4 +3,5 @@
 int avextract       PROTO((int, char **));
 int avmatch         PROTO((int, char **));
+int avperiodogram   PROTO((int, char **));
 int badimages       PROTO((int, char **));
 int calextract      PROTO((int, char **));
@@ -67,4 +68,5 @@
   {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, "badimages",   badimages,    "look for images with anomalous astrometry"},
 //  {1, "calextract",  calextract,   "extract photometry calibration"},
Index: trunk/Ohana/src/opihi/lib.data/starfuncs.c
===================================================================
--- trunk/Ohana/src/opihi/lib.data/starfuncs.c	(revision 40397)
+++ trunk/Ohana/src/opihi/lib.data/starfuncs.c	(revision 40398)
@@ -285,4 +285,5 @@
     off = j*Nx;
     for (i = Xs; i < Xe; i++) { 
+      if (!isfinite(data[i+off])) continue;
       sky[Nsky] = data[i+off];
       Nsky ++;
@@ -294,4 +295,5 @@
     off = j*Nx;
     for (i = Xs; i < Xe; i++) { 
+      if (!isfinite(data[i+off])) continue;
       sky[Nsky] = data[i+off];
       Nsky ++;
@@ -305,4 +307,5 @@
     off = j*Nx;
     for (i = Xs; i < Xe; i++) { 
+      if (!isfinite(data[i+off])) continue;
       sky[Nsky] = data[i+off];
       Nsky ++;
@@ -314,4 +317,5 @@
     off = j*Nx;
     for (i = Xs; i < Xe; i++) { 
+      if (!isfinite(data[i+off])) continue;
       sky[Nsky] = data[i+off];
       Nsky ++;
@@ -341,4 +345,5 @@
       rad2 = SQ(Xc) + SQ(Yc);
       if (rad2 > Ro2) continue;
+      // if (!isfinite(data[i+off])) continue;
       value = data[i+off] - fsky;
       Sx  += Xc*value;
Index: trunk/Ohana/src/opihi/mana/findrowpeaks.c
===================================================================
--- trunk/Ohana/src/opihi/mana/findrowpeaks.c	(revision 40397)
+++ trunk/Ohana/src/opihi/mana/findrowpeaks.c	(revision 40398)
@@ -23,4 +23,5 @@
 
   for (i = 1; i < Nrow - 1; i++) {
+    if (!isfinite(row[i])) continue; // ignore NAN values
     if (row[i] < threshold) continue;
     if (row[i] < row[i-1]) continue;