Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/Makefile =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/Makefile (revision 38611) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/Makefile (revision 38612) @@ -26,8 +26,9 @@ $(SRC)/FitChip.$(ARCH).o \ $(SRC)/FitMosaic.$(ARCH).o \ +$(SRC)/FitSimple.$(ARCH).o \ +$(SRC)/FitAstromOps.$(ARCH).o \ $(SRC)/FitPM.$(ARCH).o \ $(SRC)/FitPMandPar.$(ARCH).o \ $(SRC)/FitPosPMfixed.$(ARCH).o \ -$(SRC)/FitSimple.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ @@ -99,8 +100,9 @@ RELASTRO_CLIENT = \ $(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/FitSimple.$(ARCH).o \ +$(SRC)/FitAstromOps.$(ARCH).o \ $(SRC)/FitPM.$(ARCH).o \ -$(SRC)/FitPar.$(ARCH).o \ $(SRC)/FitPMandPar.$(ARCH).o \ -$(SRC)/FitSimple.$(ARCH).o \ +$(SRC)/FitPosPMfixed.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/include/relastro.h =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/include/relastro.h (revision 38611) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/include/relastro.h (revision 38612) @@ -147,4 +147,5 @@ double C_blue; double C_red; + int measure; } FitAstromPoint; @@ -709,5 +710,5 @@ int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog); int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange); -int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff); +int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff, int isStack); int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats); int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int cat, off_t measOff); Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitAstromOps.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitAstromOps.c (revision 38611) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitAstromOps.c (revision 38612) @@ -4,4 +4,7 @@ FitStats *FitStatsInit (int Nmax, int Nboot) { + + FitStats *fitStats = NULL; + ALLOCATE (fitStats, FitStats, 1); // counters to record successful fits or failures @@ -27,6 +30,6 @@ fitStats->NpointsAlloc = Nmax; if (Nmax > 0) { - ALLOCATE (fitStats->points, FitAstromPoints, Nmax); - ALLOCATE (fitStats->sample, FitAstromPoints, Nmax); + ALLOCATE (fitStats->points, FitAstromPoint, Nmax); + ALLOCATE (fitStats->sample, FitAstromPoint, Nmax); } @@ -48,5 +51,5 @@ // use J2000 as a reference time fitStats->T2000 = ohana_date_to_sec ("2000/01/01,12:00:00"); - return; + return fitStats; } @@ -107,5 +110,5 @@ } -void FitAstromPointInit (FitAstromObject *object) { +void FitAstromPointInit (FitAstromPoint *object) { object->X = 0.0; object->Y = 0.0; @@ -122,4 +125,5 @@ object->C_blue = 0.0; object->C_red = 0.0; + object->measure= -1; return; } Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c (revision 38611) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c (revision 38612) @@ -1,5 +1,5 @@ # include "relastro.h" -int FitPMandPAR (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) { +int FitPMandPar (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) { int i; Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c (revision 38611) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c (revision 38612) @@ -66,6 +66,4 @@ int k; - int setRefColor = areImagesMatched(); - /* calculate the average value of R,D for a single star */ @@ -92,5 +90,5 @@ // select the measurements to be used in this analysis - UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, cat, measOff); + UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, cat, measOff, FALSE); // if there are no exposure detections, use the stack position @@ -107,5 +105,5 @@ double Tmean, Trange, parRange; - FitAstromPoints_Project (fitStats->points, fitStats->Npoints, &Tmean, &Trange, &parRange); + FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange); // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq @@ -120,5 +118,5 @@ goto justPosition; } - if (Npoints <= PM_TOOFEW) { + if (fitStats->Npoints <= PM_TOOFEW) { mode = FIT_AVERAGE; goto justPosition; @@ -161,9 +159,9 @@ // to fit PM, we do not attempt PAR. thus failure to fit PAR falls back to PM-only if (mode == FIT_PM_AND_PAR) { - if (N <= PAR_TOOFEW) { + if (parRange < PAR_FACTOR_MIN) { mode = FIT_PM_ONLY; goto justPosition; } - if (parRange < PAR_FACTOR_MIN) { + if (fitStats->Npoints <= PAR_TOOFEW) { mode = FIT_PM_ONLY; goto justPosition; @@ -220,4 +218,5 @@ for (k = 0; k < fitStats->Npoints; k++) { int Nm = fitStats->points[k].measure; + myAssert (Nm >= 0, "oops"); measure[Nm].dbFlags |= ID_MEAS_USED_OBJ; if (measureBig) { measureBig[Nm].dbFlags |= ID_MEAS_USED_OBJ; } @@ -225,4 +224,5 @@ # if (0) + int setRefColor = areImagesMatched(); if (setRefColor) { // need to reassign here if isfinite() @@ -345,6 +345,4 @@ int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats) { - off_t k; - // set the default values average[0].Rstk = NAN; // RA in degrees @@ -354,11 +352,8 @@ /* calculate the average value of R,D for a single star */ - PMFit fitAve; - memset (&fitAve, 0, sizeof(fitAve)); - fitAve.chisq = NAN; + FitAstromResult fitPos; + FitAstromResultInit (&fitPos); if (average[0].Nmeasure == 0) return TRUE; - - int N = 0; int XVERB = FALSE; @@ -366,69 +361,18 @@ XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST); - // find the basic properties of the detections for this object (Tmin, Tmax, Tmean) - for (k = 0; k < average[0].Nmeasure; k++) { - - if (XVERB) { - char *date = ohana_sec_to_date (measure[k].t); - int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0; - fprintf (stderr, "stack: "OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n", k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig); - free (date); - } - - // SKIP everything except gpc1 stack data - if (!isGPC1stack(measure[k].photcode)) continue; - - // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst - int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE); - if (!keepMeasure) { - continue; - } - - R[N] = getMeanR (&measure[k], average, secfilt); - D[N] = getMeanD (&measure[k], average, secfilt); - - // dX, dY : error in arcsec -- - dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA); - dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC); - - // allow a given photcode or measurement to be - // ignored if the error is NAN (for photcode, set astromErrSys to NaN) - if (isnan(dX[N])) continue; - if (isnan(dY[N])) continue; - - // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in - // arcseconds. dR, dD are the errors in those directions in degrees. IF we have - // non-circular errors (different values for X and Y), then dR and dD will be - // incorrect: they would need to be rotated to take out the position angle - dR[N] = dX[N] / 3600.0; - dD[N] = dY[N] / 3600.0; - - // XXX use a different flag for stack measurements? - // measure[k].dbFlags |= ID_MEAS_USED_OBJ; - // if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; } - - N++; - } // loop over measurements : average[0].Nmeasure - - // if we have too few good detections for the desired fit, or too limited a - // baseline, use a fit with fewer parameters. XXX if we have too few measurements - // for even the average position, consider including the lower-quality detections? + // select the measurements to be used in this analysis + UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, 0, 0, TRUE); // too few measurements for average position (require 2 values) - if (N < 1) return FALSE; // XXX ?? - - // find the mean position - StatType statsR, statsD; - liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here) - liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here) - - 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; + if (fitStats->Npoints < 1) return FALSE; // XXX ?? + + double Tmean, Trange, parRange; + FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange); + + FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints); + FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE); + + // 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; @@ -438,14 +382,14 @@ average[0].R, average[0].D, - fitAve.Ro, fitAve.Do, - 3600*(average[0].R - fitAve.Ro), - 3600*(average[0].D - fitAve.Do)); + fitPos.Ro, fitPos.Do, + 3600*(average[0].R - fitPos.Ro), + 3600*(average[0].D - fitPos.Do)); // make sure that the fit succeeded int status = TRUE; - status &= finite(fitAve.Ro); - status &= finite(fitAve.Do); - status &= finite(fitAve.dRo); - status &= finite(fitAve.dDo); + status &= finite(fitPos.Ro); + status &= finite(fitPos.Do); + status &= finite(fitPos.dRo); + status &= finite(fitPos.dDo); if (!status) { fitStats->Nskip ++; @@ -458,9 +402,9 @@ double dXoff, dYoff; - RD_to_XY (&dXoff, &dYoff, fitAve.Ro, fitAve.Do, &fitStats->coords); + RD_to_XY (&dXoff, &dYoff, fitPos.Ro, fitPos.Do, &fitStats->coords); float dPos = hypot (dXoff, dYoff); if (dPos > MaxMeanOffset) { if (fitStats->Noffset < 100) { - fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fitAve.Ro, fitAve.Do, dXoff, dYoff); + fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fitPos.Ro, fitPos.Do, dXoff, dYoff); } fitStats->Noffset ++; @@ -469,13 +413,13 @@ // set the stack position values - average[0].Rstk = fitAve.Ro; // RA in degrees - average[0].Dstk = fitAve.Do; // DEC in degrees - average[0].dRstk = fitAve.dRo; // RA scatter in arcsec - average[0].dDstk = fitAve.dDo; // DEC scatter in arcsec + average[0].Rstk = fitPos.Ro; // RA in degrees + average[0].Dstk = fitPos.Do; // DEC in degrees + average[0].dRstk = fitPos.dRo; // RA scatter in arcsec + average[0].dDstk = fitPos.dDo; // DEC scatter in arcsec return (TRUE); } -int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff) { +int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff, int isStack) { // I've already allocated fit->points (and fit->sample) with space for fit->NpointsAlloc entries @@ -495,9 +439,13 @@ } - // SKIP gpc1 stack data - if (isGPC1stack(measure[k].photcode)) continue; - // SKIP gpc1 forced-warp data if (isGPC1warp(measure[k].photcode)) continue; + + // SKIP gpc1 stack data + if (isStack) { + if (!isGPC1stack(measure[k].photcode)) continue; + } else { + if ( isGPC1stack(measure[k].photcode)) continue; + } // reset the bit to note that a detection was used (or not) @@ -518,8 +466,8 @@ // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius if (ExcludeBogus) { - fitStats->coords.crval1 = average[0].R; - fitStats->coords.crval2 = average[0].D; + fit->coords.crval1 = average[0].R; + fit->coords.crval2 = average[0].D; double Xi, Yi; - RD_to_XY (&Xi, &Yi, Ri, Di, &fitStats->coords); + RD_to_XY (&Xi, &Yi, Ri, Di, &fit->coords); double radius = hypot(Xi, Yi); if (radius > ExcludeBogusRadius) { @@ -565,4 +513,6 @@ points[Npoints].dT = measure[k].dt; + // XXX + int setRefColor = 0; if (setRefColor) { points[Npoints].C_blue = getColorBlue (measOff+k, cat); @@ -584,4 +534,5 @@ int k; + int Npoints = fitStats->Npoints; FitAstromPoint *points = fitStats->points; @@ -598,5 +549,5 @@ Tmin = MIN(Tmin, points[k].T); Tmax = MAX(Tmax, points[k].T); - Tmean += points[k].T; + *Tmean += points[k].T; // at this point, T is in years since J2000 @@ -627,5 +578,5 @@ for (k = 0; k < Npoints; k++) { RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords); - points[k].T -= Tmean; + points[k].T -= *Tmean; } return TRUE; @@ -636,5 +587,5 @@ int i, j; - Nmax = 0; + int Nmax = 0; for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { @@ -642,5 +593,4 @@ } } - return Nmax; } @@ -664,4 +614,6 @@ ALLOCATE (values, double, Nfit); + int i; + for (i = 0; i < Nfit; i++) { switch (mode) { @@ -745,4 +697,6 @@ int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode) { + int i; + // add up the chi square for the fit double chisq = 0.0; @@ -755,11 +709,11 @@ switch (mode) { case FIT_AVERAGE: - fit->chisq = chisq / (2.0*Npts - 2.0); + fit->chisq = chisq / (2.0*Npoints - 2.0); break; case FIT_PM_ONLY: - fit->chisq = chisq / (2.0*Npts - 4.0); + fit->chisq = chisq / (2.0*Npoints - 4.0); break; case FIT_PM_AND_PAR: - fit->chisq = chisq / (2.0*Npts - 5.0); + fit->chisq = chisq / (2.0*Npoints - 5.0); break; default: