Index: /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.astro/fitplx.c
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.astro/fitplx.c	(revision 38573)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.astro/fitplx.c	(revision 38574)
@@ -14,4 +14,5 @@
   double *dX;
   double *dY;
+  int *index;
   int Npts;
 } PlxFitData;
@@ -31,4 +32,11 @@
 } 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);
@@ -70,4 +78,11 @@
     remove_argument (N, &argc, argv);
     dPsigMax = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  int Nresample = 0;
+  if ((N = get_argument (argc, argv, "-bootstrap-resample"))) {
+    remove_argument (N, &argc, argv);
+    Nresample = atoi(argv[N]);
     remove_argument (N, &argc, argv);
   }
@@ -112,32 +127,11 @@
   }
 
-  int Ntotal = tvec->Nelements; // XXX check other lengths
-
-  // find mean values to remove
-  double Nmean = 0;
-  double Tmean = 0;
-  double Rmean = 0;
-  double 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;
-
   // Ntotal : all points supplied by user
   // Nsubset : unmasked points
-  // Nmean : unmasked points (a double only used above)
-
-  float Trange = Tmax - Tmin;
-  // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", Rmean, Dmean, Tmean, Trange, Tmin, Tmax);
+  int Ntotal = tvec->Nelements; // XXX check other lengths
+  if (dPvec) ResetVector (dPvec, OPIHI_FLT, Ntotal);
+
+  double Rmean, Dmean, Tmean;
+  PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
 
   /* project coordinates to a plane centered on the object with units of arcsec */
@@ -148,144 +142,27 @@
   coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
 
-  float pXmin = +2.0;
-  float pXmax = -2.0;
-  float pYmin = +2.0;
-  float pYmax = -2.0;
-
-  PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
   PlxFitData fitdata;
   PlxFitDataAlloc (&fitdata, Ntotal);
-
-  // generate the fit values (projected X,Y; parallax factors; 
-
-  // save an index so we can supplied dPsig for the unmasked points
-  int *index;  
-  ALLOCATE (index, int, Ntotal);
-  
-  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 (&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]);
-    index[Nsubset] = i;
-    Nsubset++;
-  }
-  fitdata.Npts = Nsubset;
-  float dXRange = pXmax - pXmin;
-  float dYRange = pYmax - pYmin;
-  float parRange = hypot (dXRange, dYRange);
-	
-  // fprintf (stderr, "par factor range: %f\n", parRange);
+  PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
+
+  PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
 
   // determine dPsig for detections based on Noutlier attempts
   if (Noutlier) {
-    PlxFit testfit;
-    PlxFitData sample;
-    PlxFitDataAlloc (&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);
-    }
-
-    testfit.getChisq = FALSE;
-
-    int n;
-    int Nsamples = 0;
-    for (n = 0; n < Noutlier; n++) {
-      // bootstrap resample (fitdata -> sample)
-      PlxBootstrapResample (&fitdata, &sample);
-      
-      if (n % 1000 == 999) 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]);
-    }
-
-    int Nout = 0;
-    for (i = 0; i < fitdata.Npts; i++) {
-      if (dPsig[i] > dPsigMax) {
-	fprintf (stderr, "clip %d: %f : %f\n", i, fitdata.t[i], dPsig[i]);
-	continue;
-      }
-      sample.X [Nout] = fitdata.X [i];
-      sample.Y [Nout] = fitdata.Y [i];
-      sample.dX[Nout] = fitdata.dX[i];
-      sample.dY[Nout] = fitdata.dY[i];
-      sample.t [Nout] = fitdata.t [i];
-      sample.pX[Nout] = fitdata.pX[i];
-      sample.pY[Nout] = fitdata.pY[i];
-      Nout ++;
-    }
-    sample.Npts = Nout;
-    fprintf (stderr, "keep %d of %d\n", sample.Npts, fitdata.Npts);
-
-    if (dPvec) {
-        ResetVector (dPvec, OPIHI_FLT, Ntotal);
-	for (i = 0; i < Ntotal; i++) {
- 	  dPvec->elements.Flt[i] = NAN;
-	}
-	for (i = 0; i < fitdata.Npts; i++) {
-	  int n = index[i];
- 	  dPvec->elements.Flt[n] = dPsig[i];
-	}
-	dPvec->Nelements = Ntotal;
-    }
-    free (dPsig);
-
-    fitdata = sample;
-  }
-
-  for (i = 0; VERBOSE && (i < fitdata.Npts); i++) {
-    int n = index[i];
-    fprintf (stderr, "%f %f : %f %d : %f %f %f\n", R[n], D[n], T[n], mask[n], fitdata.t[i], fitdata.X[i], fitdata.Y[i]);
+    int clipRetry = TRUE;
+    for (i = 0; clipRetry && (i < 3); i++) {
+      clipRetry = !PlxOutlierClip (&fitdata, mask, Noutlier, dPsigMax, dPvec, VERBOSE);
+
+      // using the new mask values, reset fitdata
+      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);
+    }
+  }
+
+  for (i = 0; (VERBOSE == 2) && (i < fitdata.Npts); i++) {
+    int n = fitdata.index[i];
+    int maskValue = mask ? mask[n] : 0;
+    fprintf (stderr, "%f %f : %f %d : %f %f %f\n", R[n], D[n], T[n], maskValue, fitdata.t[i], fitdata.X[i], fitdata.Y[i]);
   }
 
@@ -298,4 +175,57 @@
   }
 
+  if (Nresample){
+    PlxFitData sample;
+    PlxFitDataAlloc (&sample, fitdata.Npts);
+
+    PlxFit *testfit = NULL;
+    ALLOCATE (testfit, PlxFit, Nresample);
+
+    int Ngood = 0;
+    for (i = 0; i < Nresample; i++) {
+      PlxBootstrapResample (&fitdata, &sample);
+      
+      if (i % 100000 == 99999) fprintf (stderr, ".");
+
+      // fit the sample
+      testfit[Ngood].getChisq = FALSE;
+      if (!FitPMandPar (&testfit[Ngood], 
+			sample.X, sample.dX, 
+			sample.Y, sample.dY, sample.t, 
+			sample.pX, sample.pY, sample.Npts, VERBOSE)) continue;
+      Ngood ++;
+    }
+
+    Vector *pvec, *uRvec, *uDvec, *Rvec, *Dvec;
+
+    // save the Nresample histograms
+    if ((pvec  = SelectVector ("plxVector", ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "plxVector");
+    if ((uRvec = SelectVector ("uRVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "uDVector");
+    if ((uDvec = SelectVector ("uDVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "uRVector");
+    if ((Rvec  = SelectVector ("RoVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "RoVector");
+    if ((Dvec  = SelectVector ("DoVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "DoVector");
+    
+    ResetVector ( pvec, OPIHI_FLT, Ngood);
+    ResetVector (uRvec, OPIHI_FLT, Ngood);
+    ResetVector (uDvec, OPIHI_FLT, Ngood);
+    ResetVector ( Rvec, OPIHI_FLT, Ngood);
+    ResetVector ( Dvec, OPIHI_FLT, Ngood);
+    
+    for (i = 0; i < Ngood; i++) {
+      pvec->elements.Flt[i]  = testfit[i].p;
+      uRvec->elements.Flt[i] = testfit[i].uR;
+      uDvec->elements.Flt[i] = testfit[i].uD;
+      Rvec->elements.Flt[i]  = testfit[i].Ro;
+      Dvec->elements.Flt[i]  = testfit[i].Do;
+    }
+
+    // now calculate median and sigma for each vector
+    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);
+  }
+
   // fprintf (stderr, "%f +/- %f | %f %f\n", fit.p, fit.dp, fit.uR, fit.uD);
 
@@ -326,12 +256,10 @@
   set_variable ("uR",   fit.uR);
   set_variable ("uD",   fit.uD);
-  set_variable ("duR",   fit.duR);
-  set_variable ("duD",   fit.duD);
+  set_variable ("duR",  fit.duR);
+  set_variable ("duD",  fit.duD);
   set_variable ("plx",  fit.p);
   set_variable ("dplx", fit.dp);
   
   set_variable ("Tmean",  Tmean);
-  set_variable ("Trange", Trange);
-  set_variable ("Prange", parRange);
 
   set_variable ("chisq", fit.chisq);
@@ -525,4 +453,5 @@
   ALLOCATE (data->pX, double, N);
   ALLOCATE (data->pY, double, N);
+  ALLOCATE (data->index, int, N);
   return TRUE;
 }
@@ -536,4 +465,5 @@
   FREE (data->pX);
   FREE (data->pY);
+  FREE (data->index);
 }
 
@@ -555,29 +485,245 @@
 }
 
-# if (0)
-int PlxSetMeanEpoch () {
+int PlxSetMeanEpoch (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;
-  double Tmean = 0;
-  double Rmean = 0;
-  double Dmean = 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];
+    *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;
-
-  float Trange = Tmax - Tmin;
-
-}
-# endif
+  *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 (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 (&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->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 (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE) {
+
+  int i, n;
+
+  PlxFit testfit;
+  testfit.getChisq = FALSE;
+
+  PlxFitData sample;
+  PlxFitDataAlloc (&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 (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 (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 (values, 0.158655, Npts);
+  double Shi = VectorFractionInterpolate (values, 0.841345, Npts);
+
+  *sigma = (Shi - Slo) / 2.0;
+
+  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;
+}