Changeset 36150

Timestamp:

Sep 24, 2013, 12:16:12 PM (13 years ago)

Author:

eugene

Message:

handle invalid weights in plx fits, list data for gauss-jordan inversion failures

File:

• branches/eam_branches/ipp-20130904/Ohana/src/opihi/cmd.astro/fitplx.c (modified) (5 diffs)

branches/eam_branches/ipp-20130904/Ohana/src/opihi/cmd.astro/fitplx.c

@@ -41 +41 @@
     remove_argument (N, &argc, argv);
     VERBOSE = TRUE;
+  }
+  if ((N = get_argument (argc, argv, "-vv"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = 2;
   }
 
@@ -145 +149 @@
 
   PlxFit fit;
-  FitPMandPar (&fit, X, dX, Y, dY, t, pX, pY, n, VERBOSE);
+  if (!FitPMandPar (&fit, X, dX, Y, dY, t, pX, pY, n, VERBOSE)) {
+    return FALSE;
+  }
 
   // fprintf (stderr, "Roff, Doff: %f, %f; dRo, dDo: %f, %f\n", fit.Ro, fit.Do, fit.dRo, fit.dDo);
@@ -204 +210 @@
   for (i = 0; i < Npts; i++) {
 
-    if (VERBOSE) fprintf (stderr, "%f %f : %f %f : %f : %f %f\n", X[i], dX[i], Y[i], dY[i], T[i], pR[i], pD[i]);
+    if (VERBOSE == 2) fprintf (stderr, "%f %f : %f %f : %f : %f %f\n", X[i], dX[i], Y[i], dY[i], T[i], pR[i], pD[i]);
 
     /* handle case where dX or dY = 0.0 */
-    wx = 1.0 / SQ(dX[i]);
-    wy = 1.0 / SQ(dY[i]);
+    wx = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
+    wy = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
 
     Wx += wx;
@@ -266 +272 @@
   B[4][0] = PRX + PDY;
 
-  dgaussjordan ((double **)A, (double **)B, 5, 1);
+  if (!dgaussjordan ((double **)A, (double **)B, 5, 1)) {
+    if (VERBOSE) fprintf (stderr, "error in fit\n");
+    if (VERBOSE == 2) {
+      int j;
+      for (i = 0; i < 5; i++) {
+        for (j = 0; j < 5; j++) {
+          fprintf (stderr, "%e ", A[i][j]);
+        }
+        fprintf (stderr, " : %e\n", A[i][0]);
+      }
+    }
+    return FALSE;
+  }
 
   fit[0].Ro = B[0][0];
@@ -285 +303 @@
     Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
     Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i];
-    chisq += SQ(X[i] - Xf) / SQ(dX[i]);
-    chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
-    if (VERBOSE) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
+    wx = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
+    wy = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
+    chisq += SQ(X[i] - Xf) * wx;
+    chisq += SQ(Y[i] - Yf) * wy;
+    // if (VERBOSE) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
   }
   fit[0].Nfit = Npts;