Index: /trunk/Ohana/src/opihi/cmd.astro/Makefile =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/Makefile (revision 39226) +++ /trunk/Ohana/src/opihi/cmd.astro/Makefile (revision 39227) @@ -37,4 +37,5 @@ $(SRC)/fitpm_irls.$(ARCH).o \ $(SRC)/fitplx_irls.$(ARCH).o \ +$(SRC)/astrom_ops.$(ARCH).o \ $(SRC)/fixwrap.$(ARCH).o \ $(SRC)/fixcols.$(ARCH).o \ Index: /trunk/Ohana/src/opihi/cmd.astro/astrom_ops.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/astrom_ops.c (revision 39227) +++ /trunk/Ohana/src/opihi/cmd.astro/astrom_ops.c (revision 39227) @@ -0,0 +1,72 @@ +# include "astro.h" +# define J2000 51544.5 /* Modified Julian date at standard epoch J2000 */ + +/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */ +int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius) { + + double n = mjd - J2000; // day number relative to standard epoch + double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration) + double g = (357.528 + 0.9856003 * n)*RAD_DEG; // Mean anomaly + + *lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees + *beta = 0.0; // approx latitude + *epsilon = (23.439 - 0.0000004 * n); // obliquity of ecliptic in degrees + *Radius = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU + return TRUE; +} + +/* given RA, DEC, Time, calculate the parallax factor */ +// RA,DEC are decimal degrees +// Time is MJD +int ParFactor (double *pR, double *pD, double RA, double DEC, double Time) { + + double lambda, beta, epsilon, Radius; + + // Time must be mjd + sun_ecliptic (Time, &lambda, &beta, &epsilon, &Radius); + + double lambda_rad = lambda*RAD_DEG; + double epsilon_rad = epsilon*RAD_DEG; + double RA_rad = RA*RAD_DEG; + double DEC_rad = DEC*RAD_DEG; + + double x = Radius*cos(lambda_rad); + double y = Radius*cos(epsilon_rad)*sin(lambda_rad); + double z = Radius*sin(epsilon_rad)*sin(lambda_rad); + + *pR = +(y*cos(RA_rad) - x*sin(RA_rad)); + *pD = -(y*sin(RA_rad) + x*cos(RA_rad))*sin(DEC_rad) + z*cos(DEC_rad); + + return TRUE; +} + +// allocate arrays but not the container +int PlxFitDataAlloc (PlxFitData *data, int N) { + + data->Npts = N; + ALLOCATE (data->X, double, N); + ALLOCATE (data->Y, double, N); + ALLOCATE (data->dX, double, N); + ALLOCATE (data->dY, double, N); + ALLOCATE (data->t, double, N); + ALLOCATE (data->pX, double, N); + ALLOCATE (data->pY, double, N); + ALLOCATE (data->Wx, double, N); + ALLOCATE (data->Wy, double, N); + ALLOCATE (data->index, int, N); + return TRUE; +} + +void PlxFitDataFree (PlxFitData *data) { + FREE (data->X); + FREE (data->Y); + FREE (data->dX); + FREE (data->dY); + FREE (data->t); + FREE (data->pX); + FREE (data->pY); + FREE (data->Wx); + FREE (data->Wy); + FREE (data->index); +} + Index: /trunk/Ohana/src/opihi/cmd.astro/fitplx.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/fitplx.c (revision 39226) +++ /trunk/Ohana/src/opihi/cmd.astro/fitplx.c (revision 39227) @@ -1,49 +1,3 @@ # include "astro.h" -# define J2000 51544.5 /* Modified Julian date at standard epoch J2000 */ - -# define ESCAPE(MSG,...) { \ - gprint (GP_ERR, MSG, __VA_ARGS__); \ - return FALSE; } - -typedef struct { - double *X; - double *Y; - double *t; - double *pX; - double *pY; - double *dX; - double *dY; - int *index; - int Npts; -} PlxFitData; - -typedef struct { - double Ro, dRo; - double Do, dDo; - - double uR, duR; - double uD, duD; - - double p, dp; - - double chisq; - int Nfit; - int getChisq; -} PlxFit; - -int VectorRobustStats (Vector *vector, double *median, double *sigma); -double VectorFractionInterpolate (double *values, float fraction, int Npts); - -int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal); -int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean); -int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE); - -int PlxFitDataAlloc (PlxFitData *data, int N); -void PlxFitDataFree (PlxFitData *data); -int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt); - -int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE); -int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius); -int ParFactor (double *pR, double *pD, double RA, double DEC, double Time); int fitplx (int argc, char **argv) { @@ -401,69 +355,4 @@ fit[0].Nfit = Npts; return (TRUE); -} - -/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */ -int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius) { - - double n = mjd - J2000; // day number relative to standard epoch - double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration) - double g = (357.528 + 0.9856003 * n)*RAD_DEG; // Mean anomaly - - *lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees - *beta = 0.0; // approx latitude - *epsilon = (23.439 - 0.0000004 * n); // obliquity of ecliptic in degrees - *Radius = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU - return TRUE; -} - -/* given RA, DEC, Time, calculate the parallax factor */ -// RA,DEC are decimal degrees -// Time is MJD -int ParFactor (double *pR, double *pD, double RA, double DEC, double Time) { - - double lambda, beta, epsilon, Radius; - - // Time must be mjd - sun_ecliptic (Time, &lambda, &beta, &epsilon, &Radius); - - double lambda_rad = lambda*RAD_DEG; - double epsilon_rad = epsilon*RAD_DEG; - double RA_rad = RA*RAD_DEG; - double DEC_rad = DEC*RAD_DEG; - - double x = Radius*cos(lambda_rad); - double y = Radius*cos(epsilon_rad)*sin(lambda_rad); - double z = Radius*sin(epsilon_rad)*sin(lambda_rad); - - *pR = +(y*cos(RA_rad) - x*sin(RA_rad)); - *pD = -(y*sin(RA_rad) + x*cos(RA_rad))*sin(DEC_rad) + z*cos(DEC_rad); - - return TRUE; -} - -// allocate arrays but not the container -int PlxFitDataAlloc (PlxFitData *data, int N) { - - data->Npts = N; - ALLOCATE (data->X, double, N); - ALLOCATE (data->Y, double, N); - ALLOCATE (data->dX, double, N); - ALLOCATE (data->dY, double, N); - ALLOCATE (data->t, double, N); - ALLOCATE (data->pX, double, N); - ALLOCATE (data->pY, double, N); - ALLOCATE (data->index, int, N); - return TRUE; -} - -void PlxFitDataFree (PlxFitData *data) { - FREE (data->X); - FREE (data->Y); - FREE (data->dX); - FREE (data->dY); - FREE (data->t); - FREE (data->pX); - FREE (data->pY); - FREE (data->index); } @@ -481,4 +370,8 @@ tgt->pX[i] = src->pX[N]; tgt->pY[i] = src->pY[N]; + + // *** make this optional? + tgt->Wx[i] = src->Wx[N]; + tgt->Wy[i] = src->Wy[N]; } return TRUE; @@ -539,4 +432,7 @@ pYmin = MIN (pYmin, fitdata->pY[Nsubset]); pYmax = MAX (pYmax, fitdata->pY[Nsubset]); + + fitdata->Wx[Nsubset] = 1.0; + fitdata->Wy[Nsubset] = 1.0; fitdata->index[Nsubset] = i; Nsubset++; Index: /trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c (revision 39226) +++ /trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c (revision 39227) @@ -1,66 +1,3 @@ # include "astro.h" -# define J2000 51544.5 /* Modified Julian date at standard epoch J2000 */ - -# define ESCAPE(MSG,...) { \ - gprint (GP_ERR, MSG, __VA_ARGS__); \ - return FALSE; } - -typedef struct { - double *X; - double *Y; - double *t; - double *pX; - double *pY; - double *dX; - double *dY; - - double *Wx; - double *Wy; - - int *index; - int Npts; -} PlxFitData; - -typedef struct { - double Ro, dRo; - double Do, dDo; - - double uR, duR; - double uD, duD; - - double p, dp; - - double chisq; - int Nfit; - int getChisq; -} PlxFit; - - -int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE); - -int VectorRobustStats_IRLS (Vector *vector, double *median, double *sigma); -double VectorFractionInterpolate_IRLS (double *values, float fraction, int Npts); - -int PlxSetMeanEpoch_IRLS (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal); -int PlxSetEpochPosition_IRLS (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean); -int PlxOutlierClip_IRLS (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE); - -int PlxFitDataAlloc_IRLS (PlxFitData *data, int N); -void PlxFitDataFree_IRLS (PlxFitData *data); -int PlxBootstrapResample_IRLS (PlxFitData *src, PlxFitData *tgt); - -int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *Wx, double *Wy, int Npts, int max_iterations, double outlier_limit, int VERBOSE); -int sun_ecliptic_IRLS (double mjd, double *lambda, double *beta, double *epsilon, double *Radius); -int ParFactor_IRLS (double *pR, double *pD, double RA, double DEC, double Time); - -int IRLS_converged (PlxFit *fit); -int Plx_weighted_LS (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts, - double **A, double **B, int VERBOSE); -double Plx_weight_cauchy (double x); -double Plx_dpsi_cauchy (double x); -double Plx_MAD(double *in, int N); - - - int fitplx_irls (int argc, char **argv) { @@ -176,5 +113,5 @@ double Rmean, Dmean, Tmean; - PlxSetMeanEpoch_IRLS (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal); + PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal); /* project coordinates to a plane centered on the object with units of arcsec */ @@ -186,6 +123,6 @@ PlxFitData fitdata; - PlxFitDataAlloc_IRLS (&fitdata, Ntotal); - PlxSetEpochPosition_IRLS (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean); + PlxFitDataAlloc (&fitdata, Ntotal); + PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean); PlxFit fit; memset (&fit, 0, sizeof(PlxFit)); @@ -193,13 +130,13 @@ if (0) { // This should be accomplished via the IRLS call - // determine dPsig for detections based on Noutlier attempts + // determine dPsig for detections based on Noutlier attempts if (Noutlier) { int clipRetry = TRUE; for (i = 0; clipRetry && (i < 3); i++) { - clipRetry = !PlxOutlierClip_IRLS (&fitdata, mask, Noutlier, dPsigMax, dPvec, VERBOSE); + clipRetry = !PlxOutlierClip (&fitdata, mask, Noutlier, dPsigMax, dPvec, VERBOSE); // using the new mask values, reset fitdata - PlxSetMeanEpoch_IRLS (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal); - PlxSetEpochPosition_IRLS (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean); + PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal); + PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean); if (VERBOSE) fprintf (stderr, "keep %d of %d\n", fitdata.Npts, Ntotal); } @@ -222,4 +159,5 @@ return FALSE; } + // Set the output mask based on the input mask and the outlier limits. if (outMask) { @@ -251,5 +189,5 @@ if (Nresample){ PlxFitData sample; - PlxFitDataAlloc_IRLS (&sample, fitdata.Npts); + PlxFitDataAlloc (&sample, fitdata.Npts); PlxFit *testfit = NULL; @@ -258,5 +196,5 @@ int Ngood = 0; for (i = 0; i < Nresample; i++) { - PlxBootstrapResample_IRLS (&fitdata, &sample); + PlxBootstrapResample (&fitdata, &sample); if (i % 100000 == 99999) fprintf (stderr, "."); @@ -296,9 +234,9 @@ // now calculate median and sigma for each vector - VectorRobustStats_IRLS (pvec, &fit.p, &fit.dp); - VectorRobustStats_IRLS (uRvec, &fit.uR, &fit.duR); - VectorRobustStats_IRLS (uDvec, &fit.uD, &fit.duD); - VectorRobustStats_IRLS (Rvec, &fit.Ro, &fit.dRo); - VectorRobustStats_IRLS (Dvec, &fit.Do, &fit.dDo); + VectorRobustStats (pvec, &fit.p, &fit.dp); + VectorRobustStats (uRvec, &fit.uR, &fit.duR); + VectorRobustStats (uDvec, &fit.uD, &fit.duD); + VectorRobustStats (Rvec, &fit.Ro, &fit.dRo); + VectorRobustStats (Dvec, &fit.Do, &fit.dDo); } @@ -588,340 +526,4 @@ } -/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */ -int sun_ecliptic_IRLS (double mjd, double *lambda, double *beta, double *epsilon, double *Radius) { - - double n = mjd - J2000; // day number relative to standard epoch - double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration) - double g = (357.528 + 0.9856003 * n)*RAD_DEG; // Mean anomaly - - *lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees - *beta = 0.0; // approx latitude - *epsilon = (23.439 - 0.0000004 * n); // obliquity of ecliptic in degrees - *Radius = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU - return TRUE; -} - -/* given RA, DEC, Time, calculate the parallax factor */ -// RA,DEC are decimal degrees -// Time is MJD -int ParFactor_IRLS (double *pR, double *pD, double RA, double DEC, double Time) { - - double lambda, beta, epsilon, Radius; - - // Time must be mjd - sun_ecliptic_IRLS (Time, &lambda, &beta, &epsilon, &Radius); - - double lambda_rad = lambda*RAD_DEG; - double epsilon_rad = epsilon*RAD_DEG; - double RA_rad = RA*RAD_DEG; - double DEC_rad = DEC*RAD_DEG; - - double x = Radius*cos(lambda_rad); - double y = Radius*cos(epsilon_rad)*sin(lambda_rad); - double z = Radius*sin(epsilon_rad)*sin(lambda_rad); - - *pR = +(y*cos(RA_rad) - x*sin(RA_rad)); - *pD = -(y*sin(RA_rad) + x*cos(RA_rad))*sin(DEC_rad) + z*cos(DEC_rad); - - return TRUE; -} - -// allocate arrays but not the container -int PlxFitDataAlloc_IRLS (PlxFitData *data, int N) { - - data->Npts = N; - ALLOCATE (data->X, double, N); - ALLOCATE (data->Y, double, N); - ALLOCATE (data->dX, double, N); - ALLOCATE (data->dY, double, N); - ALLOCATE (data->t, double, N); - ALLOCATE (data->pX, double, N); - ALLOCATE (data->pY, double, N); - ALLOCATE (data->Wx, double, N); - ALLOCATE (data->Wy, double, N); - ALLOCATE (data->index, int, N); - return TRUE; -} - -void PlxFitDataFree_IRLS (PlxFitData *data) { - FREE (data->X); - FREE (data->Y); - FREE (data->dX); - FREE (data->dY); - FREE (data->t); - FREE (data->pX); - FREE (data->pY); - FREE (data->Wx); - FREE (data->Wy); - FREE (data->index); -} - -int PlxBootstrapResample_IRLS (PlxFitData *src, PlxFitData *tgt) { - int i; - tgt->Npts = src->Npts; - for (i = 0; i < src->Npts; i++) { - int N = tgt->Npts * drand48(); - // int N = i; - tgt->X [i] = src->X [N]; - tgt->Y [i] = src->Y [N]; - tgt->dX[i] = src->dX[N]; - tgt->dY[i] = src->dY[N]; - tgt->t [i] = src->t [N]; - tgt->pX[i] = src->pX[N]; - tgt->pY[i] = src->pY[N]; - tgt->Wx[i] = src->Wx[N]; - tgt->Wy[i] = src->Wy[N]; - } - return TRUE; -} - -int PlxSetMeanEpoch_IRLS (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal) { - - int i; - - // find mean values to remove - double Nmean = 0; - *Tmean = 0; - *Rmean = 0; - *Dmean = 0; - double Tmin = +1000000; - double Tmax = -1000000; - for (i = 0; i < Ntotal; i++) { - if (mask && !mask[i]) continue; - *Rmean += R[i]; - *Dmean += D[i]; - *Tmean += T[i]; - Tmin = MIN(Tmin, T[i]); - Tmax = MAX(Tmax, T[i]); - Nmean += 1.0; - } - *Rmean /= Nmean; - *Dmean /= Nmean; - *Tmean /= Nmean; - - double Trange = Tmax - Tmin; - - // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", *Rmean, *Dmean, *Tmean, Trange, Tmin, Tmax); - - set_variable ("Trange", Trange); - return TRUE; -} - -// generate the fit values (projected X,Y; parallax factors; -int PlxSetEpochPosition_IRLS (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean) { - - int i; - - float pXmin = +2.0; - float pXmax = -2.0; - float pYmin = +2.0; - float pYmax = -2.0; - - int Nsubset = 0; - for (i = 0; i < Ntotal; i++) { - if (mask && !mask[i]) continue; - RD_to_XY (&fitdata->X[Nsubset], &fitdata->Y[Nsubset], R[i], D[i], coords); - fitdata->dX[Nsubset] = dR[i]; - fitdata->dY[Nsubset] = dD[i]; - fitdata->t[Nsubset] = (T[i] - Tmean) / 365.25; - ParFactor_IRLS (&fitdata->pX[Nsubset], &fitdata->pY[Nsubset], R[i], D[i], T[i]); - pXmin = MIN (pXmin, fitdata->pX[Nsubset]); - pXmax = MAX (pXmax, fitdata->pX[Nsubset]); - pYmin = MIN (pYmin, fitdata->pY[Nsubset]); - pYmax = MAX (pYmax, fitdata->pY[Nsubset]); - - fitdata->Wx[Nsubset] = 1.0; - fitdata->Wy[Nsubset] = 1.0; - fitdata->index[Nsubset] = i; - Nsubset++; - } - fitdata->Npts = Nsubset; - float dXRange = pXmax - pXmin; - float dYRange = pYmax - pYmin; - float parRange = hypot (dXRange, dYRange); - - set_variable ("Prange", parRange); - // fprintf (stderr, "par factor range: %f\n", parRange); - - return TRUE; -} - -/* Outlier clipping based on bootstrap-resampling tests of the plx path - * generate Noutlier resampled datasets - * fit the Noutlier plx paths - * determine and save the distribution of dXsig and dYsig for each point - * sort the resulting distributions and find dPsig (median point) for each measurement - * find the 90% point of dPsig : if > dPsigMax, only clip the 10% most deviant points - * set the dPvec values if desired - * -- mask is modified, dPvec values are set - * -- fitdata is unchanged - */ - -# define MAX_REJECT 0.1 - -int PlxOutlierClip_IRLS (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE) { - - int i, n; - - PlxFit testfit; - testfit.getChisq = FALSE; - - PlxFitData sample; - PlxFitDataAlloc_IRLS (&sample, fitdata->Npts); - - double **dXsig, **dYsig; - ALLOCATE (dXsig, double *, fitdata->Npts); - ALLOCATE (dYsig, double *, fitdata->Npts); - for (i = 0; i < fitdata->Npts; i++) { - ALLOCATE (dXsig[i], double, Noutlier); - ALLOCATE (dYsig[i], double, Noutlier); - } - - int Nsamples = 0; - for (n = 0; n < Noutlier; n++) { - // bootstrap resample (fitdata -> sample) - PlxBootstrapResample_IRLS (fitdata, &sample); - - if (n % 100000 == 99999) fprintf (stderr, "."); - - // fit the sample - if (!FitPMandPar (&testfit, - sample.X, sample.dX, - sample.Y, sample.dY, sample.t, - sample.pX, sample.pY, sample.Npts, VERBOSE)) continue; - - // fprintf (stderr, "%f +/- %f | %f %f\n", testfit.p, testfit.dp, testfit.uR, testfit.uD); - - // find the distances to the path - for (i = 0; i < fitdata->Npts; i++) { - double Xf = testfit.Ro + testfit.uR*fitdata->t[i] + testfit.p*fitdata->pX[i]; - double Yf = testfit.Do + testfit.uD*fitdata->t[i] + testfit.p*fitdata->pY[i]; - dXsig[i][Nsamples] = fabs(fitdata->X[i] - Xf) / fitdata->dX[i]; - dYsig[i][Nsamples] = fabs(fitdata->Y[i] - Yf) / fitdata->dY[i]; - // fprintf (stderr, "%f : %f %f : %f %f : %f %f : %f %f %f\n", T[i], Xf, Yf, fitdata->X[i], fitdata->Y[i], fitdata->dX[i], fitdata->dY[i], fitdata->t[i], fitdata->pX[i], fitdata->pY[i]); - } - Nsamples ++; - } - - double *dPsig; - ALLOCATE (dPsig, double, fitdata->Npts); - - for (i = 0; i < fitdata->Npts; i++) { - dsort (dXsig[i], Nsamples); - dsort (dYsig[i], Nsamples); - - // choose the median values - double dXsigMedian, dYsigMedian; - if (Nsamples % 2) { - int Ncenter = Nsamples / 2; - dXsigMedian = dXsig[i][Ncenter]; - dYsigMedian = dYsig[i][Ncenter]; - } else { - int Ncenter = Nsamples / 2 - 1; - dXsigMedian = 0.5*(dXsig[i][Ncenter] + dXsig[i][Ncenter + 1]); - dYsigMedian = 0.5*(dYsig[i][Ncenter] + dYsig[i][Ncenter + 1]); - } - // XXX replace with hypotenuse? - dPsig[i] = 0.5*(dXsigMedian + dYsigMedian); - // fprintf (stderr, "%d %10.6f %10.6f %10.6f %f %f : %f\n", i, R[i], D[i], T[i], dXsig[i][Ncenter], dYsig[i][Ncenter], dPsig[i]); - } - - // make a copy of dPsig[] and check if > 10% are > dPsigMax - double *dPsigSort; - ALLOCATE (dPsigSort, double, fitdata->Npts); - for (i = 0; i < fitdata->Npts; i++) { - dPsigSort[i] = dPsig[i]; - } - dsort (dPsigSort, fitdata->Npts); - int Nmax = (1.0 - MAX_REJECT)*fitdata->Npts; - - int completeClip = TRUE; - if (dPsigSort[Nmax] > dPsigMax) { - if (VERBOSE) fprintf (stderr, "too many outliers: %f at 90\n", dPsigSort[Nmax]); - dPsigMax = dPsigSort[Nmax]; - completeClip = FALSE; - } - - for (i = 0; i < fitdata->Npts; i++) { - if (dPsig[i] < dPsigMax) continue; - int n = fitdata->index[i]; - // fprintf (stderr, "clip %d: %f : %f\n", i, fitdata->t[i], dPsig[i]); - mask[n] = 0; // mask these points - } - - // only set dPvec if we have completed the clipping? - if (dPvec) { - for (i = 0; i < dPvec->Nelements; i++) { - dPvec->elements.Flt[i] = NAN; - } - for (i = 0; i < fitdata->Npts; i++) { - int n = fitdata->index[i]; - dPvec->elements.Flt[n] = dPsig[i]; - } - } - - free (dPsig); - free (dPsigSort); - - for (i = 0; i < fitdata->Npts; i++) { - free (dXsig[i]); - free (dYsig[i]); - } - free (dXsig); - free (dYsig); - - return completeClip; -} - -int VectorRobustStats_IRLS (Vector *vector, double *median, double *sigma) { - - // warn if vector->Nelements > 1000? 10000?) - // warn if vector is not float - - // we need to copy the vector to avoid changing the sort order - double *values = NULL; - ALLOCATE (values, double, vector->Nelements); - - int i; - int Npts = 0; - for (i = 0; i < vector->Nelements; i++) { - if (!isfinite(vector->elements.Flt[i])) continue; - values[Npts] = vector->elements.Flt[i]; - Npts++; - } - - dsort (values, Npts); - - if (Npts % 2) { - int Ncenter = Npts / 2; - *median = values[Ncenter]; - } else { - int Ncenter = Npts / 2 - 1; - *median = 0.5*(values[Ncenter] + values[Ncenter + 1]); - } - - double Slo = VectorFractionInterpolate_IRLS (values, 0.158655, Npts); - double Shi = VectorFractionInterpolate_IRLS (values, 0.841345, Npts); - - *sigma = (Shi - Slo) / 2.0; - - return TRUE; -} - -double VectorFractionInterpolate_IRLS (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; -} - - double Plx_weight_cauchy (double x) { double r = x / 2.385; Index: /trunk/Ohana/src/opihi/cmd.data/read_vectors.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/read_vectors.c (revision 39226) +++ /trunk/Ohana/src/opihi/cmd.data/read_vectors.c (revision 39227) @@ -348,4 +348,5 @@ } +# undef ESCAPE # define ESCAPE(...) { \ gprint (GP_ERR, __VA_ARGS__); \ Index: /trunk/Ohana/src/opihi/cmd.data/reindex.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/reindex.c (revision 39226) +++ /trunk/Ohana/src/opihi/cmd.data/reindex.c (revision 39227) @@ -1,5 +1,3 @@ # include "data.h" - -# define ESCAPE(MSG,...){ gprint (GP_ERR, MSG, __VA_ARGS__); goto error; } int reindex (int argc, char **argv) { Index: /trunk/Ohana/src/opihi/dvo/coordmosaic.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/coordmosaic.c (revision 39226) +++ /trunk/Ohana/src/opihi/dvo/coordmosaic.c (revision 39227) @@ -224,4 +224,5 @@ } +# undef ESCAPE # define ESCAPE(MSG) { \ if (src) fclose (src); \ Index: /trunk/Ohana/src/opihi/include/astro.h =================================================================== --- /trunk/Ohana/src/opihi/include/astro.h (revision 39226) +++ /trunk/Ohana/src/opihi/include/astro.h (revision 39227) @@ -12,3 +12,57 @@ void FreeAstro (void); +typedef struct { + double *X; + double *Y; + double *t; + double *pX; + double *pY; + double *dX; + double *dY; + double *Wx; + double *Wy; + int *index; + int Npts; +} PlxFitData; + +typedef struct { + double Ro, dRo; + double Do, dDo; + + double uR, duR; + double uD, duD; + + double p, dp; + + double chisq; + int Nfit; + int getChisq; +} PlxFit; + +int VectorRobustStats (Vector *vector, double *median, double *sigma); +double VectorFractionInterpolate (double *values, float fraction, int Npts); + +int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal); +int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean); +int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE); + +int PlxFitDataAlloc (PlxFitData *data, int N); +void PlxFitDataFree (PlxFitData *data); +int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt); + +int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE); +int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius); +int ParFactor (double *pR, double *pD, double RA, double DEC, double Time); + +/***** */ + +int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *Wx, double *Wy, int Npts, int max_iterations, double outlier_limit, int VERBOSE); + +// ?? int IRLS_converged (PlxFit *fit); +int Plx_weighted_LS (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts, + double **A, double **B, int VERBOSE); +double Plx_weight_cauchy (double x); +double Plx_dpsi_cauchy (double x); +double Plx_MAD(double *in, int N); + # endif Index: /trunk/Ohana/src/opihi/include/shell.h =================================================================== --- /trunk/Ohana/src/opihi/include/shell.h (revision 39226) +++ /trunk/Ohana/src/opihi/include/shell.h (revision 39227) @@ -4,4 +4,6 @@ # ifndef SHELL_H # define SHELL_H + +# define ESCAPE(MSG,...) { gprint (GP_ERR, MSG, __VA_ARGS__); return FALSE; } # define ISVAR(a) (isalnum (a) || (a == ':') || (a == '_'))