Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 39611) +++ trunk/Ohana/src/relastro/Makefile (revision 39612) @@ -1,3 +1,3 @@ -default: relastro relastro_client fitobj fitobj2 fitpm +default: relastro relastro_client fitobj fitobj2 fitobj3 fitpm help: @echo "make options: relastro (default)" @@ -238,4 +238,26 @@ $(BIN)/fitobj2.$(ARCH): $(FITOBJ2) +fitobj3: $(BIN)/fitobj3.$(ARCH) + +FITOBJ3 = \ +$(SRC)/fitobj3.$(ARCH).o \ +$(SRC)/ParFactor.$(ARCH).o \ +$(SRC)/FitAstromOps.$(ARCH).o \ +$(SRC)/FitPosPMfixed.$(ARCH).o \ +$(SRC)/FitPM.$(ARCH).o \ +$(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitPosStack.$(ARCH).o \ +$(SRC)/UpdateObjects.$(ARCH).o \ +$(SRC)/BootstrapOps.$(ARCH).o \ +$(SRC)/mkpolyterm.$(ARCH).o \ +$(SRC)/GetAstromError.$(ARCH).o \ +$(SRC)/MeasFilterTest.$(ARCH).o \ +$(SRC)/extra.$(ARCH).o \ +$(SRC)/plotstuff.$(ARCH).o \ +$(SRC)/fitpoly.$(ARCH).o + +$(FITOBJ3): $(INC)/relastro.h +$(BIN)/fitobj3.$(ARCH): $(FITOBJ3) + fitstk: $(BIN)/fitstk.$(ARCH) Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 39611) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 39612) @@ -5,4 +5,6 @@ # include # include + +# define OLD_METHOD 0 // choose off_t or int depending on full-scale relphot analysis resources Index: trunk/Ohana/src/relastro/src/FitPM.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPM.c (revision 39611) +++ trunk/Ohana/src/relastro/src/FitPM.c (revision 39612) @@ -1,3 +1,4 @@ # include "relastro.h" + int FitPM_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); int FitPM_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); @@ -16,6 +17,10 @@ // Convert the measurement errors into initial weights. for (i = 0; i < Npoints; i++) { - points[i].Wx = 1.0 / SQ(points[i].dX); - points[i].Wy = 1.0 / SQ(points[i].dY); + points[i].Wx = 1.0; + points[i].Wy = 1.0; + points[i].Qx = 1.0 / SQ(points[i].dX); + points[i].Qy = 1.0 / SQ(points[i].dY); + points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0 + points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0 } @@ -33,6 +38,10 @@ // Convert the measurement errors into initial weights. for (i = 0; i < Npoints; i++) { - points[i].Wx = 1.0 / SQ(points[i].dX); - points[i].Wy = 1.0 / SQ(points[i].dY); + points[i].Wx = 1.0; + points[i].Wy = 1.0; + points[i].Qx = 1.0 / SQ(points[i].dX); + points[i].Qy = 1.0 / SQ(points[i].dY); + points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0 + points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0 } @@ -76,4 +85,6 @@ points[i].Wx = weight_cauchy(points[i].rx / points[i].dX); points[i].Wy = weight_cauchy(points[i].ry / points[i].dY); + points[i].qx = points[i].Wx * points[i].Qx; + points[i].qy = points[i].Wy * points[i].Qy; } @@ -201,6 +212,6 @@ Nfit ++; - wx = points[i].Wx; - wy = points[i].Wy; + wx = points[i].qx; + wy = points[i].qy; Wx += wx; @@ -298,5 +309,4 @@ // 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/FitPMandPar.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 39611) +++ trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 39612) @@ -1,4 +1,3 @@ # include "relastro.h" -# define OLD_METHOD 0 int FitPMandPar_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); @@ -18,8 +17,4 @@ // Convert the measurement errors into initial weights. for (i = 0; i < Npoints; i++) { -# if (OLD_METHOD) - points[i].Wx = 1.0 / SQ(points[i].dX); - points[i].Wy = 1.0 / SQ(points[i].dY); -# else points[i].Wx = 1.0; points[i].Wy = 1.0; @@ -28,9 +23,7 @@ points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0 points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0 -# endif } if (!FitPMandPar_MinChisq (fit, data, points, Npoints)) return FALSE; - if (!FitPMandPar_SetChisq (fit, data, points, Npoints)) return FALSE; return TRUE; @@ -45,15 +38,10 @@ // Convert the measurement errors into initial weights. for (i = 0; i < Npoints; i++) { -# if (OLD_METHOD) - points[i].Wx = 1.0 / SQ(points[i].dX); - points[i].Wy = 1.0 / SQ(points[i].dY); - points[i].qx = points[i].Wx; // Wx, Wy start out at 1.0 - points[i].qy = points[i].Wy; // Wx, Wy start out at 1.0 -# else + points[i].Wx = 1.0; + points[i].Wy = 1.0; points[i].Qx = 1.0 / SQ(points[i].dX); points[i].Qy = 1.0 / SQ(points[i].dY); points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0 points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0 -# endif } @@ -231,11 +219,6 @@ Nfit ++; -# if (OLD_METHOD) - wx = points[i].Wx; - wy = points[i].Wy; -# else wx = points[i].qx; wy = points[i].qy; -# endif Wx += wx; Index: trunk/Ohana/src/relastro/src/FitPosPMfixed.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 39611) +++ trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 39612) @@ -1,3 +1,4 @@ # include "relastro.h" + int FitPosPMfixed_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); int FitPosPMfixed_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); @@ -45,6 +46,10 @@ // Convert the measurement errors into initial weights. for (i = 0; i < Npoints; i++) { - points[i].Wx = 1.0 / SQ(points[i].dX); - points[i].Wy = 1.0 / SQ(points[i].dY); + points[i].Wx = 1.0; + points[i].Wy = 1.0; + points[i].Qx = 1.0 / SQ(points[i].dX); + points[i].Qy = 1.0 / SQ(points[i].dY); + points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0 + points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0 } @@ -69,6 +74,10 @@ // Convert the measurement errors into initial weights. for (i = 0; i < Npoints; i++) { - points[i].Wx = 1.0 / SQ(points[i].dX); - points[i].Wy = 1.0 / SQ(points[i].dY); + points[i].Wx = 1.0; + points[i].Wy = 1.0; + points[i].Qx = 1.0 / SQ(points[i].dX); + points[i].Qy = 1.0 / SQ(points[i].dY); + points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0 + points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0 } @@ -112,4 +121,6 @@ points[i].Wx = weight_cauchy(points[i].rx / points[i].dX); points[i].Wy = weight_cauchy(points[i].ry / points[i].dY); + points[i].qx = points[i].Wx * points[i].Qx; + points[i].qy = points[i].Wy * points[i].Qy; } @@ -231,6 +242,6 @@ Nfit ++; - wx = points[i].Wx; - wy = points[i].Wy; + wx = points[i].qx; + wy = points[i].qy; Wx += wx; Index: trunk/Ohana/src/relastro/src/fitobj.c =================================================================== --- trunk/Ohana/src/relastro/src/fitobj.c (revision 39611) +++ trunk/Ohana/src/relastro/src/fitobj.c (revision 39612) @@ -69,5 +69,6 @@ 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, " mode: pos, pm, plx\n"); + fprintf (stderr, " test astrometry fitting using UpdateObjects_Chip code\n"); fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -Noutliers [10], -Nbootstrap [100]\n"); exit (2); @@ -98,5 +99,5 @@ // load reference table here if (!LoadPhotcodesText ("dvo.photcodes")) { - fprintf (stderr, "failed to load photcodes\n"); + fprintf (stderr, "failed to load photcode file dvo.photcodes (copy from ipp/ippconfig)\n"); exit (2); } @@ -121,4 +122,5 @@ for (i = 0; i < Nstars; i++) { + // generate random stars uniformly distributed on the sphere (ra = 0,36) double Ro = 360.0*drand48(); double Phi = 2.0*(drand48() - 0.5); Index: trunk/Ohana/src/relastro/src/fitobj2.c =================================================================== --- trunk/Ohana/src/relastro/src/fitobj2.c (revision 39611) +++ trunk/Ohana/src/relastro/src/fitobj2.c (revision 39612) @@ -48,4 +48,5 @@ int N_STARS = 3000; int N_POINTS = 100; + int N_MIN_POINTS = 10; int N_BOOTSTRAP = 100; float F_OUTLIERS = 0.1; @@ -62,4 +63,9 @@ remove_argument (N, &argc, argv); } + if ((N = get_argument (argc, argv, "-Nminpoints"))) { + remove_argument (N, &argc, argv); + N_MIN_POINTS = atoi (argv[N]); + remove_argument (N, &argc, argv); + } if ((N = get_argument (argc, argv, "-foutliers"))) { remove_argument (N, &argc, argv); @@ -76,6 +82,7 @@ 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"); + fprintf (stderr, " mode: pos, pm, plx\n"); + fprintf (stderr, " test astrometry fitting using UpdateObjects_Chip code (random number of measurements)\n"); + fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -foutliers [0.1], -Nbootstrap [100]\n"); exit (2); } @@ -94,5 +101,5 @@ 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"); + fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -foutliers [0.1], -Nbootstrap [100]\n"); exit (2); } @@ -105,5 +112,5 @@ // load reference table here if (!LoadPhotcodesText ("dvo.photcodes")) { - fprintf (stderr, "failed to load photcodes\n"); + fprintf (stderr, "failed to load photcode file dvo.photcodes (copy from ipp/ippconfig)\n"); exit (2); } @@ -128,5 +135,6 @@ for (i = 0; i < Nstars; i++) { - double Ro = 360.0*(drand48_cnt() - 0.5); + // generate random stars uniformly distributed on the sphere (ra = 0,36) + double Ro = 360.0*drand48_cnt(); double Phi = 2.0*(drand48_cnt() - 0.5); double Do = DEG_RAD * asin(Phi); @@ -138,5 +146,5 @@ if (FIT_MODE != FIT_PM_AND_PAR) plx = 0.0; - int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 points + int Npoints = (N_POINTS - N_MIN_POINTS)*drand48_cnt() + N_MIN_POINTS; // randomly choose N_MIN_POINTS to N_POINTS points int Noutliers_max = Npoints * F_OUTLIERS; int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points Index: trunk/Ohana/src/relastro/src/fitobj3.c =================================================================== --- trunk/Ohana/src/relastro/src/fitobj3.c (revision 39612) +++ trunk/Ohana/src/relastro/src/fitobj3.c (revision 39612) @@ -0,0 +1,293 @@ +# include "relastro.h" + +static int Nrand = 0; +double drand48_cnt () { + double value = drand48(); + Nrand ++; + return value; +} + +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_MIN_POINTS = 10; + 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, "-Nminpoints"))) { + remove_argument (N, &argc, argv); + N_MIN_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, pm, plx\n"); + fprintf (stderr, " test astrometry fitting using UpdateObjects_Chip code (random number of measurements)\n"); + fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -foutliers [0.1], -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], -foutliers [0.1], -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 photcode file dvo.photcodes (copy from ipp/ippconfig)\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++) { + + // generate random stars uniformly distributed on the sphere (ra = 0,36) + double Ro = 360.0*drand48_cnt(); + double Phi = 2.0*(drand48_cnt() - 0.5); + double Do = DEG_RAD * asin(Phi); + + double uR = 0.5*(drand48_cnt() - 0.5); + double uD = 0.5*(drand48_cnt() - 0.5); + + double plx = 0.5*(drand48_cnt() + 0.0); + if (FIT_MODE != FIT_PM_AND_PAR) plx = 0.0; + + int Npoints = (N_POINTS - N_MIN_POINTS)*drand48_cnt() + N_MIN_POINTS; // randomly choose N_MIN_POINTS to N_POINTS points + int Noutliers_max = Npoints * F_OUTLIERS; + int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points + + if (Npoints + Noutliers == 0) continue; + + // 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 | %6.2f %6.2f %6.2f 0x%08x\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, catalog->average->flags); + if (i > 0) { + // fprintf (stderr, ".\n"); + } + + 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) { + OHANA_UNUSED_PARAM(fitStats); + + 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_cnt()*(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 posError = POS_ERROR + POS_ERROR*3*drand48_cnt(); + + double dR = ohana_gaussdev_rnd(0.0, posError); + double dD = ohana_gaussdev_rnd(0.0, posError); + + 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(posError); + catalog->measureT[i].dYccd = ToShortPixels(posError); + + 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_cnt()*(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_cnt() - 0.5); + double dD = OUT_ERROR*(drand48_cnt() - 0.5); + + double posError = POS_ERROR + POS_ERROR*3*drand48_cnt(); + + 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(posError); + catalog->measureT[i].dYccd = ToShortPixels(posError); + + 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); + } + + if (Ntotal > 0) { + catalog->average->R = catalog->measureT[0].R; + catalog->average->D = catalog->measureT[0].D; + catalog->average->uRgal = uR; + catalog->average->uDgal = uD; + } + + catalog->average->Nmeasure = Ntotal; + + return catalog; +} + +// these dummy functions are used to avoid including ImageOps, etc +int areImagesMatched () { + return FALSE; +} +float getColorBlue (off_t meas, int cat) { + OHANA_UNUSED_PARAM(meas); + OHANA_UNUSED_PARAM(cat); + return (NAN); +} +float getColorRed (off_t meas, int cat) { + OHANA_UNUSED_PARAM(meas); + OHANA_UNUSED_PARAM(cat); + return (NAN); +} Index: trunk/Ohana/src/relastro/test/mana.sh =================================================================== --- trunk/Ohana/src/relastro/test/mana.sh (revision 39611) +++ trunk/Ohana/src/relastro/test/mana.sh (revision 39612) @@ -40,4 +40,20 @@ end +macro testobj3 + if ($0 != 3) + echo "USAGE: testobj3 (Nstars) (foutliers)" + break + end + + $Nmin:pos = 99 + $Nmin:pm = 99 + $Nmin:plx = 99 + + foreach mode pos pm plx + exec ../bin/fitobj3.lin64 $mode -Nstars $1 -foutliers $2 -Nminpoints $Nmin:$mode > test.$mode.dat + checkstats test.$mode.dat $mode + end +end + macro testobj2 if ($0 != 3) @@ -46,6 +62,12 @@ end - exec ../bin/fitobj2.lin64 plx -Nstars $1 -foutliers $2 > test.plx.dat - checkstats test.plx.dat + $Nmin:pos = 1 + $Nmin:pm = 5 + $Nmin:plx = 7 + + foreach mode pos pm plx + exec ../bin/fitobj2.lin64 $mode -Nstars $1 -foutliers $2 -Nminpoints $Nmin:$mode > test.$mode.dat + checkstats test.$mode.dat $mode + end end @@ -58,11 +80,5 @@ 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 - + checkstats.obj test.$mode.dat $mode 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 @@ -108,5 +124,5 @@ data $1 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 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 @@ -145,14 +161,20 @@ end -macro teststats - if ($0 != 2) - echo "USAGE: teststats (file)" +macro checkstats.obj + if ($0 != 3) + echo "USAGE: checkstats.obj (file) (mode)" 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 + if ("$2" == "pos") + 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 + if ("$2" == "pm") + 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 23 + end + if ("$2" == "plx") + 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 + end foreach f R D uD uR P @@ -185,47 +207,27 @@ # 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 - -macro checkstats - if ($0 != 2) - echo "USAGE: checkstats (file)" - break - end - - teststats $1 - - delete dRbin oRbin - for i 0 100 - subset tmp = dRf if (Npts == $i) - vstat tmp -q - concat $SIGMA dRbin - concat $MEAN oRbin - end - - create n 0 dRbin[] - set dRmod = 0.010 / sqrt (n) - lim n -0.001 0.015; clear; box; plot n dRbin - plot -c red -pt 7 n dRmod -end - -macro load.stack - if ($0 != 2) - echo "USAGE: load.stack (file)" +end + +macro teststats + if ($0 != 3) + echo "USAGE: teststats (file) (mode)" break end data $1 - read Npts 1 Nout 2 Nstk 3 Ro 4 Do 5 uRo 6 uDo 7 Po 8 Rx 10 Dx 11 uRx 12 uDx 3 Px 14 Tx 16 dRo 18 dDo 19 duRo 20 duDo 21 dPo 22 Nfit 23 Rstk 25 Dstk 26 dRstk 27 dDstk 28 chisq 32 - # 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 uD uR P + + if ("$2" == "pos") + 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 24 + end + if ("$2" == "pm") + 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 25 + end + if ("$2" == "plx") + 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 + end + + foreach f R D uD uR P set d$f = $f\o - $f\x end - set sdR = 3600*(Ro - Rstk)*dcos(Do) - set sdD = 3600*(Do - Dstk) set Rf = Rx - uRx*(Tx - $Tref)/3600.0/dcos(Dx) @@ -242,7 +244,4 @@ vstat -q dP ; $value_dP = $SIGMA*1000 - vstat -q sdR; $value_sdR = $SIGMA*1000 - vstat -q sdD; $value_sdD = $SIGMA*1000 - # echo "fit errors" vstat -q dRo; $value_dRo = $MEAN*1000; $sigma_dRo = $SIGMA*1000 @@ -259,6 +258,78 @@ 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 + +macro checkstats + if ($0 != 3) + echo "USAGE: checkstats (file) (mode)" + break + end + + teststats $1 $2 - delete dRbin oRbin + delete -q dRbin oRbin + for i 0 100 + subset tmp = dRf if (Npts == $i) + vstat tmp -q + concat $SIGMA dRbin + concat $MEAN oRbin + end + + create n 0 dRbin[] + set dRmod = 0.010 / sqrt (n) + lim n -0.001 0.015; clear; box; plot n dRbin + plot -c red -pt 7 n dRmod +end + +macro load.stack + if ($0 != 2) + echo "USAGE: load.stack (file)" + break + end + + data $1 + read Npts 1 Nout 2 Nstk 3 Ro 4 Do 5 uRo 6 uDo 7 Po 8 Rx 10 Dx 11 uRx 12 uDx 3 Px 14 Tx 16 dRo 18 dDo 19 duRo 20 duDo 21 dPo 22 Nfit 23 Rstk 25 Dstk 26 dRstk 27 dDstk 28 chisq 32 + # 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 uD uR P + set d$f = $f\o - $f\x + end + set sdR = 3600*(Ro - Rstk)*dcos(Do) + set sdD = 3600*(Do - Dstk) + + 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 + + vstat -q sdR; $value_sdR = $SIGMA*1000 + vstat -q sdD; $value_sdD = $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 + + delete -q dRbin oRbin for i 0 100 subset tmp = dRf if (Npts == $i)