Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/include/relastro.h =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/include/relastro.h (revision 38603) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/include/relastro.h (revision 38604) @@ -28,4 +28,12 @@ typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget; + +typedef enum { + FIT_RESULT_RA, + FIT_RESULT_DEC, + FIT_RESULT_uR, + FIT_RESULT_uD, + FIT_RESULT_PLX, +} FitAstromResultMode; typedef enum { @@ -119,5 +127,4 @@ double p, dp; - int getChisq; double chisq; int Nfit; @@ -140,5 +147,5 @@ double C_blue; double C_red; -} FitAstromPoints; +} FitAstromPoint; typedef struct { @@ -149,16 +156,17 @@ off_t Noffset; + double *values; FitAstromResult *fit; // use bootstrap resampling to generate Nfit fits to measure the stats int Nfit; int NfitAlloc; - FitAstromPoints *points; - FitAstromPoints *sample; + FitAstromPoint *points; + FitAstromPoint *sample; int Npoints; int NpointsAlloc; - FitAstromData *fitPos; - FitAstromData *fitPM; - FitAstromData *fitPar; + FitAstromData *fitdataPos; + FitAstromData *fitdataPM; + FitAstromData *fitdataPar; Coords coords; @@ -491,7 +499,7 @@ int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius); int ParFactor (double *pR, double *pD, double RA, double Dec, double Time, double Tmean); -int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB); -int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts); -int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int XVERB); +int FitPM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); +int FitPMandPar (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); +int FitPosPMfixed (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); Mosaic *getMosaicForImage (off_t N); @@ -694,2 +702,21 @@ int client_logger_init (char *dirname); int client_logger_message (char *format,...); + +int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode); +double VectorFractionInterpolate (double *values, float fraction, int Npts); +int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode); +int BootstrapResample (FitAstromPoints *sample, FitAstromPoints *points, int Npoints); +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_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); + +void FitAstromResultInit (FitAstromResult *fit); +void FitAstromPointInit (FitAstromObject *object); +void FitAstromDataFree (FitAstromData *fit); +FitAstromData *FitAstromDataInit (int Nterms); +void FitStatsFree (FitStats *fitStats); +void FitStatsSum (FitStats *src, FitStats *tgt); +void FitStatsReset (FitStats *tgt); +FitStats *FitStatsInit (int Nmax, int Nboot); Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitAstromOps.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitAstromOps.c (revision 38603) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitAstromOps.c (revision 38604) @@ -12,10 +12,11 @@ fitStats->Noffset = 0; - // container to hold the fit results (Nboot > 1 for bootstrap resampling) - fitStats->fit = NULL; + fitStats->values = NULL; // pre-allocated array for median & robust sigma + fitStats->fit = NULL; // container to hold the fit results (Nboot > 1 for bootstrap resampling) fitStats->Nfit = 0; fitStats->NfitAlloc = Nboot; if (Nboot > 0) { ALLOCATE (fitStats->fit, FitAstromResult, Nboot); + ALLOCATE (fitStats->values, double, Nboot); } @@ -31,12 +32,12 @@ // pre-allocated fit matrices for the 3 fit options - fitStats->fitPos = NULL; - fitStats->fitPM = NULL; - fitStats->fitPar = NULL; + fitStats->fitdataPos = NULL; + fitStats->fitdataPM = NULL; + fitStats->fitdataPar = NULL; if (Nmax > 0) { - fitStats->fitPos = FitAstromDataInit (2); - fitStats->fitPM = FitAstromDataInit (4); - fitStats->fitPar = FitAstromDataInit (5); + fitStats->fitdataPos = FitAstromDataInit (2); + fitStats->fitdataPM = FitAstromDataInit (4); + fitStats->fitdataPar = FitAstromDataInit (5); } @@ -72,10 +73,11 @@ FREE (fitStats->fit); + FREE (fitStats->values); FREE (fitStats->points); FREE (fitStats->sample); - FitAstromDataFree (fitStats->fitPos); - FitAstromDataFree (fitStats->fitPM); - FitAstromDataFree (fitStats->fitPar); + FitAstromDataFree (fitStats->fitdataPos); + FitAstromDataFree (fitStats->fitdataPM); + FitAstromDataFree (fitStats->fitdataPar); free (fitStats); @@ -105,7 +107,5 @@ } -} - -void FitAstromObjectInit (FitAstromObject *object) { +void FitAstromPointInit (FitAstromObject *object) { object->X = 0.0; object->Y = 0.0; @@ -125,23 +125,20 @@ } -void freeObjectData () { +void FitAstromResultInit (FitAstromResult *fit) { - free (R); - free (D); - free (T); - free (X); - free (Y); + fit->Ro = 0.0; + fit->dRo = 0.0; + fit->Do = 0.0; + fit->dDo = 0.0; + fit->uR = 0.0; + fit->duR = 0.0; + fit->uD = 0.0; + fit->duD = 0.0; + fit->p = 0.0; + fit->dp = 0.0; - free (dR); - free (dD); - free (dT); - free (dX); - free (dY); - - free (pX); - free (pY); - - free (C_blue); - free (C_red); -} - + fit->chisq = NAN; + fit->Nfit = 0; + + return; +} Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPM.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPM.c (revision 38603) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPM.c (revision 38604) @@ -65,17 +65,4 @@ fit->Nfit = Npoints; - // add up the chi square for the fit - if (fit->getChisq) { - double chisq = 0.0; - for (i = 0; i < Npoints; i++) { - double Xf = fit->Ro + fit->uR*points[i].T; - double Yf = fit->Do + fit->uD*points[i].T; - chisq += SQ(points[i].X - Xf) / SQ(points[i].dX); - chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY); - } - // the reduced chisq is divided by (Ndof = 2*Npts - 4) - fit->chisq = chisq / (2.0*Npts - 4.0); - } - return (TRUE); } Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c (revision 38603) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c (revision 38604) @@ -90,17 +90,4 @@ fit->Nfit = Npoints; - // add up the chi square for the fit - if (fit->getChisq) { - double chisq = 0.0; - for (i = 0; i < Npoints; i++) { - double Xf = fit->Ro + fit->uR*points[i].T + fit->p*points[i].pR; - double Yf = fit->Do + fit->uD*points[i].T + fit->p*points[i].pD; - chisq += SQ(points[i].X - Xf) / SQ(points[i].dX); - chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY); - } - // the reduced chisq is divided by (Ndof = 2*Npts - 5) - fit->chisq = chisq / (2.0*Npts - 5.0); - } - return (TRUE); } Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPosPMfixed.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPosPMfixed.c (revision 38603) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPosPMfixed.c (revision 38604) @@ -37,17 +37,4 @@ fit->Nfit = Npoints; - // add up the chi square for the fit - if (fit->getChisq) { - double chisq = 0.0; - for (i = 0; i < Npoints; i++) { - double Xf = fit->Ro + fit->uR*points[i].T; - double Yf = fit->Do + fit->uD*points[i].T; - chisq += SQ(points[i].X - Xf) / SQ(points[i].dX); - chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY); - } - // the reduced chisq is divided by (Ndof = 2*Npts - 4) - fit->chisq = chisq / (2.0*Npts - 4.0); - } - return (TRUE); } Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c (revision 38603) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c (revision 38604) @@ -2,7 +2,4 @@ # define PAR_TOOFEW 5 # define NBOOT 100 - -int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset); -int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats); // This function operates on both Measure and MeasureTiny. In the big stages, this should @@ -71,8 +68,8 @@ /* calculate the average value of R,D for a single star */ - PMFit fit; memset (&fit, 0, sizeof(fit)); - PMFit fitAve; memset (&fitAve, 0, sizeof(fitAve)); fitAve.chisq = NAN; - PMFit fitPM; memset (&fitPM, 0, sizeof(fitPM)); fitPM.chisq = NAN; - PMFit fitPAR; memset (&fitPAR, 0, sizeof(fitPAR)); fitPAR.chisq = NAN; + FitAstromResult fitPos, fitPM, fitPar; + FitAstromResultInit (&fitPos); + FitAstromResultInit (&fitPM); + FitAstromResultInit (&fitPar); // if we fail to fit the astrometry for some reason, we need to set/reset these @@ -96,8 +93,8 @@ // select the measurements to be used in this analysis - FitAstromObject *points = UpdateObjects_SelectMeasures (average, secfilt, measure, measureBig, cat, measOff, &Npoints); + UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, cat, measOff); // if there are no exposure detections, use the stack position - if (Npoints < 1) { + if (fitStats->Npoints < 1) { if (isfinite(average[0].Rstk) && isfinite(average[0].Dstk)) { average[0].R = average[0].Rstk; @@ -111,5 +108,5 @@ double Tmean, Trange, parRange; - UpdateObjects_Project (points, Npoints, &Tmean, &Trange, &parRange); + UpdateObjects_Project (fitStats->points, fitStats->Npoints, &Tmean, &Trange, &parRange); // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq @@ -129,13 +126,19 @@ } + // XXX handle the case with only 1 or few points + fitStats->Nfit = 0; for (k = 0; k < fitStats->NfitAlloc; k++) { - BootstrapResample (sample, points, Npoints); - if (!FitPM (&fitStats->fit[k], fitdata, sample, Npoints)) continue; + BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints); + // XXX init fitStats->fit[k]? set fitStats->fit[k].getChisq? + if (!FitPM (&fitStats->fit[k], fitStats->fitPM, fitStats->sample, fitStats->Npoints)) continue; fitStats->Nfit ++; } - // XXX convert the fit distributions of uR,uD to (uRo, uDo) +/- (duRo, duDo) - - if (XVERB) fprintf (stderr, "fitted PM: %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); + + FitAstromSetChisq (&fitPM, fitStats->points, fitStats->Npoints, FIT_PM_ONLY); // project Ro, Do back to RA,DEC @@ -171,20 +174,24 @@ fitStats->Nfit = 0; for (k = 0; k < fitStats->NfitAlloc; k++) { - BootstrapResample (sample, points, Npoints); - FitPMandPar (&fitStats->fit[k], fitsdata, sample, Npoints); + BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints); + FitPMandPar (&fitStats->fit[k], fitStats->fitPar, fitStats->sample, fitStats->Npoints); fitStats->Nfit ++; } - // XXX convert the fit distributions of uR,uD to (uRo, uDo) +/- (duRo, duDo) - - if (XVERB) fprintf (stderr, "fitted PM+PAR: %f - %f : %f %f : %f %f : %f %f : %f vs %f vs %f\n", Tmin, Tmax, fitPAR.Ro, fitPAR.Do, fitPAR.uR, fitPAR.uD, fitPAR.p, fitPAR.dp, fitPAR.chisq, fitPM.chisq, fitAve.chisq); - - XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX); + + FitAstromSetChisq (&fitPar, fitStats->points, fitStats->Npoints, FIT_PM_AND_PAR); + + XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &coords); average[0].flags |= ID_STAR_FIT_PAR; fitStats->Npar ++; - if (fabs(fitPAR.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); + if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); // XXX a hard-wired hack... - if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) { + if ((fabs(fitPar.uR) > 2.0) || (fabs(fitPar.uD) > 2.0)) { mode = FIT_PM_ONLY; } @@ -198,34 +205,26 @@ fitStats->Nfit = 0; for (k = 0; k < fitStats->NfitAlloc; k++) { - BootstrapResample (sample, points, Npoints); - FitPosPMfixed (&fitStats->fit[k], fitsdata, sample, Npoints); + BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints); + FitPosPMfixed (&fitStats->fit[k], fitStats->fitPos, fitStats->sample, fitStats->Npoints); fitStats->Nfit ++; } - - if (XVERB) fprintf (stderr, "average: %f %f\n", fitAve.Ro, fitAve.Do); + BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE); // project Ro, Do back to RA,DEC - XY_to_RD (&fitAve.Ro, &fitAve.Do, fitAve.Ro, fitAve.Do, &coords); + XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &coords); average[0].flags |= ID_STAR_FIT_AVE; fitStats->Nave ++; - - // XXX fitAve.dRo = statsR.sigma; - // XXX fitAve.dDo = statsD.sigma; - // XXX - // XXX fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN; - // XXX fitAve.Nfit = N; - // XXX - // XXX fitAve.uR = fitAve.duR = 0.0; - // XXX fitAve.uD = fitAve.duD = 0.0; - // XXX fitAve.p = fitAve.dp = 0.0; } // XXX update the bit flags of which points were used - for (k = 0; k < Npoints; k++) { - int Nm = points[k].measure; + for (k = 0; k < fitStats->Npoints; k++) { + int Nm = fitStats->points[k].measure; measure[Nm].dbFlags |= ID_MEAS_USED_OBJ; if (measureBig) { measureBig[Nm].dbFlags |= ID_MEAS_USED_OBJ; } } +# if (0) if (setRefColor) { // need to reassign here if isfinite() @@ -238,13 +237,14 @@ average[0].refColorRed = colorMedian; } +# endif /* choose the result based on the chisq values */ // XXXX for now, just use the mode as the result: int result = mode; - + FitAstromResult fit; switch (result) { case FIT_AVERAGE: average[0].flags |= ID_STAR_USE_AVE; - fit = fitAve; + fit = fitPos; break; case FIT_PM_ONLY: @@ -254,5 +254,5 @@ case FIT_PM_AND_PAR: average[0].flags |= ID_STAR_USE_PAR; - fit = fitPAR; + fit = fitPar; break; } @@ -314,5 +314,5 @@ average[0].uR, average[0].uD, - fitAve.chisq, fitPM.chisq, fitPAR.chisq); + fitPos.chisq, fitPM.chisq, fitPar.chisq); average[0].R = fit.Ro; // RA in degrees @@ -329,7 +329,7 @@ average[0].dP = fit.dp; // parallax error in arcsec - average[0].ChiSqAve = fitAve.chisq; + average[0].ChiSqAve = fitPos.chisq; average[0].ChiSqPM = fitPM.chisq; - average[0].ChiSqPar = fitPAR.chisq; + average[0].ChiSqPar = fitPar.chisq; average[0].Tmean = (Tmean * 86400 * 365.25) + T2000; @@ -479,22 +479,10 @@ } -int BootstrapResample (FitAstromPoints *sample, FitAstromPoints *points, int Npoints) { - int i; - - // I need to draw Npoints random entries from 'points' with replacement: - for (i = 0; i < Npoints; i++) { - int N = Npoints * drand48(); - sample[i] = points[N]; - } - return TRUE; -} - -FitAstromObject *UpdateObjects_SelectMeasures (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff, int *npoints) { - - // XXX move this above, but just for re-working now: - int Npoints = 0; - int NPOINTS = average->Nmeasure; - FitAstromObject *points = NULL; - ALLOCATE (points, FitAstromObject, NPOINTS); +int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff) { + + // I've already allocated fit->points (and fit->sample) with space for fit->NpointsAlloc entries + + int Npoints = fit->Npoints = 0; + FitAstromPoint *points = fit->points; // find the basic properties of the detections for this object (Tmin, Tmax, Tmean) @@ -551,5 +539,5 @@ } - FitAstromObjectInit (&points[Npoints]); + FitAstromPointInit (&points[Npoints]); points[Npoints].R = Ri; @@ -586,13 +574,17 @@ points[Npoints].measure = k; Npoints++; + + myAssert (Npoints <= fit->NpointsAlloc, "oops"); } // loop over measurements : average[0].Nmeasure - *npoints = Npoints; - return points; -} - -int FitAstromObject_Project (FitAstromObject *points, int Npoints, double *Tmean, double *Trange, double *parRange) { + fit->Npoints = Npoints; + return TRUE; +} + +int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange) { int k; + + FitAstromPoint *points = fitStats->points; // find Tmin & Tmax from the list of accepted measurements @@ -632,11 +624,11 @@ /* we need to do the fit in a locally linear space; choose a ref coordinate */ - coords.crval1 = points[0].R; - coords.crval2 = points[0].D; + fitStats->coords.crval1 = points[0].R; + fitStats->coords.crval2 = points[0].D; // project all of the R,D coordinates to a plane centered on this coordinate. set // the times to be relative to Tmean for (k = 0; k < Npoints; k++) { - RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &coords); + RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords); points[k].T -= Tmean; } @@ -658,27 +650,124 @@ } -/* fitting proper-motion and parallax: - - given a source at position r,d, at a time t, we need to calculate a vector (pr,pd) - - let x,y be the coordinate in the linearized frame with y parallel to DEC lines - - L,B are the ecliptic longitude and latitude of the object, - dL and dB are the offsets in the L and B directions - - dL = sin(t - topp) - dB = cos(t - topp)*sin(B) - - these need to be rotated to the R,D frame to yield pR,pD. Then, the equation of motion - for the source in the x,y frame is: - - x = Ro + uR * (t - to) + p * pR - y = Do + uD * (t - to) + p * pD - - the unknowns in these equations are Ro, uR, Do, uD, and p - - XXX think through the concepts for the pole a bit better. all objects near the pole - move the same way with the same phase. choose a projection center and define dL,dB relative - to that center point coordinate system? - -*/ +int BootstrapResample (FitAstromPoints *sample, FitAstromPoints *points, int Npoints) { + int i; + + // I need to draw Npoints random entries from 'points' with replacement: + for (i = 0; i < Npoints; i++) { + int N = Npoints * drand48(); + sample[i] = points[N]; + } + return TRUE; +} + +// calculate mean and sigma points for the 5 fit parameter +int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) { + + // generate a histogram for the selected element + double *values = NULL; + ALLOCATE (values, double, Nfit); + + for (i = 0; i < Nfit; i++) { + switch (mode) { + case FIT_RESULT_RA: + values[i] = fit[i].Ro; + break; + case FIT_RESULT_DEC: + values[i] = fit[i].Do; + break; + case FIT_RESULT_uR: + values[i] = fit[i].uR; + break; + case FIT_RESULT_uD: + values[i] = fit[i].uD; + break; + case FIT_RESULT_PLX: + values[i] = fit[i].p; + break; + default: + myAbort ("invalid option"); + } + } + + dsort (values, Nfit); + + double median; + if (Nfit % 2) { + int Ncenter = Nfit / 2; + median = values[Ncenter]; + } else { + int Ncenter = Nfit / 2 - 1; + median = 0.5*(values[Ncenter] + values[Ncenter + 1]); + } + + double Slo = VectorFractionInterpolate (values, 0.158655, Nfit); + double Shi = VectorFractionInterpolate (values, 0.841345, Nfit); + double sigma = (Shi - Slo) / 2.0; + + switch (mode) { + case FIT_RESULT_RA: + result->Ro = median; + result->dRo = sigma; + break; + case FIT_RESULT_DEC: + result->Do = median; + result->dDo = sigma; + break; + case FIT_RESULT_uR: + result->uR = median; + result->duR = sigma; + break; + case FIT_RESULT_uD: + result->uD = median; + result->duD = sigma; + break; + case FIT_RESULT_PLX: + result->p = median; + result->dp = sigma; + break; + default: + myAbort ("invalid option"); + } + + return TRUE; +} + +double VectorFractionInterpolate (double *values, float fraction, int Npts) { + + float F = fraction * Npts; + int N = fraction * Npts; + + if (N < 0 ) return NAN; + if (N >= Npts - 2) return NAN; + + // interpolate between N,N+1 + + double S = (F - N) * (values[N+1] - values[N]) + values[N]; + return S; +} + +int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode) { + + // add up the chi square for the fit + double chisq = 0.0; + for (i = 0; i < Npoints; i++) { + double Xf = fit->Ro + fit->uR*points[i].T + fit->p*points[i].pR; + double Yf = fit->Do + fit->uD*points[i].T + fit->p*points[i].pD; + chisq += SQ(points[i].X - Xf) / SQ(points[i].dX); + chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY); + } + switch (mode) { + case FIT_AVERAGE: + fit->chisq = chisq / (2.0*Npts - 2.0); + break; + case FIT_PM_ONLY: + fit->chisq = chisq / (2.0*Npts - 4.0); + break; + case FIT_PM_AND_PAR: + fit->chisq = chisq / (2.0*Npts - 5.0); + break; + default: + myAbort ("invalid mode"); + } + return (TRUE); +}