Changeset 40145

Timestamp:

Sep 26, 2017, 6:22:26 AM (9 years ago)

Author:

eugene

Message:

working on periodogram-fm.c

File:

• branches/eam_branches/ohana.20170822/src/opihi/cmd.data/periodogram-fm.c (modified) (4 diffs)

branches/eam_branches/ohana.20170822/src/opihi/cmd.data/periodogram-fm.c

@@ -8 +8 @@
   float mean, var, w, tau, P, Pc, Ps, Po;
   float C, S, cs, sn, cs2, sn2, ratio;
-  Vector *time, *flux, *power, *period;
+  Vector *time, *flux, *dflux, *power, *period;
 
   VERBOSE = FALSE;
@@ -16 +16 @@
   }
 
-  if (argc != 7) {
-    gprint (GP_ERR, "USAGE: periodogram (time) (flux) (minP) (maxP) (period) (power)\n");
+  if (argc != 8) {
+    gprint (GP_ERR, "USAGE: periodogram-fm (time) (flux) (dflux) (minP) (maxP) (period) (power)\n");
     return (FALSE);
   }
   
-  if ((time = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((flux = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  minP = atof(argv[3]);
-  maxP = atof(argv[4]);
-  if ((period = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((power = SelectVector (argv[6], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  // XXX allow dflux to be dropped?
+
+  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);
+  minP = atof(argv[4]);
+  maxP = atof(argv[5]);
+  if ((period = SelectVector (argv[6], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((power  = SelectVector (argv[7], ANYVECTOR, TRUE)) == NULL) return (FALSE);
 
   REQUIRE_VECTOR_FLT (time, FALSE); 
@@ -62 +65 @@
   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
   P = minP;
   while (P < maxP) {
@@ -68 +76 @@
     /* find the period offset tau  */
     tv = time[0].elements.Flt;
-    cs = sn = 0;
+    df = dflux[0].elements.Flt;
+
+    cs = sn = cs2 = sn2 = 0;
+
     for (i = 0; i < Npt; i++, tv++) {
+      opihi_flt wt = 
       cs += cos (*tv*w*2);
       sn += sin (*tv*2*2);