Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 39369) +++ trunk/Ohana/src/relastro/Makefile (revision 39370) @@ -200,4 +200,26 @@ $(BIN)/fitobj.$(ARCH): $(FITOBJ) +fitobj2: $(BIN)/fitobj2.$(ARCH) + +FITOBJ2 = \ +$(SRC)/fitobj2.$(ARCH).o \ +$(SRC)/ParFactor.$(ARCH).o \ +$(SRC)/FitAstromOps.$(ARCH).o \ +$(SRC)/FitPosPMfixed.$(ARCH).o \ +$(SRC)/FitPM.$(ARCH).o \ +$(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/UpdateObjects.$(ARCH).o \ +$(SRC)/BootstrapOps.$(ARCH).o \ +$(SRC)/mkpolyterm.$(ARCH).o \ +$(SRC)/GetAstromError.$(ARCH).o \ +$(SRC)/ImageOps.$(ARCH).o \ +$(SRC)/MosaicOps.$(ARCH).o \ +$(SRC)/extra.$(ARCH).o \ +$(SRC)/plotstuff.$(ARCH).o \ +$(SRC)/fitpoly.$(ARCH).o + +$(FITOBJ2): $(INC)/relastro.h +$(BIN)/fitobj2.$(ARCH): $(FITOBJ2) + # $(SRC)/FitPM_IRLS.$(ARCH).o \ # $(SRC)/FitPMandPar_IRLS.$(ARCH).o \ Index: trunk/Ohana/src/relastro/src/FitPM.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPM.c (revision 39369) +++ trunk/Ohana/src/relastro/src/FitPM.c (revision 39370) @@ -186,5 +186,5 @@ } - FitPM_SetChisq (fit, data, points, Npoints); + if (!FitPM_SetChisq (fit, data, points, Npoints)) return FALSE; return (TRUE); } @@ -198,6 +198,8 @@ Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; + int Nfit = 0; for (i = 0; i < Npoints; i++) { if (points[i].mask) continue; // respect the mask if set + Nfit ++; wx = points[i].Wx; @@ -222,4 +224,5 @@ YT += points[i].Y*TWy; } + if (Nfit < 3) return FALSE; // X^T W X @@ -290,4 +293,9 @@ } + if (fit[0].Nfit < 3) { + fit[0].chisq = NAN; + return FALSE; + } + // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms) fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms); Index: trunk/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 39369) +++ trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 39370) @@ -189,5 +189,5 @@ } - FitPMandPar_SetChisq (fit, data, points, Npoints); + if (!FitPMandPar_SetChisq (fit, data, points, Npoints)) return FALSE; return (TRUE); } @@ -205,6 +205,8 @@ Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; + int Nfit = 0; for (i = 0; i < Npoints; i++) { if (points[i].mask) continue; // respect the mask if set + Nfit ++; wx = points[i].Wx; @@ -244,4 +246,5 @@ YT += points[i].Y*TWy; } + if (Nfit < 3) return FALSE; data->A[0][0] = Wx; @@ -321,7 +324,13 @@ } + if (fit[0].Nfit < 3) { + fit[0].chisq = NAN; + return FALSE; + } + // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms) fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms); + return TRUE; } Index: trunk/Ohana/src/relastro/src/FitPosPMfixed.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 39369) +++ trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 39370) @@ -180,5 +180,5 @@ } - FitPosPMfixed_SetChisq (fit, data, points, Npoints); + if (!FitPosPMfixed_SetChisq (fit, data, points, Npoints)) return FALSE; return (TRUE); } @@ -192,6 +192,8 @@ Wx = Wy = Tx = Ty = Xs = Ys = 0.0; + int Nfit = 0; for (i = 0; i < Npoints; i++) { if (points[i].mask) continue; // respect the mask if set + Nfit ++; wx = points[i].Wx; @@ -207,4 +209,5 @@ Ys += points[i].Y*wy; } + if (Nfit < 2) return FALSE; // X^T W X @@ -257,4 +260,9 @@ } + if (fit[0].Nfit < 2) { + fit[0].chisq = NAN; + return FALSE; + } + // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms) fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms); Index: trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 39369) +++ trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 39370) @@ -766,11 +766,11 @@ if (USE_GALAXY_MODEL) { for (i = 0; i < Nfit; i++) { - fit->uR = average->uRgal; - fit->uD = average->uDgal; + fit[i].uR = average->uRgal; + fit[i].uD = average->uDgal; } } else { for (i = 0; i < Nfit; i++) { - fit->uR = 0.0; - fit->uD = 0.0; + fit[i].uR = 0.0; + fit[i].uD = 0.0; } } Index: trunk/Ohana/src/relastro/src/fitobj.c =================================================================== --- trunk/Ohana/src/relastro/src/fitobj.c (revision 39369) +++ trunk/Ohana/src/relastro/src/fitobj.c (revision 39370) @@ -9,5 +9,5 @@ # define MJD_REF_YRS 12.4148 /* 2012/06/01,00:00:00 */ -# define MJD_J2000 51544.0 /* 2010/01/01,00:00:00 */ +# define MJD_J2000 51544.5 /* 2000/01/01,12:00:00 */ # define POS_ERROR 0.010 /* arcsec */ @@ -109,4 +109,5 @@ long A = now.tv_sec + now.tv_usec * 1000000; srand48(A); + // srand48(1); } @@ -179,10 +180,12 @@ catalog->average->uRgal = uR; catalog->average->uDgal = uD; + + // int T2000 = ohana_date_to_sec ("2000/01/01,12:00:00"); for (i = 0; i < Npoints; i++) { // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) - double Tmjd = MJD_MIN_MJD + drand48()*(MJD_MAX_MJD - MJD_MIN_MJD); + double Tyears = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS); + double Tmjd = 365.25*Tyears + MJD_J2000; catalog->measureT[i].t = ohana_mjd_to_sec (Tmjd); - double Tyears = (Tmjd - MJD_J2000) / 365.25; double pR, pD; @@ -208,8 +211,8 @@ for (i = Npoints; i < Ntotal; i++) { // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) - double Tmjd = MJD_MIN_MJD + drand48()*(MJD_MAX_MJD - MJD_MIN_MJD); + double Tyears = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS); + double Tmjd = 365.25*Tyears + MJD_J2000; catalog->measureT[i].t = ohana_mjd_to_sec (Tmjd); - double Tyears = (Tmjd - MJD_J2000) / 365.25; double pR, pD; Index: trunk/Ohana/src/relastro/src/fitobj2.c =================================================================== --- trunk/Ohana/src/relastro/src/fitobj2.c (revision 39370) +++ trunk/Ohana/src/relastro/src/fitobj2.c (revision 39370) @@ -0,0 +1,243 @@ +# include "relastro.h" + +Catalog *mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad); + +# define MJD_MIN_MJD 55197 /* 2010/01/01,00:00:00 */ +# define MJD_MAX_MJD 57023 /* 2015/01/01,00:00:00 */ +# define MJD_MIN_YRS 10.00 /* 2010/01/01,00:00:00 */ +# define MJD_MAX_YRS 14.9993 /* 2015/01/01,00:00:00 */ +# define MJD_REF_YRS 12.4148 /* 2012/06/01,00:00:00 */ + +# define MJD_J2000 51544.5 /* 2000/01/01,12:00:00 */ + +# define POS_ERROR 0.010 /* arcsec */ +# define OUT_ERROR 0.500 /* arcsec */ + +// # define N_STARS 3000 +// # define N_POINTS 100 +// # define N_OUTLIERS 5 +// # define N_BOOTSTRAP 100 + +# define dcos(THETA) cos(RAD_DEG*THETA) +# define dsin(THETA) sin(RAD_DEG*THETA) + +enum { + FIT_POS_NONE, + FIT_POS_IRLS, + FIT_POS_BOOT, + FIT_POS_NOCLIP, + FIT_PM_IRLS, + FIT_PM_BOOT, + FIT_PM_NOCLIP, + FIT_PLX_IRLS, + FIT_PLX_BOOT, + FIT_PLX_NOCLIP, +}; + +static int Nsecfilt = 0; + +int main (int argc, char **argv) { + + int N_STARS = 3000; + int N_POINTS = 100; + int N_BOOTSTRAP = 100; + float F_OUTLIERS = 0.1; + + int N; + if ((N = get_argument (argc, argv, "-Nstars"))) { + remove_argument (N, &argc, argv); + N_STARS = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-Npoints"))) { + remove_argument (N, &argc, argv); + N_POINTS = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-foutliers"))) { + remove_argument (N, &argc, argv); + F_OUTLIERS = atof (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-Nbootstrap"))) { + remove_argument (N, &argc, argv); + N_BOOTSTRAP = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + + if (argc != 2) { + fprintf (stderr, "USAGE: %s (mode)\n", argv[0]); + fprintf (stderr, " mode: pos, pos-irls, pos-boot, pm, pm-irls, pm-boot, plx, plx-irls, plx-boot\n"); + fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -Noutliers [10], -Nbootstrap [100]\n"); + exit (2); + } + + FIT_MODE = FIT_NONE; + if (!strcasecmp(argv[1], "pos")) { + FIT_MODE = FIT_AVERAGE; + } + if (!strcasecmp(argv[1], "pm")) { + FIT_MODE = FIT_PM_ONLY; + } + if (!strcasecmp(argv[1], "plx")) { + FIT_MODE = FIT_PM_AND_PAR; + } + if (FIT_MODE == FIT_NONE) { + fprintf (stderr, "USAGE: %s (mode)\n", argv[0]); + fprintf (stderr, " mode: pos, pm, plx\n"); + fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -Noutliers [10], -Nbootstrap [100]\n"); + exit (2); + } + + // some relastro globals need to be set + N_BOOTSTRAP_SAMPLES = N_BOOTSTRAP; + USE_GALAXY_MODEL = TRUE; + MaxMeanOffset = 10.0; + + // load reference table here + if (!LoadPhotcodesText ("dvo.photcodes")) { + fprintf (stderr, "failed to load photcodes\n"); + exit (2); + } + Nsecfilt = GetPhotcodeNsecfilt(); + + // init the random seed + { + struct timeval now; + gettimeofday (&now, NULL); + long A = now.tv_sec + now.tv_usec * 1000000; + srand48(A); + // srand48(1); + } + + ohana_gaussdev_init (); + + int i; + + FitStats *fitStats = FitStatsInit (N_POINTS * (1 + 2*F_OUTLIERS), N_BOOTSTRAP); + + int Nstars = N_STARS; + for (i = 0; i < Nstars; i++) { + + double Ro = 360.0*(drand48() - 0.5); + double Phi = 2.0*(drand48() - 0.5); + double Do = DEG_RAD * asin(Phi); + + double uR = 0.5*(drand48() - 0.5); + double uD = 0.5*(drand48() - 0.5); + + double plx = (FIT_MODE == FIT_PM_AND_PAR) ? 0.5*(drand48() + 0.0) : 0.0; + + int Npoints = (N_POINTS - 2)*drand48() + 2; // minimum of 2 points + int Noutliers = Npoints * F_OUTLIERS; + + // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points + Catalog *catalog = mkstar (fitStats, Ro, Do, uR, uD, plx, Npoints, Noutliers); + + if (!UpdateObjects_Chips (catalog->average, catalog->secfilt, catalog->measureT, NULL, Nsecfilt, fitStats, 0, 0)) goto escape; + + double Tmean = ohana_sec_to_mjd(catalog->average->Tmean); + double Tyear = (Tmean - MJD_J2000) / 365.25; + fprintf (stdout, "%3d %3d %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f | %8.6f %8.6f %8.6f %8.6f %8.6f %d | %f %f %f\n", + Npoints, Noutliers, Ro, Do, uR, uD, plx, + catalog->average->R, catalog->average->D, catalog->average->uR, catalog->average->uD, catalog->average->P, Tyear, + catalog->average->dR, catalog->average->dD, catalog->average->duR, catalog->average->duD, catalog->average->dP, catalog->average->Npos, + catalog->average->ChiSqAve, catalog->average->ChiSqPM, catalog->average->ChiSqPar); + + continue; + + escape: + fprintf (stdout, "%3d %3d %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f | %8.6f %8.6f %8.6f %8.6f %8.6f %d | %f\n", + Npoints, Noutliers, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, 0.0, NAN, NAN, NAN, NAN, NAN, 0, 0.0); + } + // return exit_status(); + ohana_gaussdev_free(); + + exit (0); +} + +Catalog *mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad) { + + int i; + + int Ntotal = Npoints + Nbad; + + ALLOCATE_PTR (catalog, Catalog, 1); + ALLOCATE (catalog->average, Average, 1); + ALLOCATE (catalog->secfilt, SecFilt, Nsecfilt); + ALLOCATE (catalog->measureT, MeasureTiny, Npoints + Nbad); + + dvo_average_init (catalog->average); + for (i = 0; i < Ntotal; i++) { + dvo_measureT_init (&catalog->measureT[i]); + } + + catalog->average->objID = 1; + catalog->average->catID = 1; + + catalog->average->R = Ro; + catalog->average->D = Do; + catalog->average->uR = uR; + catalog->average->uD = uD; + catalog->average->uRgal = uR; + catalog->average->uDgal = uD; + + // int T2000 = ohana_date_to_sec ("2000/01/01,12:00:00"); + for (i = 0; i < Npoints; i++) { + // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) + double Tyears = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS); + double Tmjd = 365.25*Tyears + MJD_J2000; + + catalog->measureT[i].t = ohana_mjd_to_sec (Tmjd); + + double pR, pD; + ParFactor (&pR, &pD, Ro, Do, Tyears); + + double dR = ohana_gaussdev_rnd(0.0, POS_ERROR); + double dD = ohana_gaussdev_rnd(0.0, POS_ERROR); + + catalog->measureT[i].R = Ro + (dR + plx*pR + uR*(Tyears - MJD_REF_YRS))/3600/dcos(Do); + catalog->measureT[i].D = Do + (dD + plx*pD + uD*(Tyears - MJD_REF_YRS))/3600; + + catalog->measureT[i].dXccd = ToShortPixels(POS_ERROR); + catalog->measureT[i].dYccd = ToShortPixels(POS_ERROR); + + catalog->measureT[i].photcode = 10233; + + // nominal values to avoid dropping in MeasFilterTest + catalog->measureT[i].M = 20.0; + catalog->measureT[i].dM = 0.02; + catalog->measureT[i].dt = 2.5*log10(30.0); + } + + for (i = Npoints; i < Ntotal; i++) { + // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) + double Tyears = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS); + double Tmjd = 365.25*Tyears + MJD_J2000; + + catalog->measureT[i].t = ohana_mjd_to_sec (Tmjd); + + double pR, pD; + ParFactor (&pR, &pD, Ro, Do, Tyears); + + double dR = OUT_ERROR*(drand48() - 0.5); + double dD = OUT_ERROR*(drand48() - 0.5); + + catalog->measureT[i].R = Ro + (dR + plx*pR + uR*(Tyears - MJD_REF_YRS))/3600/dcos(Do); + catalog->measureT[i].D = Do + (dD + plx*pD + uD*(Tyears - MJD_REF_YRS))/3600; + + catalog->measureT[i].dXccd = ToShortPixels(POS_ERROR); + catalog->measureT[i].dYccd = ToShortPixels(POS_ERROR); + + catalog->measureT[i].photcode = 10233; + + // nominal values to avoid dropping in MeasFilterTest + catalog->measureT[i].M = 20.0; + catalog->measureT[i].dM = 0.02; + catalog->measureT[i].dt = 2.5*log10(30.0); + } + catalog->average->Nmeasure = Ntotal; + + return catalog; +} Index: trunk/Ohana/src/relastro/src/fitpm.c =================================================================== --- trunk/Ohana/src/relastro/src/fitpm.c (revision 39369) +++ trunk/Ohana/src/relastro/src/fitpm.c (revision 39370) @@ -33,4 +33,8 @@ }; +float ERROR_BIAS = 1.0; +float RA_BIAS = 0.0; +float DEC_BIAS = 0.0; + int main (int argc, char **argv) { @@ -59,4 +63,19 @@ remove_argument (N, &argc, argv); N_BOOTSTRAP = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-error-bias"))) { + remove_argument (N, &argc, argv); + ERROR_BIAS = atof (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-ra-bias"))) { + remove_argument (N, &argc, argv); + RA_BIAS = atof (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-dec-bias"))) { + remove_argument (N, &argc, argv); + DEC_BIAS = atof (argv[N]); remove_argument (N, &argc, argv); } @@ -108,4 +127,5 @@ long A = now.tv_sec + now.tv_usec * 1000000; srand48(A); + srand48(1); } @@ -184,5 +204,5 @@ fitStats->fit[k].uR = uR; fitStats->fit[k].uD = uD; - if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue; + if (!FitPosPMfixed_Basic (&fitStats->fit[fitStats->Nfit], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue; fitStats->Nfit ++; } @@ -204,5 +224,5 @@ for (k = 0; k < fitStats->NfitAlloc; k++) { BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); - if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, Nnomask)) continue; + if (!FitPM_Basic (&fitStats->fit[fitStats->Nfit], fitStats->fitdataPM, fitStats->sample, Nnomask)) continue; fitStats->Nfit ++; } @@ -226,5 +246,5 @@ for (k = 0; k < fitStats->NfitAlloc; k++) { BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); - if (!FitPMandPar_Basic (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask)) continue; + if (!FitPMandPar_Basic (&fitStats->fit[fitStats->Nfit], fitStats->fitdataPar, fitStats->sample, Nnomask)) continue; fitStats->Nfit ++; } @@ -282,6 +302,6 @@ points[i].D = Do + (dD + plx*pD + uD*(points[i].T - MJD_REF_YRS))/3600; - points[i].dX = POS_ERROR; - points[i].dY = POS_ERROR; + points[i].dX = POS_ERROR*ERROR_BIAS; + points[i].dY = POS_ERROR*ERROR_BIAS; } @@ -295,12 +315,12 @@ ParFactor (&pR, &pD, Ro, Do, points[i].T); - double dR = OUT_ERROR*(drand48() - 0.5); - double dD = OUT_ERROR*(drand48() - 0.5); + double dR = OUT_ERROR*(drand48() - 0.5) + RA_BIAS; + double dD = OUT_ERROR*(drand48() - 0.5) + DEC_BIAS; points[i].R = Ro + (dR + plx*pR + uR*(points[i].T - MJD_REF_YRS))/3600/dcos(Do); points[i].D = Do + (dD + plx*pD + uD*(points[i].T - MJD_REF_YRS))/3600; - points[i].dX = POS_ERROR; - points[i].dY = POS_ERROR; + points[i].dX = POS_ERROR*ERROR_BIAS; + points[i].dY = POS_ERROR*ERROR_BIAS; } fitStats->Npoints = Ntotal; Index: trunk/Ohana/src/relastro/test/mana.sh =================================================================== --- trunk/Ohana/src/relastro/test/mana.sh (revision 39369) +++ trunk/Ohana/src/relastro/test/mana.sh (revision 39370) @@ -30,4 +30,23 @@ end +macro testobj + if ($0 != 1) + echo "USAGE: testobj" + break + end + + foreach mode pos pm plx + exec ../bin/fitobj.lin64 $mode > test.$mode.dat + if ("$mode" == "plx") + echo "is plx" + reload test.$mode.dat 1 + else + reload test.$mode.dat 0 + end + + fprintf "%5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq + end +end + macro testnboot if ($0 != 2) @@ -38,5 +57,5 @@ foreach nboot 20 30 100 300 1000 exec ../bin/fitpm.lin64 -Nstars 3000 -Noutliers 20 -Npoints 100 -Nbootstrap $nboot $1 > test.$1.dat - reload test.$1.dat + reload test.$1.dat 0 fprintf "%4d : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $nboot $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq @@ -50,7 +69,9 @@ end - foreach outliers 0 3 10 30 100 + # foreach outliers 0 3 10 30 100 + + foreach outliers 10 exec ../bin/fitpm.lin64 -Nstars 3000 -Noutliers $outliers -Npoints 100 -Nbootstrap 100 $1 > test.$1.dat - reload test.$1.dat + reload test.$1.dat 0 fprintf "%4d : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $outliers $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq @@ -59,11 +80,16 @@ macro reload - if ($0 != 2) - echo "USAGE: reload (file)" + if ($0 != 3) + echo "USAGE: reload (file) (isObjPLX)" break end data $1 - read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22 + if ($2) + read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26 + # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24 + else + read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22 + end foreach f R D uD uR P @@ -98,4 +124,49 @@ end +macro teststats + if ($0 != 2) + echo "USAGE: teststats (file)" + break + end + + data $1 + # read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26 + # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24 + read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22 + + foreach f R D uD uR P + set d$f = $f\o - $f\x + end + + set Rf = Rx - uRx*(Tx - $Tref)/3600.0/dcos(Dx) + set Df = Dx - uDx*(Tx - $Tref)/3600.0 + + set dRf = 3600*(Ro - Rf)*dcos(Dx) + set dDf = 3600*(Do - Df) + + # echo "fit scatter" + vstat -q dRf; $value_dRf = $SIGMA*1000 + vstat -q dDf; $value_dDf = $SIGMA*1000 + vstat -q duR; $value_duR = $SIGMA*1000 + vstat -q duD; $value_duD = $SIGMA*1000 + vstat -q dP ; $value_dP = $SIGMA*1000 + + # echo "fit errors" + vstat -q dRo; $value_dRo = $MEAN*1000; $sigma_dRo = $SIGMA*1000 + vstat -q dDo; $value_dDo = $MEAN*1000; $sigma_dDo = $SIGMA*1000 + vstat -q duRo; $value_duRo = $MEAN*1000; $sigma_duRo = $SIGMA*1000 + vstat -q duDo; $value_duDo = $MEAN*1000; $sigma_duDo = $SIGMA*1000 + vstat -q dPo ; $value_dPo = $MEAN*1000; $sigma_dPo = $SIGMA*1000 + + # echo "Nfit" + vstat -q Nfit; $mean_Nfit = $MEAN + + # echo chisq + vstat -q chisq; $mean_chisq = $MEAN + + fprintf "%5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq + +end + # determine the parallax factor (time) (ra) (dec) (parfR) (parfD) are vectors, macro getParFactor_vectors