Index: /trunk/Ohana/src/relastro/Makefile =================================================================== --- /trunk/Ohana/src/relastro/Makefile (revision 39374) +++ /trunk/Ohana/src/relastro/Makefile (revision 39375) @@ -1,3 +1,3 @@ -default: relastro relastro_client +default: relastro relastro_client fitobj fitobj2 fitpm help: @echo "make options: relastro (default)" @@ -33,4 +33,5 @@ $(SRC)/FitPMandPar.$(ARCH).o \ $(SRC)/FitPosPMfixed.$(ARCH).o \ +$(SRC)/FitPosStack.$(ARCH).o \ $(SRC)/BootstrapOps.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ @@ -108,4 +109,5 @@ $(SRC)/FitPMandPar.$(ARCH).o \ $(SRC)/FitPosPMfixed.$(ARCH).o \ +$(SRC)/FitPosStack.$(ARCH).o \ $(SRC)/BootstrapOps.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ @@ -171,4 +173,5 @@ $(SRC)/FitPM.$(ARCH).o \ $(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitPosStack.$(ARCH).o \ $(SRC)/BootstrapOps.$(ARCH).o \ $(SRC)/mkpolyterm.$(ARCH).o \ @@ -187,4 +190,5 @@ $(SRC)/FitPM.$(ARCH).o \ $(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitPosStack.$(ARCH).o \ $(SRC)/UpdateObjects.$(ARCH).o \ $(SRC)/BootstrapOps.$(ARCH).o \ @@ -209,4 +213,5 @@ $(SRC)/FitPM.$(ARCH).o \ $(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitPosStack.$(ARCH).o \ $(SRC)/UpdateObjects.$(ARCH).o \ $(SRC)/BootstrapOps.$(ARCH).o \ @@ -221,4 +226,27 @@ $(FITOBJ2): $(INC)/relastro.h $(BIN)/fitobj2.$(ARCH): $(FITOBJ2) + +fitstk: $(BIN)/fitstk.$(ARCH) + +FITSTK = \ +$(SRC)/fitstk.$(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)/ImageOps.$(ARCH).o \ +$(SRC)/MosaicOps.$(ARCH).o \ +$(SRC)/extra.$(ARCH).o \ +$(SRC)/plotstuff.$(ARCH).o \ +$(SRC)/fitpoly.$(ARCH).o + +$(FITSTK): $(INC)/relastro.h +$(BIN)/fitstk.$(ARCH): $(FITSTK) # $(SRC)/FitPM_IRLS.$(ARCH).o \ Index: /trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- /trunk/Ohana/src/relastro/include/relastro.h (revision 39374) +++ /trunk/Ohana/src/relastro/include/relastro.h (revision 39375) @@ -173,4 +173,9 @@ int Nfit; int NfitAlloc; + + double * Xstack; + double *dXstack; + double * Ystack; + double *dYstack; FitAstromPoint *points; @@ -755,4 +760,5 @@ int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average); void AstromErrorSetLoop (int Nloop, int isImageMode); +void FitPointsClearMasks (FitAstromPoint *points, int Npoints); int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); @@ -760,4 +766,5 @@ int FitPosPMfixed_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); int FitPosPMfixed_Single (FitAstromResult *fit, FitAstromPoint *points, int Npoints); +int FitPosStack (FitAstromResult *fit, FitStats *fitstats); double MedianAbsDeviation(FitAstromPoint *points, int Npoints); Index: /trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 39374) +++ /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 39375) @@ -84,4 +84,11 @@ SetZeroPoint (25.0); + if (USE_GALAXY_MODEL) { + if (!InitGalaxyModel ("FEAST-HIPPARCOS")) { + fprintf (stderr, "failed to init galaxy model\n"); + exit (2); + } + } + free (config); free (file); Index: /trunk/Ohana/src/relastro/src/FitAstromOps.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitAstromOps.c (revision 39374) +++ /trunk/Ohana/src/relastro/src/FitAstromOps.c (revision 39375) @@ -85,4 +85,8 @@ // containers to hold the measurements for a given star + fitStats-> Xstack = NULL; // pre-allocated array for median + fitStats->dXstack = NULL; // pre-allocated array for median + fitStats-> Ystack = NULL; // pre-allocated array for median + fitStats->dYstack = NULL; // pre-allocated array for median fitStats->points = NULL; fitStats->sample = NULL; @@ -94,4 +98,8 @@ ALLOCATE (fitStats->sample, FitAstromPoint, Nmax); ALLOCATE (fitStats->nomask, FitAstromPoint, Nmax); + ALLOCATE (fitStats-> Xstack, double, Nmax); + ALLOCATE (fitStats->dXstack, double, Nmax); + ALLOCATE (fitStats-> Ystack, double, Nmax); + ALLOCATE (fitStats->dYstack, double, Nmax); } @@ -143,4 +151,9 @@ FREE (fitStats->nomask); + FREE (fitStats-> Xstack); + FREE (fitStats->dXstack); + FREE (fitStats-> Ystack); + FREE (fitStats->dYstack); + FitAstromDataFree (fitStats->fitdataPos); FitAstromDataFree (fitStats->fitdataPM); @@ -168,4 +181,13 @@ return fit; +} + +void FitPointsClearMasks (FitAstromPoint *points, int Npoints) { + + int i; + for (i = 0; i < Npoints; i++) { + points[i].mask = FALSE; + } + return; } Index: /trunk/Ohana/src/relastro/src/FitPosPMfixed.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 39374) +++ /trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 39375) @@ -26,4 +26,5 @@ fit->duD = 0.0; fit->dp = 0.0; + fit->Nfit = 1; found = i; Index: /trunk/Ohana/src/relastro/src/ICRF.c =================================================================== --- /trunk/Ohana/src/relastro/src/ICRF.c (revision 39374) +++ /trunk/Ohana/src/relastro/src/ICRF.c (revision 39375) @@ -63,5 +63,5 @@ for (i = 0; i < catalog[c].Naverage; i++) { - if (!(catalog[c].average[i].flags & ID_ICRF_QSO)) continue; + if (!(catalog[c].average[i].flags & ID_OBJ_ICRF_QSO)) continue; // only save a single value Index: /trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 39374) +++ /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 39375) @@ -116,5 +116,5 @@ // if we fail to fit the astrometry for some reason, we need to set/reset these - average[0].flags |= ID_STAR_NO_ASTROM; + average[0].flags |= ID_OBJ_NO_MEAN_ASTROM; average[0].ChiSqAve = NAN; average[0].ChiSqPM = NAN; @@ -142,5 +142,5 @@ average[0].dR = average[0].dRstk; average[0].dD = average[0].dDstk; - average[0].flags |= ID_STACK_ASTROM; + average[0].flags |= ID_OBJ_STACK_FOR_MEAN; } return FALSE; @@ -207,5 +207,5 @@ if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); - average[0].flags |= ID_STAR_FIT_PAR; + average[0].flags |= ID_OBJ_FIT_PAR; fitStats->Npar ++; @@ -225,5 +225,5 @@ mode = FIT_PM_ONLY; } else { - average[0].flags |= ID_STAR_USE_PAR; + average[0].flags |= ID_OBJ_USE_PAR; } } @@ -242,8 +242,8 @@ } - if (average[0].flags & ID_STAR_USE_PAR) { + if (average[0].flags & ID_OBJ_USE_PAR) { // get the best pm fit and chisq given the set of points (mask is respected) if (!FitPM_Basic (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) { - average[0].flags |= ID_STAR_BAD_PM; + average[0].flags |= ID_OBJ_BAD_PM; goto justPosition; } @@ -277,5 +277,5 @@ if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); - average[0].flags |= ID_STAR_FIT_PM; + average[0].flags |= ID_OBJ_FIT_PM; fitStats->Npm ++; @@ -290,8 +290,8 @@ if (!valid) { mode = FIT_AVERAGE; - average[0].flags |= ID_STAR_BAD_PM; + average[0].flags |= ID_OBJ_BAD_PM; } else { - if (!(average[0].flags & ID_STAR_USE_PAR)) { - average[0].flags |= ID_STAR_USE_PM; + if (!(average[0].flags & ID_OBJ_USE_PAR)) { + average[0].flags |= ID_OBJ_USE_PM; } } @@ -305,6 +305,7 @@ FitAstromResultSetPM (&fitPos, 1, average); // fprintf (stderr, "fit 1: %f %f : %f %f\n", fitPos.Ro, fitPos.Do, fitPos.uR, fitPos.uD); - if (average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM)) { + if (average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM)) { if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { + // if this fails, stick with the PM and/or PAR fit from above goto doneWithFit; } @@ -313,12 +314,22 @@ // I will not try to outlier-reject, just calculate the weighted average if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { + // if we have tried this, we have not masked any points; we will find a single unmasked point below goto doneWithFit; } - // XXX Maybe add a flag here? - average[0].flags |= ID_STAR_USE_AVE; + average[0].flags |= ID_OBJ_USE_AVE; + average[0].flags |= ID_OBJ_RAW_AVE; goto useBasic; } if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { - goto doneWithFit; + // if the above fails, we need to clear the masks and try again below + FitPointsClearMasks (fitStats->points, fitStats->Npoints); + // just calculate the weighted average + if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { + // if this fails, find a single unmasked point below + goto doneWithFit; + } + average[0].flags |= ID_OBJ_USE_AVE; + average[0].flags |= ID_OBJ_RAW_AVE; + goto useBasic; } // fprintf (stderr, "fit 2: %f %f : %f %f\n", fitPos.Ro, fitPos.Do, fitPos.uR, fitPos.uD); @@ -327,9 +338,12 @@ int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); if (Nnomask < POS_MIN_NPTS_BOOT) { + // if the above fails, we need to clear the masks and try again below + FitPointsClearMasks (fitStats->points, fitStats->Npoints); if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { + // if this fails, find a single unmasked point below goto doneWithFit; } - // XXX Maybe add a flag here? - average[0].flags |= ID_STAR_USE_AVE; + average[0].flags |= ID_OBJ_USE_AVE; + average[0].flags |= ID_OBJ_RAW_AVE; goto useBasic; } @@ -343,5 +357,5 @@ BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); } - average[0].flags |= ID_STAR_USE_AVE; + average[0].flags |= ID_OBJ_USE_AVE; } @@ -351,5 +365,5 @@ if (fabs(fitPos.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); - average[0].flags |= ID_STAR_FIT_AVE; + average[0].flags |= ID_OBJ_FIT_AVE; fitStats->Nave ++; } @@ -357,5 +371,5 @@ doneWithFit: // if no valid fit has been found, try to use a single unmasked point: - if (!(average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM | ID_STAR_USE_AVE))) { + if (!(average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM | ID_OBJ_USE_AVE))) { if (!FitPosPMfixed_Single (&fitPos, fitStats->points, fitStats->Npoints)) { fitStats->Nskip ++; @@ -366,6 +380,6 @@ if (fabs(fitPos.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); - average[0].flags |= ID_STAR_USE_AVE; - average[0].flags |= ID_STAR_FIT_AVE; + average[0].flags |= ID_OBJ_USE_AVE; + average[0].flags |= ID_OBJ_FIT_AVE; fitStats->Nave ++; } @@ -422,14 +436,14 @@ FitAstromResultInit (&fit); - if (average[0].flags & ID_STAR_USE_PAR) { - myAssert ((average[0].flags & (ID_STAR_USE_PM | ID_STAR_USE_AVE)) == 0, "programming error"); + if (average[0].flags & ID_OBJ_USE_PAR) { + myAssert ((average[0].flags & (ID_OBJ_USE_PM | ID_OBJ_USE_AVE)) == 0, "programming error"); fit = fitPar; } - if (average[0].flags & ID_STAR_USE_PM) { - myAssert ((average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_AVE)) == 0, "programming error"); + if (average[0].flags & ID_OBJ_USE_PM) { + myAssert ((average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_AVE)) == 0, "programming error"); fit = fitPM; } - if (average[0].flags & ID_STAR_USE_AVE) { - myAssert ((average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM)) == 0, "programming error"); + if (average[0].flags & ID_OBJ_USE_AVE) { + myAssert ((average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM)) == 0, "programming error"); fit = fitPos; } @@ -516,5 +530,5 @@ // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM) - average[0].flags &= ~ID_STAR_NO_ASTROM; + average[0].flags &= ~ID_OBJ_NO_MEAN_ASTROM; return (TRUE); @@ -554,4 +568,5 @@ average[0].dRstk = average[0].dR; // RA scatter in arcsec average[0].dDstk = average[0].dD; // DEC scatter in arcsec + average[0].flags |= ID_OBJ_MEAN_FOR_STACK; } @@ -563,11 +578,10 @@ // the stack positions are not statistically independent... - FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints); - FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE); + FitPosStack (&fitPos, fitStats); // project Ro, Do back to RA,DEC XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords); - // XXX choose stack flag? average[0].flags |= ID_STAR_FIT_AVE; + // XXX choose stack flag? average[0].flags |= ID_OBJ_FIT_AVE; fitStats->Nave ++; @@ -610,4 +624,6 @@ average[0].dRstk = fitPos.dRo; // RA scatter in arcsec average[0].dDstk = fitPos.dDo; // DEC scatter in arcsec + average[0].Tmean = (Tmean * 86400 * 365.25) + fitStats->T2000; + average[0].flags &= ~ID_OBJ_MEAN_FOR_STACK; return (TRUE); Index: /trunk/Ohana/src/relastro/src/bcatalog.c =================================================================== --- /trunk/Ohana/src/relastro/src/bcatalog.c (revision 39374) +++ /trunk/Ohana/src/relastro/src/bcatalog.c (revision 39375) @@ -25,4 +25,9 @@ int Nsecfilt; Coords coords; + + CoordTransform *transform = NULL; + if (USE_GALAXY_MODEL) { + transform = InitTransform (COORD_CELESTIAL, COORD_GALACTIC); + } /* for outlier rejection, project coordinates to a plane centered on the object with units of arcsec */ @@ -43,12 +48,17 @@ Nmeasure = Naverage = 0; - DVOAverageFlags averageBits = - ID_STAR_FIT_AVE | - ID_STAR_FIT_PM | - ID_STAR_FIT_PAR | - ID_STAR_USE_AVE | - ID_STAR_USE_PM | - ID_STAR_USE_PAR | - ID_STAR_NO_ASTROM ; + DVOAverageFlags astromBits = + ID_OBJ_LARGE_PM | // star with large proper motion + ID_OBJ_RAW_AVE | // simple weighted average position was used (no IRLS fitting) + ID_OBJ_FIT_AVE | // average position was fitted + ID_OBJ_FIT_PM | // proper motion model was fitted + ID_OBJ_FIT_PAR | // parallax model was fitted + ID_OBJ_USE_AVE | // average position used (not PM or PAR) + ID_OBJ_USE_PM | // proper motion used (not AVE or PAR) + ID_OBJ_USE_PAR | // parallax used (not AVE or PM) + ID_OBJ_NO_MEAN_ASTROM | // mean astrometry could not be measured + ID_OBJ_STACK_FOR_MEAN | // stack position used for mean astrometry + ID_OBJ_MEAN_FOR_STACK | // mean astrometry could not be measured + ID_OBJ_BAD_PM; // failure to measure proper-motion model if (VERBOSE2 && ExcludeBogus && (fbogus == NULL)) { @@ -78,9 +88,9 @@ if (RESET) { - subcatalog[0].average[Naverage].flags &= ~averageBits; + subcatalog[0].average[Naverage].flags &= ~astromBits; } // if we want to do static correction to ICRF, keep the ICRF QSOs (regardless of other data quality) - if (USE_ICRF_CORRECT && (subcatalog[0].average[Naverage].flags & ID_ICRF_QSO)) { + if (USE_ICRF_CORRECT && (subcatalog[0].average[Naverage].flags & ID_OBJ_ICRF_QSO)) { Nm = 0; for (j = 0; j < catalog[0].average[i].Nmeasure; j++) { @@ -115,19 +125,49 @@ } - // if we want to correct to the galay model, we need to use objects with distances determined by Green & Schlafly + // if we want to correct to the galaxy model, we need to use objects with distances determined by Green & Schlafly if (USE_GALAXY_MODEL) { // if we use the galaxy model, we need to require the model pm to exist - if (!catalog[0].average[i].Nstarpar) continue; - if (!isfinite(catalog[0].average[i].uRgal)) continue; - if (!isfinite(catalog[0].average[i].uDgal)) continue; + + int gotStarpar = TRUE; + gotStarpar = gotStarpar && (catalog[0].average[i].Nstarpar > 0); + gotStarpar = gotStarpar && isfinite(catalog[0].average[i].uRgal); + gotStarpar = gotStarpar && isfinite(catalog[0].average[i].uDgal); + if (!gotStarpar) { + // if we do not have starpar for some objects, generate an estimated starpar using d = 500pc + + // NOTE: distance is in kiloparsec + double distance = 0.5; // pow(10.0, 0.2*(starpar->DistMag + 5.0)) / 1000.0; + + double galLon, galLat; + ApplyTransform (&galLon, &galLat, catalog[0].average[i].R, catalog[0].average[i].D, transform); + + double Lrad = galLon * RAD_DEG; + double Brad = galLat * RAD_DEG; + + double uL_gal, uB_gal; + GalaxyMotionModel_radians(&uL_gal, &uB_gal, Lrad, Brad); + + double uL_sol, uB_sol; + SolarMotionModel_radians(&uL_sol, &uB_sol, Lrad, Brad, distance); + + double uL = uL_gal + uL_sol; + double uB = uB_gal + uB_sol; + + double uR, uD; + TransformProperMotion (&uR, &uD, uL, uB, catalog[0].average[i].R, catalog[0].average[i].D, transform); + + catalog[0].average[i].uRgal = uR; + catalog[0].average[i].uDgal = uD; + } else { - // filter stars which are not well-handled by the galactic motion model - int mStarpar = catalog[0].average[i].starparOffset; - myAssert (mStarpar < catalog[0].Nstarpar, "oops"); + // filter stars which are not well-handled by the galactic motion model + int mStarpar = catalog[0].average[i].starparOffset; + myAssert (mStarpar < catalog[0].Nstarpar, "oops"); - StarPar *starpar = &catalog[0].starpar[mStarpar]; - if (starpar->DistMag < MIN_DISTANCE_MOD) continue; - if (starpar->DistMag > MAX_DISTANCE_MOD) continue; - if (starpar->dDistMag > MAX_DISTANCE_MOD_ERR) continue; + StarPar *starpar = &catalog[0].starpar[mStarpar]; + if (starpar->DistMag < MIN_DISTANCE_MOD) continue; + if (starpar->DistMag > MAX_DISTANCE_MOD) continue; + if (starpar->dDistMag > MAX_DISTANCE_MOD_ERR) continue; + } } @@ -264,4 +304,5 @@ } + if (transform) free (transform); return (TRUE); } Index: /trunk/Ohana/src/relastro/src/fitobj2.c =================================================================== --- /trunk/Ohana/src/relastro/src/fitobj2.c (revision 39374) +++ /trunk/Ohana/src/relastro/src/fitobj2.c (revision 39375) @@ -136,9 +136,11 @@ double plx = 0.5*(drand48_cnt() + 0.0); - if (FIT_MODE == FIT_PM_AND_PAR) plx = 0.0; + if (FIT_MODE != FIT_PM_AND_PAR) plx = 0.0; int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 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 @@ -189,10 +191,10 @@ 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; + // 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"); @@ -250,4 +252,12 @@ 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; Index: /trunk/Ohana/src/relastro/src/fitstk.c =================================================================== --- /trunk/Ohana/src/relastro/src/fitstk.c (revision 39375) +++ /trunk/Ohana/src/relastro/src/fitstk.c (revision 39375) @@ -0,0 +1,270 @@ +# 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, int Nstack); + +# 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) + +static int Nsecfilt = 0; + +int main (int argc, char **argv) { + + int N_STARS = 3000; + int N_POINTS = 100; + int N_STACK = 6; + 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, "-Nstack"))) { + remove_argument (N, &argc, argv); + N_STACK = 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, " 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_cnt() - 0.5); + 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*drand48_cnt(); // minimum of 0 points + int Noutliers_max = Npoints * F_OUTLIERS; + int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points + int Nstack = N_STACK*drand48_cnt(); // minimum of 0 points + + if (Npoints + Nstack + 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, Nstack); + + UpdateObjects_Stack (catalog->average, catalog->secfilt, catalog->measureT, NULL, Nsecfilt, fitStats); + UpdateObjects_Chips (catalog->average, catalog->secfilt, catalog->measureT, NULL, Nsecfilt, fitStats, 0, 0); + + double Tmean = ohana_sec_to_mjd(catalog->average->Tmean); + double Tyear = (Tmean - MJD_J2000) / 365.25; + fprintf (stdout, "%3d %3d %3d %12.7f %12.7f %8.5f %8.5f %8.5f | %12.7f %12.7f %8.5f %8.5f %8.5f | %12.7f | %8.5f %8.5f %8.5f %8.5f %8.5f %3d | %10.5f %10.5f %8.5f %8.5f | %6.2f %6.2f %6.2f 0x%08x\n", + Npoints, Noutliers, Nstack, 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->Rstk, catalog->average->Dstk, catalog->average->dRstk, catalog->average->dDstk, + catalog->average->ChiSqAve, catalog->average->ChiSqPM, catalog->average->ChiSqPar, catalog->average->flags); + } + ohana_gaussdev_free(); + + exit (0); +} + +Catalog *mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad, int Nstack) { + + int i; + + int Ntotal = Npoints + Nbad + Nstack; + + ALLOCATE_PTR (catalog, Catalog, 1); + ALLOCATE (catalog->average, Average, 1); + ALLOCATE (catalog->secfilt, SecFilt, Nsecfilt); + ALLOCATE (catalog->measureT, MeasureTiny, Ntotal); + + 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; + + // 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 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 < Npoints + Nbad; 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); + + 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); + } + + // double oRstk = POS_ERROR*(drand48_cnt() - 0.5); + // double oDstk = POS_ERROR*(drand48_cnt() - 0.5); + + double oRstk = 0.0*(drand48_cnt() - 0.5); + double oDstk = 0.0*(drand48_cnt() - 0.5); + + for (i = Npoints + Nbad; 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); + + catalog->measureT[i].R = Ro + oRstk/3600.0/dcos(Do); + catalog->measureT[i].D = Do + oDstk/3600.0; + + catalog->measureT[i].dXccd = ToShortPixels(POS_ERROR); + catalog->measureT[i].dYccd = ToShortPixels(POS_ERROR); + + catalog->measureT[i].photcode = 11200; + + // 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; +} Index: /trunk/Ohana/src/relastro/src/hpm_objects.c =================================================================== --- /trunk/Ohana/src/relastro/src/hpm_objects.c (revision 39374) +++ /trunk/Ohana/src/relastro/src/hpm_objects.c (revision 39375) @@ -264,5 +264,5 @@ good |= (testcat.average[0].ChiSqPM < MAX_CHISQ_PM); good |= (testcat.average[0].Npos > MIN_NPOS_FRAC * testcat.average[0].Nmeasure); - good &= ((testcat.average[0].flags & ID_STAR_FIT_PM) > 0); + good &= ((testcat.average[0].flags & ID_OBJ_FIT_PM) > 0); if (good) { Index: /trunk/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- /trunk/Ohana/src/relastro/src/relastro_objects.c (revision 39374) +++ /trunk/Ohana/src/relastro/src/relastro_objects.c (revision 39375) @@ -15,12 +15,17 @@ ID_MEAS_USED_CHIP; - DVOAverageFlags averageBits = - ID_STAR_FIT_AVE | - ID_STAR_FIT_PM | - ID_STAR_FIT_PAR | - ID_STAR_USE_AVE | - ID_STAR_USE_PM | - ID_STAR_USE_PAR | - ID_STAR_NO_ASTROM ; + DVOAverageFlags astromBits = + ID_OBJ_LARGE_PM | // star with large proper motion + ID_OBJ_RAW_AVE | // simple weighted average position was used (no IRLS fitting) + ID_OBJ_FIT_AVE | // average position was fitted + ID_OBJ_FIT_PM | // proper motion model was fitted + ID_OBJ_FIT_PAR | // parallax model was fitted + ID_OBJ_USE_AVE | // average position used (not PM or PAR) + ID_OBJ_USE_PM | // proper motion used (not AVE or PAR) + ID_OBJ_USE_PAR | // parallax used (not AVE or PM) + ID_OBJ_NO_MEAN_ASTROM | // mean astrometry could not be measured + ID_OBJ_STACK_FOR_MEAN | // stack position used for mean astrometry + ID_OBJ_MEAN_FOR_STACK | // mean astrometry could not be measured + ID_OBJ_BAD_PM; // failure to measure proper-motion model // XXX need to decide how to determine PARALLEL mode... @@ -69,5 +74,5 @@ if (RESET) { for (j = 0; j < catalog.Naverage; j++) { - catalog.average[j].flags &= ~averageBits; + catalog.average[j].flags &= ~astromBits; m = catalog.average[j].measureOffset; for (k = 0; k < catalog.average[j].Nmeasure; k++) { Index: /trunk/Ohana/src/relastro/test/dvo.photcodes =================================================================== --- /trunk/Ohana/src/relastro/test/dvo.photcodes (revision 39374) +++ /trunk/Ohana/src/relastro/test/dvo.photcodes (revision 39375) @@ -32,2 +32,3 @@ 10233 GPC1.i.XY33 dep 24.611 -0.044 0.000 - - 0.0000 0 3 0.000 1.000 0.000 0.000 0x0000 0x0000 0xe0440130 0x1003bc88 + 11200 GPC1.i.SkyChip dep 25.000 0.000 0.000 - - 0.0000 0 3 0.000 1.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 Index: /trunk/Ohana/src/relastro/test/mana.sh =================================================================== --- /trunk/Ohana/src/relastro/test/mana.sh (revision 39374) +++ /trunk/Ohana/src/relastro/test/mana.sh (revision 39375) @@ -28,4 +28,24 @@ testnboot pm-boot testnboot plx-boot +end + +macro teststk + if ($0 != 3) + echo "USAGE: testobj2 (Nstars) (foutliers)" + break + end + + exec ../bin/fitstk.lin64 plx -Nstars $1 -foutliers $2 > test.stk.dat + load.stack test.stk.dat +end + +macro testobj2 + if ($0 != 3) + echo "USAGE: testobj2 (Nstars) (foutliers)" + break + end + + exec ../bin/fitobj2.lin64 plx -Nstars $1 -foutliers $2 > test.plx.dat + checkstats test.plx.dat end @@ -191,4 +211,66 @@ 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 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 + # determine the parallax factor (time) (ra) (dec) (parfR) (parfD) are vectors, macro getParFactor_vectors