Index: /branches/eam_branches/relastro.20100326/src/UpdateObjects.c =================================================================== --- /branches/eam_branches/relastro.20100326/src/UpdateObjects.c (revision 27547) +++ /branches/eam_branches/relastro.20100326/src/UpdateObjects.c (revision 27548) @@ -44,5 +44,5 @@ StatType statsR, statsD; Coords coords; - PMFit fit; + PMFit fitAve, fitPM, fitPar; time_t To; off_t Nave, Npm, Npar, Nskip; @@ -64,5 +64,5 @@ // 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,7 +91,9 @@ 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,4 +127,5 @@ Tmin = MIN(Tmin, T[N]); Tmax = MAX(Tmax, T[N]); + Tmean += T[N]; dR[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA); @@ -153,49 +156,59 @@ 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,31 +217,10 @@ 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,8 +249,24 @@ } + /* 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,5 +286,10 @@ 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; }