Index: trunk/Ohana/src/opihi/cmd.data/Makefile
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 40015)
+++ trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 40165)
@@ -102,4 +102,5 @@
 $(SRC)/peak.$(ARCH).o		\
 $(SRC)/periodogram.$(ARCH).o	\
+$(SRC)/periodogram-fm.$(ARCH).o	\
 $(SRC)/plot.$(ARCH).o		\
 $(SRC)/dot.$(ARCH).o		\
Index: trunk/Ohana/src/opihi/cmd.data/dot.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/dot.c	(revision 40015)
+++ trunk/Ohana/src/opihi/cmd.data/dot.c	(revision 40165)
@@ -29,5 +29,5 @@
 
   /* set point style and errorbar mode (these are NOT sticky) */
-  graphmode.style = 2;
+  graphmode.style = KAPA_PLOT_POINTS;
   graphmode.etype = 0;
 
Index: trunk/Ohana/src/opihi/cmd.data/grid.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/grid.c	(revision 40015)
+++ trunk/Ohana/src/opihi/cmd.data/grid.c	(revision 40165)
@@ -175,6 +175,6 @@
 
   Xvec.Nelements = Yvec.Nelements = N;
-  graphmode.style = 2; /* points */
-  graphmode.ptype = 100; /* connect a pair */
+  graphmode.style = KAPA_PLOT_POINTS; /* points */
+  graphmode.ptype = KAPA_POINT_PAIR_CONNECT; /* connect pairs of points */
   graphmode.etype = 0;
   PlotVectorPair (kapa, &Xvec, &Yvec, NULL, &graphmode);
Index: trunk/Ohana/src/opihi/cmd.data/init.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/init.c	(revision 40015)
+++ trunk/Ohana/src/opihi/cmd.data/init.c	(revision 40165)
@@ -91,4 +91,5 @@
 int peak             PROTO((int, char **));
 int periodogram      PROTO((int, char **));
+int periodogram_fm   PROTO((int, char **));
 int plot             PROTO((int, char **));
 int dot              PROTO((int, char **));
@@ -270,5 +271,6 @@
   {1, "parity",       parity,           "set image parity"},
   {1, "peak",         peak,             "find vector peak in range"},
-  {1, "periodogram",  periodogram,      "measure periods in unevenly sampled data"},
+  {1, "periodogram",    periodogram,    "measure periods in unevenly sampled data (Lomb-Scargle)"},
+  {1, "periodogram_fm", periodogram_fm, "measure periods in unevenly sampled data (generalized Lomb-Scargle; floating mean)"},
   {1, "plot",         plot,             "plot a pair of vectors"},
   {1, "png",          jpeg,             "convert display graphic to PNG"},
Index: trunk/Ohana/src/opihi/cmd.data/line.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/line.c	(revision 40015)
+++ trunk/Ohana/src/opihi/cmd.data/line.c	(revision 40165)
@@ -33,5 +33,5 @@
 
   /* set point style and errorbar mode (these are NOT sticky) */
-  graphmode.style = 0;
+  graphmode.style = KAPA_PLOT_CONNECT;
   graphmode.etype = 0;
 
Index: trunk/Ohana/src/opihi/cmd.data/periodogram-fm.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/periodogram-fm.c	(revision 40165)
+++ trunk/Ohana/src/opihi/cmd.data/periodogram-fm.c	(revision 40165)
@@ -0,0 +1,146 @@
+# include "data.h"
+# define MIN_VAR 1e-8
+int periodogram_fm (int argc, char **argv) {
+  
+  int N;
+
+  int VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    VERBOSE = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 8) {
+    gprint (GP_ERR, "USAGE: periodogram_fm (time) (flux) (dflux) (minP) (maxP) (period) (power)\n");
+    return (FALSE);
+  }
+  
+  // XXX allow dflux to be dropped?
+
+  Vector *time, *flux, *dflux, *power, *period;
+  if ((time   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((flux   = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((dflux  = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((period = SelectVector (argv[6], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((power  = SelectVector (argv[7], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  double minP = atof(argv[4]);
+  double maxP = atof(argv[5]);
+
+  REQUIRE_VECTOR_FLT (time, FALSE); 
+  REQUIRE_VECTOR_FLT (flux, FALSE); 
+  REQUIRE_VECTOR_FLT (dflux, FALSE); 
+
+  /* find the max baseline, sum the inverse variances */
+  double minT = time[0].elements.Flt[0];
+  double maxT = time[0].elements.Flt[0];
+  opihi_flt Weight = 0;
+  int Npt = time[0].Nelements;
+
+  opihi_flt *tv =  time[0].elements.Flt;
+  opihi_flt *df = dflux[0].elements.Flt;
+  for (int i = 0; i < Npt; 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 (FALSE);
+  }
+
+  int Np = 0;
+  int NP = 100;
+  ResetVector (power,  OPIHI_FLT, NP);
+  ResetVector (period, OPIHI_FLT, NP);
+
+  // for testing, we are going to write out the terms explicitly 
+  // for optimization, some of the tri functions can be calculated more efficiently 
+  // by storing cos,sin(w t) and using some recurrence rules
+
+  // XXX consider choice of period or frequency step
+  double P = minP;
+  while (P <= maxP) {
+    double w = 2*M_PI/P;
+    
+    /* find the period offset tau  */
+    tv = time[0].elements.Flt;
+    df = dflux[0].elements.Flt;
+
+    double cs = 0.0, sn = 0.0, sn2 = 0.0, cs2 = 0.0;
+
+    for (int i = 0; i < Npt; 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[0].elements.Flt;
+    df = dflux[0].elements.Flt;
+    opihi_flt *fv = flux[0].elements.Flt;
+
+    // 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 < Npt; 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[0].elements.Flt[Np] = Po;
+    period[0].elements.Flt[Np] = P;
+    Np ++;
+    if (Np >= NP) {
+      NP += 100;
+      REALLOCATE (power[0].elements.Flt, opihi_flt, NP);
+      REALLOCATE (period[0].elements.Flt, opihi_flt, NP);
+    }
+
+    double ratio = 1 + 0.1*P/dTime;
+
+    if (VERBOSE) gprint (GP_ERR, "tau: %f, P: %f, ratio: %f, dTime: %f, nextP: %f\n", tau, P, ratio, dTime, P*ratio);
+
+    P *= ratio;
+  }
+
+  ResetVector (power,  OPIHI_FLT, Np);
+  ResetVector (period, OPIHI_FLT, Np);
+ 
+  return (TRUE);
+}
Index: trunk/Ohana/src/opihi/cmd.data/test/periodogram-fm.sh
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/test/periodogram-fm.sh	(revision 40165)
+++ trunk/Ohana/src/opihi/cmd.data/test/periodogram-fm.sh	(revision 40165)
@@ -0,0 +1,313 @@
+
+if (not($?PLOT)) set PLOT = 0
+
+list tests
+ test1
+ test2
+ test3
+ memtest1
+end
+
+# test using even samples
+macro test1
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q t f period power
+
+ create t 0 100
+ set f = sin(2*$PI*t/$P)
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f df 5 50 period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.5)
+   $PASS = 0
+   echo "OFFSET: {$peakpos - $P}"
+ end
+
+ if ($PLOT)
+  lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line
+ end
+end
+
+# test using random samples
+macro test2
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 0 100
+ set t = 100 * rnd(x) 
+ set f = sin(2*$PI*t/$P)
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f df 5 50 period power
+#periodogram t f 5 50 period power
+
+ # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+ # lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.5)
+   $PASS = 0
+ end
+
+ if ($PLOT)
+  lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line
+ end
+end
+
+# test using random samples, higher frequency
+macro test3
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 2.0
+
+ delete -q x t f period power
+
+ create x 0 100
+ set t = 100 * rnd(x) 
+ set f = sin(2*$PI*t/$P)
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f df 1 10 period power
+#periodogram t f 1 10 period power
+
+ # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+ # lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using random samples, offset start
+macro test4
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 500 800
+ set t = 300 * rnd(x) + 500
+ set f = sin(2*$PI*t/$P)
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f df 2 30 period power
+#periodogram t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using random samples, offset start, non-zero DC
+macro test5
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 500 800
+ set t = 300 * rnd(x) + 500
+ set f = sin(2*$PI*t/$P) + 0.5
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using random samples, offset start, non-zero DC, some noise
+macro test6
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 500 800
+ set t = 300 * rnd(x) + 500
+ set fraw = sin(2*$PI*t/$P) + 0.5
+
+ # 0.05 : peakpos = 14.95
+ # 0.10 : peakpos = 15.04 (
+ gaussdev df t[] 0.0 0.25
+ set f = fraw + df
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using fewer random samples, offset start, non-zero DC, some noise
+macro test7
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 0 100
+ set t = 100 * rnd(x)
+ set fraw = sin(2*$PI*t/$P) + 0.5
+
+ # 0.05 : peakpos = 14.95
+ # 0.10 : peakpos = 15.04 (
+ gaussdev df t[] 0.0 0.25
+ set f = fraw + df
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using fewer random samples, high frequency, non-zero DC, some noise
+macro test8
+ if ($0 != 2)
+   echo "USAGE: test8: Ndays")
+   break
+ end
+ 
+ local Ndays 
+ $Ndays = $1
+
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 0.8*rnd(0) + 0.2
+ $trueP = $P
+
+ delete -q x t f period power
+
+ create x 0 $Ndays
+
+ # t is a time in days, but we always have 4 within 1 hour:
+ set t0 = int(100 * rnd(x))
+ set dtx = (3/24) * rnd(x)
+ set t0 = t0 + dtx
+
+ set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440
+ set dt2 = (15.0 / 1440) * rnd(x) + (15 + 7.5) / 1440
+ set dt3 = (15.0 / 1440) * rnd(x) + (30 + 7.5) / 1440
+
+ delete -q t
+ concat t0 t
+ set tmp = t0 + dt1; concat tmp t
+ set tmp = t0 + dt2; concat tmp t
+ set tmp = t0 + dt3; concat tmp t
+
+ set fraw = sin(2*$PI*t/$P) + 0.5
+
+ # 0.05 : peakpos = 14.95
+ # 0.10 : peakpos = 15.04 (
+ gaussdev df t[] 0.0 0.25
+ set f = fraw + df
+ set df = 0.01 + zero(f)
+
+ periodogram_fm t f 0.1 2.0 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# Memory test
+macro memtest1
+
+ local i
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 500 800
+ set t = 300 * rnd(x) + 500
+ set f = sin(2*$PI*t/$P)
+
+ list word -x "ps -p $PID -o rss"
+ $startmem = $word:1
+
+ for i 0 100
+  periodogram_fm t f 2 30 period power
+ end
+  
+ list word -x "ps -p $PID -o rss"
+ $endmem = $word:1
+
+ $PASS = 1
+
+ if ($endmem - $startmem > 180)
+   $PASS = 0
+   echo "growth: {$endmem-$startmem}"
+   echo "kB/loop: {($endmem-$startmem)/100}"
+ end
+end
Index: trunk/Ohana/src/opihi/cmd.data/test/periodogram.sh
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/test/periodogram.sh	(revision 40015)
+++ trunk/Ohana/src/opihi/cmd.data/test/periodogram.sh	(revision 40165)
@@ -1,2 +1,4 @@
+
+if (not($?PLOT)) set PLOT = 0
 
 list tests
@@ -28,4 +30,8 @@
    echo "OFFSET: {$peakpos - $P}"
  end
+
+ if ($PLOT)
+  lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line
+ end
 end
 
@@ -99,4 +105,151 @@
 
  periodogram t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using random samples, offset start, non-zero DC
+macro test5
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 500 800
+ set t = 300 * rnd(x) + 500
+ set f = sin(2*$PI*t/$P) + 0.5
+
+ periodogram t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using random samples, offset start, non-zero DC, some noise
+macro test6
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 500 800
+ set t = 300 * rnd(x) + 500
+ set fraw = sin(2*$PI*t/$P) + 0.5
+
+ # 0.05 : peakpos = 14.95
+ # 0.10 : peakpos = 15.04 (
+ gaussdev df t[] 0.0 0.25
+ set f = fraw + df
+
+ periodogram t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using fewer random samples, offset start, non-zero DC, some noise
+macro test7
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 15.0
+
+ delete -q x t f period power
+
+ create x 0 100
+ set t = 100 * rnd(x)
+ set fraw = sin(2*$PI*t/$P) + 0.5
+
+ # 0.05 : peakpos = 14.95
+ # 0.10 : peakpos = 15.04 (
+ gaussdev df t[] 0.0 0.25
+ set f = fraw + df
+
+ periodogram t f 2 30 period power
+
+#  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
+#  lim -n 1 period power; clear; box; plot period power
+
+ peak -q period power
+
+ if (abs ($peakpos - $P) > 0.05)
+   $PASS = 0
+ end
+end
+
+# test using fewer random samples, high frequency, non-zero DC, some noise
+macro test8
+ if ($0 != 2)
+   echo "USAGE: test8: Ndays")
+   break
+ end
+ 
+ local Ndays 
+ $Ndays = $1
+
+ $PASS = 1
+ break -auto off
+
+ local P PI
+ $PI = 3.14159265359
+ $P  = 0.8*rnd(0) + 0.2
+ $trueP = $P
+
+ delete -q x t f period power
+
+ create x 0 $Ndays
+
+ # t is a time in days, but we always have 4 within 1 hour:
+ set t0 = int(100 * rnd(x))
+ set dtx = (3/24) * rnd(x)
+ set t0 = t0 + dtx
+
+ set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440
+ set dt2 = (15.0 / 1440) * rnd(x) + (15 + 7.5) / 1440
+ set dt3 = (15.0 / 1440) * rnd(x) + (30 + 7.5) / 1440
+
+ delete -q t
+ concat t0 t
+ set tmp = t0 + dt1; concat tmp t
+ set tmp = t0 + dt2; concat tmp t
+ set tmp = t0 + dt3; concat tmp t
+
+ set fraw = sin(2*$PI*t/$P) + 0.5
+
+ # 0.05 : peakpos = 14.95
+ # 0.10 : peakpos = 15.04 (
+ gaussdev df t[] 0.0 0.25
+ set f = fraw + df
+
+ periodogram t f 0.1 2.0 period power
 
 #  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
Index: trunk/Ohana/src/opihi/cmd.data/zplot.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/zplot.c	(revision 40015)
+++ trunk/Ohana/src/opihi/cmd.data/zplot.c	(revision 40165)
@@ -117,5 +117,5 @@
 
   /* point size determined by Zvec */
-  graphmode.style = 2; /* plot points */
+  graphmode.style = KAPA_PLOT_POINTS; /* points */
   graphmode.size = -1; /* point size determined by Zvec */
   PlotVectorTriplet (kapa, xvec, yvec, &Zvec, mask, &graphmode);
@@ -226,5 +226,5 @@
 
   /* point size determined by Zvec */
-  graphmode.style = 2; /* plot points */
+  graphmode.style = KAPA_PLOT_POINTS; /* plot points */
   graphmode.color = -1; /* point color determined by Zvec */
   graphmode.etype = 0; /* no errorbars */