Changeset 39227 for trunk/Ohana/src/opihi/cmd.astro/fitplx.c

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

File:

• trunk/Ohana/src/opihi/cmd.astro/fitplx.c (modified) (4 diffs)

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

@@ -1 +1 @@
 # 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 +355 @@
   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 +370 @@
     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 +432 @@
     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++;