Changeset 39558 for trunk/Ohana/src/opihi/cmd.astro

Timestamp:

Apr 28, 2016, 11:09:46 AM (10 years ago)

Author:

eugene

Message:

add datan2, set atan2 to return radians as it should

Location:

trunk/Ohana/src/opihi/cmd.astro

Files:

• specpairfit.c (modified) (2 diffs)
• transform.c (modified) (1 diff)

trunk/Ohana/src/opihi/cmd.astro/specpairfit.c

@@ -19 +19 @@
   CastVector (window, OPIHI_INT);
 
-  // minimize (flux1 - flux2*alpha) in window defined by mask
+  // minimize (flux2 - flux1*A - B) in window defined by mask
   // note that the mask is a SELECTION mask not an EXCLUSION mask
 
-  double F12 = 0.0; 
-  double F22 = 0.0;
+  double R = 0.0, F1 = 0.0, F11 = 0.0, F12 = 0.0, F2 = 0.0;
   for (i = 0; i < flux1->Nelements; i++) {
-    if ((mask & window->elements.Int[i]) == 0) continue;
+    // if ((mask & window->elements.Int[i]) == 0) continue;
     double weight = 1.0 / (SQ(dflux1->elements.Flt[i]) + SQ(dflux2->elements.Flt[i]));
+    R   += weight;
+    F1  += flux1->elements.Flt[i] * weight;
+    F2  += flux2->elements.Flt[i] * weight;
+    F11 += flux1->elements.Flt[i] * flux1->elements.Flt[i] * weight;
     F12 += flux1->elements.Flt[i] * flux2->elements.Flt[i] * weight;
-    F22 += flux2->elements.Flt[i] * flux2->elements.Flt[i] * weight;
   }
 
-  double Ao = F12 / F22;
-  double dA = sqrt(1.0 / F22);
+  int nterm = 2;
+  double **b = NULL, **c = NULL;
+  ALLOCATE (b, double *, nterm);
+  ALLOCATE (c, double *, nterm);
+  for (i = 0; i < nterm; i++) {
+    ALLOCATE (c[i], double, nterm);
+    ALLOCATE (b[i], double, 1);
+  }
+  c[0][0] = F11;
+  c[1][0] = c[0][1] = F1;
+  c[1][1] = R;
 
-  int Ndof = -1; // 1 parameter fit
+  b[0][0] = F12;
+  b[1][0] = F2;
+
+  if (!dgaussjordan (c, b, nterm, 1)) {
+    gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n");
+    goto escape;
+  }
+
+  double Ao = b[0][0];
+  double dA = sqrt(c[0][0]);
+  double Bo = b[1][0];
+  double dB = sqrt(c[1][1]);
+
+  for (i = 0; i < nterm; i++) {
+    free (b[i]);
+    free (c[i]);
+  }
+  free (b);
+  free (c);
+
+  int Ndof = -2; // 2 parameter fit
   double chisq = 0.0;
   for (i = 0; i < flux1->Nelements; i++) {
     if ((mask & window->elements.Int[i]) == 0) continue;
     double weight = 1.0 / (SQ(dflux1->elements.Flt[i]) + SQ(dflux2->elements.Flt[i]));
-    chisq += SQ(flux1->elements.Flt[i] - Ao * flux2->elements.Flt[i]) * weight;
+    chisq += SQ(flux1->elements.Flt[i] - Ao * flux2->elements.Flt[i] - Bo) * weight;
     Ndof ++;
   }
@@ -45 +76 @@
   double chisqNu = chisq / Ndof;
 
+  // fprintf (stderr, "R: %f, F1: %f, F11: %f, F2: %f, F12: %f\n", R, F1, F11, F2, F12);
   // fprintf (stderr, "Ao: %f +/- %f, chisq: %f, chisq_nu : %f for %d dof\n", Ao, dA, chisq, chisqNu, Ndof);
   set_variable ("Ao", Ao);
   set_variable ("dA", dA);
+  set_variable ("Bo", Bo);
+  set_variable ("dB", dB);
   set_variable ("Xv", chisqNu);
   set_variable ("Nd", Ndof);

trunk/Ohana/src/opihi/cmd.astro/transform.c

@@ -53 +53 @@
           X = x; Y = y;
           if ((X > -1) && (X < Nx) && (Y > -1) && (Y < Ny)) {
+            if (!isfinite(*Vin)) continue;
             Vout[X + Y*Nx] += *Vin;
             Sout[X + Y*Nx] ++;