IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39227


Ignore:
Timestamp:
Dec 3, 2015, 5:05:24 PM (11 years ago)
Author:
eugene
Message:

merge common code between fitplx and fitplx_irls; move ESCAPE to shell.h

Location:
trunk/Ohana/src/opihi
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.astro/Makefile

    r39225 r39227  
    3737$(SRC)/fitpm_irls.$(ARCH).o  \
    3838$(SRC)/fitplx_irls.$(ARCH).o  \
     39$(SRC)/astrom_ops.$(ARCH).o        \
    3940$(SRC)/fixwrap.$(ARCH).o           \
    4041$(SRC)/fixcols.$(ARCH).o           \

  • trunk/Ohana/src/opihi/cmd.astro/fitplx.c

    r38986 r39227  
    11# include "astro.h"
    2 # define J2000 51544.5       /* Modified Julian date at standard epoch J2000 */
    3 
    4 # define ESCAPE(MSG,...) {                      \
    5     gprint (GP_ERR, MSG, __VA_ARGS__);          \
    6     return FALSE; }
    7 
    8 typedef struct {
    9   double *X;
    10   double *Y;
    11   double *t;
    12   double *pX;
    13   double *pY;
    14   double *dX;
    15   double *dY;
    16   int *index;
    17   int Npts;
    18 } PlxFitData;
    19 
    20 typedef struct {
    21   double Ro, dRo;
    22   double Do, dDo;
    23 
    24   double uR, duR;
    25   double uD, duD;
    26 
    27   double p, dp;
    28 
    29   double chisq;
    30   int Nfit;
    31   int getChisq;
    32 } PlxFit;
    33 
    34 int VectorRobustStats (Vector *vector, double *median, double *sigma);
    35 double VectorFractionInterpolate (double *values, float fraction, int Npts);
    36 
    37 int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal);
    38 int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean);
    39 int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE);
    40 
    41 int PlxFitDataAlloc (PlxFitData *data, int N);
    42 void PlxFitDataFree (PlxFitData *data);
    43 int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt);
    44 
    45 int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE);
    46 int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius);
    47 int ParFactor (double *pR, double *pD, double RA, double DEC, double Time);
    482
    493int fitplx (int argc, char **argv) {
    … …  
    401355  fit[0].Nfit = Npts;
    402356  return (TRUE);
    403 }
    404 
    405 /* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */
    406 int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius) {
    407 
    408   double n = mjd - J2000;                          // day number relative to standard epoch
    409   double L = 280.460 + 0.9856474 * n;              // mean solar longitute (corr. for aberration)
    410   double g = (357.528 + 0.9856003 * n)*RAD_DEG;    // Mean anomaly
    411 
    412   *lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees
    413   *beta = 0.0;                                     // approx latitude
    414   *epsilon = (23.439 - 0.0000004 * n);             // obliquity of ecliptic in degrees
    415   *Radius = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU
    416   return TRUE;
    417 }
    418 
    419 /* given RA, DEC, Time, calculate the parallax factor */
    420 // RA,DEC are decimal degrees
    421 // Time is MJD
    422 int ParFactor (double *pR, double *pD, double RA, double DEC, double Time) {
    423 
    424   double lambda, beta, epsilon, Radius;
    425 
    426   // Time must be mjd
    427   sun_ecliptic (Time, &lambda, &beta, &epsilon, &Radius);
    428 
    429   double lambda_rad = lambda*RAD_DEG;
    430   double epsilon_rad = epsilon*RAD_DEG;
    431   double RA_rad = RA*RAD_DEG;
    432   double DEC_rad = DEC*RAD_DEG;
    433 
    434   double x = Radius*cos(lambda_rad);
    435   double y = Radius*cos(epsilon_rad)*sin(lambda_rad);
    436   double z = Radius*sin(epsilon_rad)*sin(lambda_rad);
    437 
    438   *pR = +(y*cos(RA_rad) - x*sin(RA_rad));
    439   *pD = -(y*sin(RA_rad) + x*cos(RA_rad))*sin(DEC_rad) + z*cos(DEC_rad);
    440 
    441   return TRUE;
    442 }
    443 
    444 // allocate arrays but not the container
    445 int PlxFitDataAlloc (PlxFitData *data, int N) {
    446 
    447   data->Npts = N;
    448   ALLOCATE (data->X, double, N);
    449   ALLOCATE (data->Y, double, N);
    450   ALLOCATE (data->dX, double, N);
    451   ALLOCATE (data->dY, double, N);
    452   ALLOCATE (data->t, double, N);
    453   ALLOCATE (data->pX, double, N);
    454   ALLOCATE (data->pY, double, N);
    455   ALLOCATE (data->index, int, N);
    456   return TRUE;
    457 }
    458 
    459 void PlxFitDataFree (PlxFitData *data) {
    460   FREE (data->X);
    461   FREE (data->Y);
    462   FREE (data->dX);
    463   FREE (data->dY);
    464   FREE (data->t);
    465   FREE (data->pX);
    466   FREE (data->pY);
    467   FREE (data->index);
    468357}
    469358
    … …  
    481370    tgt->pX[i] = src->pX[N];
    482371    tgt->pY[i] = src->pY[N];
     372
     373    // *** make this optional?
     374    tgt->Wx[i] = src->Wx[N];
     375    tgt->Wy[i] = src->Wy[N];
    483376  }
    484377  return TRUE;
    … …  
    539432    pYmin = MIN (pYmin, fitdata->pY[Nsubset]);
    540433    pYmax = MAX (pYmax, fitdata->pY[Nsubset]);
     434
     435    fitdata->Wx[Nsubset] = 1.0;
     436    fitdata->Wy[Nsubset] = 1.0;   
    541437    fitdata->index[Nsubset] = i;
    542438    Nsubset++;

  • trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c

    r39226 r39227  
    11# include "astro.h"
    2 # define J2000 51544.5       /* Modified Julian date at standard epoch J2000 */
    3 
    4 # define ESCAPE(MSG,...) {                      \
    5     gprint (GP_ERR, MSG, __VA_ARGS__);          \
    6     return FALSE; }
    7 
    8 typedef struct {
    9   double *X;
    10   double *Y;
    11   double *t;
    12   double *pX;
    13   double *pY;
    14   double *dX;
    15   double *dY;
    16 
    17   double *Wx;
    18   double *Wy;
    19  
    20   int *index;
    21   int Npts;
    22 } PlxFitData;
    23 
    24 typedef struct {
    25   double Ro, dRo;
    26   double Do, dDo;
    27 
    28   double uR, duR;
    29   double uD, duD;
    30 
    31   double p, dp;
    32 
    33   double chisq;
    34   int Nfit;
    35   int getChisq;
    36 } PlxFit;
    37 
    38 
    39 int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE);
    40 
    41 int VectorRobustStats_IRLS (Vector *vector, double *median, double *sigma);
    42 double VectorFractionInterpolate_IRLS (double *values, float fraction, int Npts);
    43 
    44 int PlxSetMeanEpoch_IRLS (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal);
    45 int PlxSetEpochPosition_IRLS (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean);
    46 int PlxOutlierClip_IRLS (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE);
    47 
    48 int PlxFitDataAlloc_IRLS (PlxFitData *data, int N);
    49 void PlxFitDataFree_IRLS (PlxFitData *data);
    50 int PlxBootstrapResample_IRLS (PlxFitData *src, PlxFitData *tgt);
    51 
    52 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);
    53 int sun_ecliptic_IRLS (double mjd, double *lambda, double *beta, double *epsilon, double *Radius);
    54 int ParFactor_IRLS (double *pR, double *pD, double RA, double DEC, double Time);
    55 
    56 int IRLS_converged (PlxFit *fit);
    57 int Plx_weighted_LS (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts,
    58                      double **A, double **B, int VERBOSE);
    59 double Plx_weight_cauchy (double x);
    60 double Plx_dpsi_cauchy (double x);
    61 double Plx_MAD(double *in, int N);
    62 
    63 
    64 
    652
    663int fitplx_irls (int argc, char **argv) {
    … …  
    176113 
    177114  double Rmean, Dmean, Tmean;
    178   PlxSetMeanEpoch_IRLS (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
     115  PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
    179116
    180117  /* project coordinates to a plane centered on the object with units of arcsec */
    … …  
    186123
    187124  PlxFitData fitdata;
    188   PlxFitDataAlloc_IRLS (&fitdata, Ntotal);
    189   PlxSetEpochPosition_IRLS (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
     125  PlxFitDataAlloc (&fitdata, Ntotal);
     126  PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
    190127
    191128  PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
    … …  
    193130  if (0) {
    194131    // This should be accomplished via the IRLS call
    195   // determine dPsig for detections based on Noutlier attempts
     132    // determine dPsig for detections based on Noutlier attempts
    196133    if (Noutlier) {
    197134      int clipRetry = TRUE;
    198135      for (i = 0; clipRetry && (i < 3); i++) {
    199         clipRetry = !PlxOutlierClip_IRLS (&fitdata, mask, Noutlier, dPsigMax, dPvec, VERBOSE);
     136        clipRetry = !PlxOutlierClip (&fitdata, mask, Noutlier, dPsigMax, dPvec, VERBOSE);
    200137       
    201138        // using the new mask values, reset fitdata
    202         PlxSetMeanEpoch_IRLS (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
    203         PlxSetEpochPosition_IRLS (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
     139        PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
     140        PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
    204141        if (VERBOSE) fprintf (stderr, "keep %d of %d\n", fitdata.Npts, Ntotal);
    205142      }
    … …  
    222159    return FALSE;
    223160  }
     161
    224162  // Set the output mask based on the input mask and the outlier limits.
    225163  if (outMask) {
    … …  
    251189  if (Nresample){
    252190    PlxFitData sample;
    253     PlxFitDataAlloc_IRLS (&sample, fitdata.Npts);
     191    PlxFitDataAlloc (&sample, fitdata.Npts);
    254192
    255193    PlxFit *testfit = NULL;
    … …  
    258196    int Ngood = 0;
    259197    for (i = 0; i < Nresample; i++) {
    260       PlxBootstrapResample_IRLS (&fitdata, &sample);
     198      PlxBootstrapResample (&fitdata, &sample);
    261199     
    262200      if (i % 100000 == 99999) fprintf (stderr, ".");
    … …  
    296234
    297235    // now calculate median and sigma for each vector
    298     VectorRobustStats_IRLS (pvec,  &fit.p,  &fit.dp);
    299     VectorRobustStats_IRLS (uRvec, &fit.uR, &fit.duR);
    300     VectorRobustStats_IRLS (uDvec, &fit.uD, &fit.duD);
    301     VectorRobustStats_IRLS (Rvec,  &fit.Ro, &fit.dRo);
    302     VectorRobustStats_IRLS (Dvec,  &fit.Do, &fit.dDo);
     236    VectorRobustStats (pvec,  &fit.p,  &fit.dp);
     237    VectorRobustStats (uRvec, &fit.uR, &fit.duR);
     238    VectorRobustStats (uDvec, &fit.uD, &fit.duD);
     239    VectorRobustStats (Rvec,  &fit.Ro, &fit.dRo);
     240    VectorRobustStats (Dvec,  &fit.Do, &fit.dDo);
    303241  }
    304242
    … …  
    588526}
    589527
    590 /* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */
    591 int sun_ecliptic_IRLS (double mjd, double *lambda, double *beta, double *epsilon, double *Radius) {
    592 
    593   double n = mjd - J2000;                          // day number relative to standard epoch
    594   double L = 280.460 + 0.9856474 * n;              // mean solar longitute (corr. for aberration)
    595   double g = (357.528 + 0.9856003 * n)*RAD_DEG;    // Mean anomaly
    596 
    597   *lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees
    598   *beta = 0.0;                                     // approx latitude
    599   *epsilon = (23.439 - 0.0000004 * n);             // obliquity of ecliptic in degrees
    600   *Radius = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU
    601   return TRUE;
    602 }
    603 
    604 /* given RA, DEC, Time, calculate the parallax factor */
    605 // RA,DEC are decimal degrees
    606 // Time is MJD
    607 int ParFactor_IRLS (double *pR, double *pD, double RA, double DEC, double Time) {
    608 
    609   double lambda, beta, epsilon, Radius;
    610 
    611   // Time must be mjd
    612   sun_ecliptic_IRLS (Time, &lambda, &beta, &epsilon, &Radius);
    613 
    614   double lambda_rad = lambda*RAD_DEG;
    615   double epsilon_rad = epsilon*RAD_DEG;
    616   double RA_rad = RA*RAD_DEG;
    617   double DEC_rad = DEC*RAD_DEG;
    618 
    619   double x = Radius*cos(lambda_rad);
    620   double y = Radius*cos(epsilon_rad)*sin(lambda_rad);
    621   double z = Radius*sin(epsilon_rad)*sin(lambda_rad);
    622 
    623   *pR = +(y*cos(RA_rad) - x*sin(RA_rad));
    624   *pD = -(y*sin(RA_rad) + x*cos(RA_rad))*sin(DEC_rad) + z*cos(DEC_rad);
    625 
    626   return TRUE;
    627 }
    628 
    629 // allocate arrays but not the container
    630 int PlxFitDataAlloc_IRLS (PlxFitData *data, int N) {
    631 
    632   data->Npts = N;
    633   ALLOCATE (data->X, double, N);
    634   ALLOCATE (data->Y, double, N);
    635   ALLOCATE (data->dX, double, N);
    636   ALLOCATE (data->dY, double, N);
    637   ALLOCATE (data->t, double, N);
    638   ALLOCATE (data->pX, double, N);
    639   ALLOCATE (data->pY, double, N);
    640   ALLOCATE (data->Wx, double, N);
    641   ALLOCATE (data->Wy, double, N);
    642   ALLOCATE (data->index, int, N);
    643   return TRUE;
    644 }
    645 
    646 void PlxFitDataFree_IRLS (PlxFitData *data) {
    647   FREE (data->X);
    648   FREE (data->Y);
    649   FREE (data->dX);
    650   FREE (data->dY);
    651   FREE (data->t);
    652   FREE (data->pX);
    653   FREE (data->pY);
    654   FREE (data->Wx);
    655   FREE (data->Wy);
    656   FREE (data->index);
    657 }
    658 
    659 int PlxBootstrapResample_IRLS (PlxFitData *src, PlxFitData *tgt) {
    660   int i;
    661   tgt->Npts = src->Npts;
    662   for (i = 0; i < src->Npts; i++) {
    663     int N = tgt->Npts * drand48();
    664     // int N = i;
    665     tgt->X [i] = src->X [N];
    666     tgt->Y [i] = src->Y [N];
    667     tgt->dX[i] = src->dX[N];
    668     tgt->dY[i] = src->dY[N];
    669     tgt->t [i] = src->t [N];
    670     tgt->pX[i] = src->pX[N];
    671     tgt->pY[i] = src->pY[N];
    672     tgt->Wx[i] = src->Wx[N];
    673     tgt->Wy[i] = src->Wy[N];
    674   }
    675   return TRUE;
    676 }
    677 
    678 int PlxSetMeanEpoch_IRLS (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal) {
    679 
    680   int i;
    681 
    682   // find mean values to remove
    683   double Nmean = 0;
    684   *Tmean = 0;
    685   *Rmean = 0;
    686   *Dmean = 0;
    687   double Tmin = +1000000;
    688   double Tmax = -1000000;
    689   for (i = 0; i < Ntotal; i++) {
    690     if (mask && !mask[i]) continue;
    691     *Rmean += R[i];
    692     *Dmean += D[i];
    693     *Tmean += T[i];
    694     Tmin = MIN(Tmin, T[i]);
    695     Tmax = MAX(Tmax, T[i]);
    696     Nmean += 1.0;
    697   }
    698   *Rmean /= Nmean;
    699   *Dmean /= Nmean;
    700   *Tmean /= Nmean;
    701  
    702   double Trange = Tmax - Tmin;
    703 
    704   // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", *Rmean, *Dmean, *Tmean, Trange, Tmin, Tmax);
    705 
    706   set_variable ("Trange", Trange);
    707   return TRUE;
    708 }
    709 
    710 // generate the fit values (projected X,Y; parallax factors;
    711 int PlxSetEpochPosition_IRLS (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean) {
    712 
    713   int i;
    714 
    715   float pXmin = +2.0;
    716   float pXmax = -2.0;
    717   float pYmin = +2.0;
    718   float pYmax = -2.0;
    719 
    720   int Nsubset = 0;
    721   for (i = 0; i < Ntotal; i++) {
    722     if (mask && !mask[i]) continue;
    723     RD_to_XY (&fitdata->X[Nsubset], &fitdata->Y[Nsubset], R[i], D[i], coords);
    724     fitdata->dX[Nsubset] = dR[i];
    725     fitdata->dY[Nsubset] = dD[i];
    726     fitdata->t[Nsubset] = (T[i] - Tmean) / 365.25;
    727     ParFactor_IRLS (&fitdata->pX[Nsubset], &fitdata->pY[Nsubset], R[i], D[i], T[i]);
    728     pXmin = MIN (pXmin, fitdata->pX[Nsubset]);
    729     pXmax = MAX (pXmax, fitdata->pX[Nsubset]);
    730     pYmin = MIN (pYmin, fitdata->pY[Nsubset]);
    731     pYmax = MAX (pYmax, fitdata->pY[Nsubset]);
    732 
    733     fitdata->Wx[Nsubset] = 1.0;
    734     fitdata->Wy[Nsubset] = 1.0;   
    735     fitdata->index[Nsubset] = i;
    736     Nsubset++;
    737   }
    738   fitdata->Npts = Nsubset;
    739   float dXRange = pXmax - pXmin;
    740   float dYRange = pYmax - pYmin;
    741   float parRange = hypot (dXRange, dYRange);
    742 
    743   set_variable ("Prange", parRange);
    744   // fprintf (stderr, "par factor range: %f\n", parRange);
    745 
    746   return TRUE;
    747 }
    748 
    749 /* Outlier clipping based on bootstrap-resampling tests of the plx path
    750  * generate Noutlier resampled datasets
    751  * fit the Noutlier plx paths
    752  * determine and save the distribution of dXsig and dYsig for each point
    753  * sort the resulting distributions and find dPsig (median point) for each measurement
    754  * find the 90% point of dPsig : if > dPsigMax, only clip the 10% most deviant points
    755  * set the dPvec values if desired
    756  * -- mask is modified, dPvec values are set
    757  * -- fitdata is unchanged
    758  */
    759 
    760 # define MAX_REJECT 0.1
    761 
    762 int PlxOutlierClip_IRLS (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE) {
    763 
    764   int i, n;
    765 
    766   PlxFit testfit;
    767   testfit.getChisq = FALSE;
    768 
    769   PlxFitData sample;
    770   PlxFitDataAlloc_IRLS (&sample, fitdata->Npts);
    771 
    772   double **dXsig, **dYsig;
    773   ALLOCATE (dXsig, double *, fitdata->Npts);
    774   ALLOCATE (dYsig, double *, fitdata->Npts);
    775   for (i = 0; i < fitdata->Npts; i++) {
    776     ALLOCATE (dXsig[i], double, Noutlier);
    777     ALLOCATE (dYsig[i], double, Noutlier);
    778   }
    779 
    780   int Nsamples = 0;
    781   for (n = 0; n < Noutlier; n++) {
    782     // bootstrap resample (fitdata -> sample)
    783     PlxBootstrapResample_IRLS (fitdata, &sample);
    784      
    785     if (n % 100000 == 99999) fprintf (stderr, ".");
    786 
    787     // fit the sample
    788     if (!FitPMandPar (&testfit,
    789                       sample.X, sample.dX,
    790                       sample.Y, sample.dY, sample.t,
    791                       sample.pX, sample.pY, sample.Npts, VERBOSE)) continue;
    792 
    793     // fprintf (stderr, "%f +/- %f | %f %f\n", testfit.p, testfit.dp, testfit.uR, testfit.uD);
    794 
    795     // find the distances to the path
    796     for (i = 0; i < fitdata->Npts; i++) {
    797       double Xf = testfit.Ro + testfit.uR*fitdata->t[i] + testfit.p*fitdata->pX[i];
    798       double Yf = testfit.Do + testfit.uD*fitdata->t[i] + testfit.p*fitdata->pY[i];
    799       dXsig[i][Nsamples] = fabs(fitdata->X[i] - Xf) / fitdata->dX[i];
    800       dYsig[i][Nsamples] = fabs(fitdata->Y[i] - Yf) / fitdata->dY[i];
    801       // 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]);
    802     }
    803     Nsamples ++;
    804   }
    805 
    806   double *dPsig;
    807   ALLOCATE (dPsig, double, fitdata->Npts);
    808    
    809   for (i = 0; i < fitdata->Npts; i++) {
    810     dsort (dXsig[i], Nsamples);
    811     dsort (dYsig[i], Nsamples);
    812 
    813     // choose the median values
    814     double dXsigMedian, dYsigMedian;
    815     if (Nsamples % 2) {
    816       int Ncenter = Nsamples / 2;
    817       dXsigMedian = dXsig[i][Ncenter];
    818       dYsigMedian = dYsig[i][Ncenter];
    819     } else {
    820       int Ncenter = Nsamples / 2 - 1;
    821       dXsigMedian = 0.5*(dXsig[i][Ncenter] + dXsig[i][Ncenter + 1]);
    822       dYsigMedian = 0.5*(dYsig[i][Ncenter] + dYsig[i][Ncenter + 1]);
    823     }
    824     // XXX replace with hypotenuse?
    825     dPsig[i] = 0.5*(dXsigMedian + dYsigMedian);
    826     // 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]);
    827   }
    828 
    829   // make a copy of dPsig[] and check if > 10% are > dPsigMax
    830   double *dPsigSort;
    831   ALLOCATE (dPsigSort, double, fitdata->Npts);
    832   for (i = 0; i < fitdata->Npts; i++) {
    833     dPsigSort[i] = dPsig[i];
    834   }
    835   dsort (dPsigSort, fitdata->Npts);
    836   int Nmax = (1.0 - MAX_REJECT)*fitdata->Npts;
    837 
    838   int completeClip = TRUE;
    839   if (dPsigSort[Nmax] > dPsigMax) {
    840     if (VERBOSE) fprintf (stderr, "too many outliers: %f at 90\n", dPsigSort[Nmax]);
    841     dPsigMax = dPsigSort[Nmax];
    842     completeClip = FALSE;
    843   }
    844 
    845   for (i = 0; i < fitdata->Npts; i++) {
    846     if (dPsig[i] < dPsigMax) continue;
    847     int n = fitdata->index[i];
    848     // fprintf (stderr, "clip %d: %f : %f\n", i, fitdata->t[i], dPsig[i]);
    849     mask[n] = 0; // mask these points
    850   }
    851 
    852   // only set dPvec if we have completed the clipping?
    853   if (dPvec) {
    854     for (i = 0; i < dPvec->Nelements; i++) {
    855       dPvec->elements.Flt[i] = NAN;
    856     }
    857     for (i = 0; i < fitdata->Npts; i++) {
    858       int n = fitdata->index[i];
    859       dPvec->elements.Flt[n] = dPsig[i];
    860     }
    861   }
    862 
    863   free (dPsig);
    864   free (dPsigSort);
    865  
    866   for (i = 0; i < fitdata->Npts; i++) {
    867     free (dXsig[i]);
    868     free (dYsig[i]);
    869   }
    870   free (dXsig);
    871   free (dYsig);
    872 
    873   return completeClip;
    874 }
    875 
    876 int VectorRobustStats_IRLS (Vector *vector, double *median, double *sigma) {
    877 
    878   // warn if vector->Nelements > 1000? 10000?)
    879   // warn if vector is not float
    880 
    881   // we need to copy the vector to avoid changing the sort order
    882   double *values = NULL;
    883   ALLOCATE (values, double, vector->Nelements);
    884 
    885   int i;
    886   int Npts = 0;
    887   for (i = 0; i < vector->Nelements; i++) {
    888     if (!isfinite(vector->elements.Flt[i])) continue;
    889     values[Npts] = vector->elements.Flt[i];
    890     Npts++;
    891   }
    892 
    893   dsort (values, Npts);
    894 
    895   if (Npts % 2) {
    896     int Ncenter = Npts / 2;
    897     *median = values[Ncenter];
    898   } else {
    899     int Ncenter = Npts / 2 - 1;
    900     *median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
    901   }
    902 
    903   double Slo = VectorFractionInterpolate_IRLS (values, 0.158655, Npts);
    904   double Shi = VectorFractionInterpolate_IRLS (values, 0.841345, Npts);
    905 
    906   *sigma = (Shi - Slo) / 2.0;
    907 
    908   return TRUE;
    909 }
    910 
    911 double VectorFractionInterpolate_IRLS (double *values, float fraction, int Npts) {
    912 
    913   float F = fraction * Npts;
    914   int   N = fraction * Npts;
    915 
    916   if (N < 0        ) return NAN;
    917   if (N >= Npts - 2) return NAN;
    918 
    919   // interpolate between N,N+1
    920    
    921   double S = (F - N) * (values[N+1] - values[N]) + values[N];
    922   return S;
    923 }
    924 
    925 
    926528double Plx_weight_cauchy (double x) {
    927529  double r = x / 2.385;

  • trunk/Ohana/src/opihi/cmd.data/read_vectors.c

    r39217 r39227  
    348348}
    349349
     350# undef ESCAPE
    350351# define ESCAPE(...) {          \
    351352  gprint (GP_ERR, __VA_ARGS__); \

  • trunk/Ohana/src/opihi/cmd.data/reindex.c

    r34461 r39227  
    11# include "data.h"
    2 
    3 # define ESCAPE(MSG,...){ gprint (GP_ERR, MSG, __VA_ARGS__); goto error; }
    42
    53int reindex (int argc, char **argv) {

  • trunk/Ohana/src/opihi/dvo/coordmosaic.c

    r36833 r39227  
    224224}
    225225
     226# undef ESCAPE
    226227# define ESCAPE(MSG) { \
    227228    if (src) fclose (src); \

  • trunk/Ohana/src/opihi/include/astro.h

    r37807 r39227  
    1212void FreeAstro (void);
    1313
     14typedef struct {
     15  double *X;
     16  double *Y;
     17  double *t;
     18  double *pX;
     19  double *pY;
     20  double *dX;
     21  double *dY;
     22  double *Wx;
     23  double *Wy;
     24  int *index;
     25  int Npts;
     26} PlxFitData;
     27
     28typedef struct {
     29  double Ro, dRo;
     30  double Do, dDo;
     31
     32  double uR, duR;
     33  double uD, duD;
     34
     35  double p, dp;
     36
     37  double chisq;
     38  int Nfit;
     39  int getChisq;
     40} PlxFit;
     41
     42int VectorRobustStats (Vector *vector, double *median, double *sigma);
     43double VectorFractionInterpolate (double *values, float fraction, int Npts);
     44
     45int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal);
     46int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean);
     47int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE);
     48
     49int PlxFitDataAlloc (PlxFitData *data, int N);
     50void PlxFitDataFree (PlxFitData *data);
     51int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt);
     52
     53int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE);
     54int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius);
     55int ParFactor (double *pR, double *pD, double RA, double DEC, double Time);
     56
     57/***** */
     58
     59int 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);
     60
     61// ?? int IRLS_converged (PlxFit *fit);
     62int Plx_weighted_LS (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts,
     63                     double **A, double **B, int VERBOSE);
     64double Plx_weight_cauchy (double x);
     65double Plx_dpsi_cauchy (double x);
     66double Plx_MAD(double *in, int N);
     67
    1468# endif

  • trunk/Ohana/src/opihi/include/shell.h

    r37049 r39227  
    44# ifndef SHELL_H
    55# define SHELL_H
     6
     7# define ESCAPE(MSG,...) { gprint (GP_ERR, MSG, __VA_ARGS__); return FALSE; }
    68
    79# define ISVAR(a) (isalnum (a) || (a == ':') || (a == '_'))
Note: See TracChangeset for help on using the changeset viewer.