Index: /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram-fm.sh
===================================================================
--- /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram-fm.sh	(revision 40218)
+++ /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram-fm.sh	(revision 40219)
@@ -55,5 +55,4 @@
 
  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
@@ -88,5 +87,4 @@
 
  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
@@ -97,4 +95,8 @@
  if (abs ($peakpos - $P) > 0.05)
    $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
@@ -117,5 +119,4 @@
 
  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
@@ -126,4 +127,8 @@
  if (abs ($peakpos - $P) > 0.05)
    $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
@@ -145,5 +150,5 @@
  set df = 0.01 + zero(f)
 
- periodogram_fm t f 2 30 period power
+ periodogram_fm t f df 2 30 period power
 
 #  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
@@ -155,7 +160,11 @@
    $PASS = 0
  end
-end
-
-# test using random samples, offset start, non-zero DC, some noise
+
+ 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 300 random samples, offset start, non-zero DC, some noise
 macro test6
  $PASS = 1
@@ -178,5 +187,5 @@
  set df = 0.01 + zero(f)
 
- periodogram_fm t f 2 30 period power
+ periodogram_fm t f df 2 30 period power
 
 #  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
@@ -188,7 +197,11 @@
    $PASS = 0
  end
-end
-
-# test using fewer random samples, offset start, non-zero DC, some noise
+
+ 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 100 fewer random samples, offset start, non-zero DC, some noise
 macro test7
  $PASS = 1
@@ -211,5 +224,5 @@
  set df = 0.01 + zero(f)
 
- periodogram_fm t f 2 30 period power
+ periodogram_fm t f df 2 30 period power
 
 #  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
@@ -221,22 +234,28 @@
    $PASS = 0
  end
-end
-
-# test using fewer random samples, high frequency, non-zero DC, some noise
+
+ 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 Ndays random samples, RR Lyrae-sized light curves (0.7 mag),
+# optional noise level 
 macro test8
- if ($0 != 2)
-   echo "USAGE: test8: Ndays")
+ if ($0 != 4)
+   echo "USAGE: test8: Period Ndays (df)"
    break
  end
  
  local Ndays 
- $Ndays = $1
-
- $PASS = 1
- break -auto off
-
- local P PI
- $PI = 3.14159265359
- $P  = 0.8*rnd(0) + 0.2
+ $P = $1
+ $Ndays = $2
+ $dM = $3
+
+ $PASS = 1
+ break -auto off
+
+ local PI
+ $PI = 3.14159265359
  $trueP = $P
 
@@ -246,7 +265,7 @@
 
  # 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 tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100
+ set dtx = (3/24) * rnd(x);  # choose a starting time within that night
+ set t0 = tday + dtx
 
  set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440
@@ -260,13 +279,13 @@
  set tmp = t0 + dt3; concat tmp t
 
- set fraw = sin(2*$PI*t/$P) + 0.5
+ set fraw = 0.75*sin(2*$PI*t/$P)
 
  # 0.05 : peakpos = 14.95
  # 0.10 : peakpos = 15.04 (
- gaussdev df t[] 0.0 0.25
+ gaussdev df t[] 0.0 $dM
  set f = fraw + df
- set df = 0.01 + zero(f)
-
- periodogram_fm t f 0.1 2.0 period power
+ set df = $dM + zero(f)
+
+ periodogram_fm t f df 0.1 20.0 period power
 
 #  lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f
@@ -278,5 +297,233 @@
    $PASS = 0
  end
-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 Ndays random samples, RR Lyrae-sized light curves (0.7 mag),
+# optional noise level 
+# compare periodogram and periodogram_fm
+macro test9
+ if ($0 != 4)
+   echo "USAGE: test8: Period Ndays (df)"
+   break
+ end
+ 
+ local Ndays 
+ $P = $1
+ $Ndays = $2
+ $dM = $3
+
+ $PASS = 1
+ break -auto off
+
+ local PI
+ $PI = 3.14159265359
+ $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 tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100
+ set dtx = (3/24) * rnd(x);  # choose a starting time within that night
+ set t0 = tday + 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 = 0.75*sin(2*$PI*t/$P)
+
+ # 0.05 : peakpos = 14.95
+ # 0.10 : peakpos = 15.04 (
+ gaussdev df t[] 0.0 $dM
+ set f = fraw + df
+ set df = $dM + zero(f)
+
+ periodogram_fm t f df 0.1 20.0 period_fm power_fm
+ periodogram t f 0.1 20.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_fm power_fm
+ $peakval_fm = $peakval
+
+ peak -q period power
+ vstat -q power
+ set power = power / $MAX
+
+# if (abs ($peakpos - $P) > 0.05)
+#   $PASS = 0
+# end
+
+ set freq = 1 / period
+ set freq_fm = 1 / period_fm
+ $Freq = 1 / $P
+
+ if ($PLOT)
+  if (1)
+    lim period power; clear; box
+    line -c red70 -lw 3 $P 0 to $P $peakval_fm;
+    plot period power -x line -c grey70 -lw 2
+    plot period_fm power_fm -x line -c black
+  else
+    lim freq power; clear; box
+    line -c red70 -lw 3 $Freq 0 to $Freq 1.0
+    plot freq power -x line -c grey70 -lw 2
+    plot freq_fm power_fm -x line -c black
+  end
+ end
+end
+
+
+# test using Ndays random samples, RR Lyrae-sized light curves (0.7 mag),
+# optional noise level 
+# compare periodogram and periodogram_fm
+macro test10
+ if ($0 != 4)
+   echo "USAGE: test8: Period Ndays (df)"
+   break
+ end
+ 
+ local Ndays 
+ $P = $1
+ $Ndays = $2
+ $dM = $3
+
+ $PASS = 1
+ break -auto off
+
+ local PI
+ $PI = 3.14159265359
+ $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 tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100
+ set dtx = (3/24) * rnd(x);  # choose a starting time within that night
+ set t0 = tday + 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 = 0.75*sin(2*$PI*t/$P)
+
+ # 0.05 : peakpos = 14.95
+ # 0.10 : peakpos = 15.04 (
+ gaussdev df t[] 0.0 $dM
+ set f = fraw + df
+ set df = $dM + zero(f)
+
+ periodogram_fm t f df 0.05 20.0 period_fm power_fm
+
+ gaussdev df t[] 0.0 $dM
+ set Fo = df
+ periodogram_fm t Fo df 0.05 20.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_fm power_fm
+ $peakval_fm = $peakval
+
+ peak -q period power
+
+# if (abs ($peakpos - $P) > 0.05)
+#   $PASS = 0
+# end
+
+ set freq = 1 / period
+ set freq_fm = 1 / period_fm
+ $Freq = 1 / $P
+
+ if ($PLOT)
+  if (1)
+    lim period_fm power_fm; clear; box
+    line -c red70 -lw 3 $P 0 to $P $peakval_fm;
+    plot period power -x line -c grey70 -lw 2
+    plot period_fm power_fm -x line -c black
+  else
+    lim freq power; clear; box
+    line -c red70 -lw 3 $Freq 0 to $Freq 1.0
+    plot freq power -x line -c grey70 -lw 2
+    plot freq_fm power_fm -x line -c black
+  end
+ end
+end
+
+# we have time (MJD) and mag
+# we generate the folded lightcure and measure sigma relative to the smoothed version (bins of 0.1 period)
+macro fold.one.period
+  if ($0 != 5)
+    echo "USAGE: fold.one.period (time) (mag) (magErr) (period)"
+    break
+  end
+
+  local myTime myMag myMagErr myPeriod
+  $myTime = $1
+  $myMag  = $2
+  $myMagErr  = $3
+  $myPeriod = $4
+
+  set phi = $myTime / $myPeriod - int($myTime / $myPeriod)
+
+  if ($PLOT_FOLD)
+    lim -n phi phi $myMag; clear; box; 
+  end
+
+  delete -q magResid
+
+  $dPhi = 0.05; # half of bin size
+  create nphi $dPhi {1 + $dPhi} {2*$dPhi}
+  set magR = zero(nphi)
+  set magS = zero(nphi)
+  for i 0 nphi[]
+    subset tmp_mag_sub = $myMag where (phi >= nphi[$i] - $dPhi) && (phi < nphi[$i] + $dPhi)
+    vstat -q tmp_mag_sub
+    magR[$i] = $MEDIAN
+    magS[$i] = $SIGMA
+
+    set magDelta = tmp_mag_sub - $MEDIAN
+    concat magDelta magResid 
+
+    if ($PLOT_FOLD)
+      subset tmp_phi_sub = phi where (phi >= nphi[$i] - $dPhi) && (phi < nphi[$i] + $dPhi)
+      if ($i % 2)
+        plot tmp_phi_sub tmp_mag_sub -pt 7 -sz 3 -c blue -lw 2
+      else
+        plot tmp_phi_sub tmp_mag_sub -pt 7 -sz 3 -c red -lw 2
+      end
+    end  
+  end
+
+  if ($PLOT_FOLD)  
+    plot -pt 10 -sz 1.5 phi $myMag -dy $myMagErr
+    plot -pt 2 -sz 2.0 -c red nphi magR -dy magS
+  end
+
+  vstat -q magResid
+end
+
+
 
 # Memory test
@@ -298,5 +545,5 @@
 
  for i 0 100
-  periodogram_fm t f 2 30 period power
+  periodogram_fm t f df 2 30 period power
  end
   
Index: /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram.sh
===================================================================
--- /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram.sh	(revision 40218)
+++ /branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram.sh	(revision 40219)
@@ -209,18 +209,19 @@
 # test using fewer random samples, high frequency, non-zero DC, some noise
 macro test8
- if ($0 != 2)
-   echo "USAGE: test8: Ndays")
+ if ($0 != 4)
+   echo "USAGE: test8: Period Ndays df"
    break
  end
  
  local Ndays 
- $Ndays = $1
-
- $PASS = 1
- break -auto off
-
- local P PI
- $PI = 3.14159265359
- $P  = 0.8*rnd(0) + 0.2
+ $P = $1
+ $Ndays = $2
+ $dM = $3
+
+ $PASS = 1
+ break -auto off
+
+ local PI
+ $PI = 3.14159265359
  $trueP = $P
 
@@ -230,7 +231,7 @@
 
  # 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 tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100
+ set dtx = (3/24) * rnd(x);  # choose a starting time within that night
+ set t0 = tday + dtx
 
  set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440
@@ -244,9 +245,9 @@
  set tmp = t0 + dt3; concat tmp t
 
- set fraw = sin(2*$PI*t/$P) + 0.5
+ set fraw = 0.75*sin(2*$PI*t/$P)
 
  # 0.05 : peakpos = 14.95
  # 0.10 : peakpos = 15.04 (
- gaussdev df t[] 0.0 0.25
+ gaussdev df t[] 0.0 $dM
  set f = fraw + df
 
@@ -260,4 +261,7 @@
  if (abs ($peakpos - $P) > 0.05)
    $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