Changeset 27548

Timestamp:

Mar 31, 2010, 7:03:35 PM (16 years ago)

Author:

eugene

Message:

calculate all subsets of fits

File:

• branches/eam_branches/relastro.20100326/src/UpdateObjects.c (modified) (8 diffs)

branches/eam_branches/relastro.20100326/src/UpdateObjects.c

@@ -44 +44 @@
   StatType statsR, statsD;
   Coords coords;
-  PMFit fit;
+  PMFit fitAve, fitPM, fitPar;
   time_t To;
   off_t Nave, Npm, Npar, Nskip;
@@ -64 +64 @@
 
   // use J2000 as a reference time
-  To = ohana_date_to_sec ("2000/01/01");
+  T2000 = ohana_date_to_sec ("2000/01/01");
 
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
@@ -91 +91 @@
       m = catalog[i].average[j].measureOffset;
 
-      Tmin = Tmax = (catalog[i].measure[m].t - To) / (86400*365.25);
+      Tmin = Tmax = (catalog[i].measure[m].t - T2000) / (86400*365.25);
       mode = FIT_MODE;
 
+      // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
+      Tmean = 0;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
 
@@ -125 +127 @@
         Tmin = MIN(Tmin, T[N]);
         Tmax = MAX(Tmax, T[N]);
+        Tmean += T[N];
 
         dR[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
@@ -153 +156 @@
       coords.crval1 = R[0];
       coords.crval2 = D[0];
+      Tmean /= (float) N;
       
-      /* project all of the R,D coordinates to a plane centered on this coordinate */
+      /* project all of the R,D coordinates to a plane centered on this coordinate
+         set the times to be relative to Tmean */
       for (k = 0; k < N; k++) {
         RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
         dX[k] =  dR[k];
         dY[k] =  dD[k];
+        T[k] -= Tmean;
         // fprintf (stderr, "%d %f %f %f  %f %f\n", k, T[k], R[k], D[k], X[k], Y[k]);
       }   
 
-      /* fit the model components as needed */
-      switch (mode) {
-        case FIT_AVERAGE:
+      // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
+
+      // fit the average model
+      if ((mode == FIT_AVERAGE) || (mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
           liststats (R, dR, N, &statsR);
           liststats (D, dD, N, &statsD);
 
-          fit.Ro = statsR.mean;
-          fit.dRo = 3600.0*statsR.sigma;
-
-          fit.Do = statsD.mean;
-          fit.dDo = 3600.0*statsD.sigma;
-
-          fit.chisq = 0.5*(statsR.chisq + statsD.chisq);
-          fit.Nfit = N;
-
-          fit.uR = fit.duR = 0.0;
-          fit.uD = fit.duD = 0.0;
-          fit.p  = fit.dp  = 0.0;
-
+          fitAve.Ro = statsR.mean;
+          fitAve.dRo = 3600.0*statsR.sigma;
+
+          fitAve.Do = statsD.mean;
+          fitAve.dDo = 3600.0*statsD.sigma;
+
+          fitAve.chisq = 0.5*(statsR.chisq + statsD.chisq);
+          fitAve.Nfit = N;
+
+          fitAve.uR = fitAve.duR = 0.0;
+          fitAve.uD = fitAve.duD = 0.0;
+          fitAve.p  = fitAve.dp  = 0.0;
+          catalog[i].average[j].flags |= ID_STAR_MEAS_AVE;
           Nave ++;
-          break;
-
-        case FIT_PM_ONLY:
-          FitPM (&fit, X, dX, Y, dY, T, N);
-          // fprintf (stderr, "fitted:  %f - %f : %f %f : %f %f : %f\n", Tmin, Tmax, fit.Ro, fit.Do, fit.uR, fit.uD, fit.p);
+      }
+
+      if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
+          FitPM (&fitPM, X, dX, Y, dY, T, N);
+          // fprintf (stderr, "fitted:  %f - %f : %f %f : %f %f : %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p);
           // project Ro, Do back to RA,DEC
-          XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
-          // fprintf (stderr, "project: %f %f : %f %f : %f\n", fit.Ro, fit.Do, fit.uR, fit.uD, fit.p);
+          XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords);
+          // fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p);
           // continue;
 
-          fit.p  = fit.dp  = 0.0;
-
+          fitPM.p  = fitPM.dp  = 0.0;
+          catalog[i].average[j].flags |= ID_STAR_MEAS_PAR;
           Npm ++;
-          break;
-
-        case FIT_PAR_ONLY:
+      }
+
+      if (mode == FIT_PAR_ONLY) {
+          fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
+          exit (2);
+      }
+
+# if (0)
+      if (mode == FIT_PM_AND_PAR) {
           fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
           exit (2);
@@ -204 +217 @@
             ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
           }
-          FitPar (&fit, X, dX, Y, dY, pX, pY, N);
-
-          // project Ro, Do back to RA,DEC
-          XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
-
-          fit.uR = fit.duR = 0.0;
-          fit.uD = fit.duD = 0.0;
-
+          FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N);
+          XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
+          catalog[i].average[j].flags |= ID_STAR_MEAS_PAR;
           Npar ++;
-          break;
-
-        case FIT_PM_AND_PAR:
-          fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
-          exit (2);
-
-          for (k = 0; k < N; k++) {
-            ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
-          }
-          FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N);
-          XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
-          Npar ++;
-          break;
-
-        default:
-          fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
-          exit (2);
       }   
+# endif
 
       if (0 && (j < 100)) {
@@ -257 +249 @@
       }
 
+      /* choose the result based on the chisq values */
+
+      switch (result) {
+        case FIT_AVERAGE:
+          catalog[i].average[j].flags |= ID_STAR_USE_AVE;
+          fit = fitAve;
+          break;
+        case FIT_PM_ONLY:
+          catalog[i].average[j].flags |= ID_STAR_USE_PM;
+          fit = fitPM;
+          break;
+        case FIT_PM_AND_PAR:
+          catalog[i].average[j].flags |= ID_STAR_USE_PAR;
+          fit = fitPAR;
+          break;
+      }
+
       // the measure fields must be updated before the average fields
       m = catalog[i].average[j].measureOffset;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-        // XXX why was this here?? if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
         setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
         setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
@@ -278 +286 @@
       catalog[i].average[j].dP  = fit.dp; // parallax error in arcsec
 
-      catalog[i].average[j].Xp  = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NAN_S_SHORT;
+      // Xp is supposed to be the position scatter, not the chisq : fix this:
+      // catalog[i].average[j].Xp  = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NAN_S_SHORT;
+      catalog[i].average[j].chiSqAve  = fitAve.ChiSq;
+      catalog[i].average[j].chiSqPM  = fitPM.ChiSq;
+      catalog[i].average[j].chiSqPar  = fitPar.ChiSq;
+      catalog[i].average[j].Tmean = Tmean;
     }