Index: trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c	(revision 39235)
+++ trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c	(revision 39236)
@@ -25,14 +25,4 @@
   }
 
-  int Noutlier = 0;
-  float dPsigMax = FLT_MAX;
-  if ((N = get_argument (argc, argv, "-outlier-tests"))) {
-    remove_argument (N, &argc, argv);
-    Noutlier = atoi(argv[N]);
-    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"))) {
@@ -42,12 +32,4 @@
   }
 
-  Vector *dPvec = NULL;
-  if ((N = get_argument (argc, argv, "-dPsig"))) {
-    if (!Noutlier) { gprint (GP_ERR, "-dPsig requires -outlier-tests to be non-zero\n"); return FALSE; }
-    remove_argument (N, &argc, argv);
-    if (!(dPvec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
-    remove_argument (N, &argc, argv);
-  }
-
   int max_iterations = 10;
   if ((N = get_argument (argc, argv, "-max-iterations"))) {
@@ -57,11 +39,4 @@
   }
 
-  Vector *outMask = NULL;
-  if ((N = get_argument (argc, argv, "-out-mask"))) {
-    remove_argument (N, &argc, argv);
-    if (!(outMask = SelectVector(argv[N], ANYVECTOR, TRUE))) return FALSE;
-    remove_argument (N, &argc, argv);
-  }
-
   double outlier_limit = 0.1;
   if ((N = get_argument (argc, argv, "-outlier-limit"))) {
@@ -74,8 +49,5 @@
   if (argc != 6) {
     gprint (GP_ERR, "USAGE: fitplx (ra) (dR) (dec) (dD) (mjd) [-mask mask] [-v] [-vv]\n");
-    gprint (GP_ERR, "  -outlier-tests Nsamples dPsigMax : run Nsample bootstrap-resamples to define the path deviations and reject based on dPsigMax\n");
-    gprint (GP_ERR, "  -dPsig vec : save path deviations in vec\n");
     gprint (GP_ERR, "  -max-iterations : maximum number of IRLS iterations to run (default 10)\n");
-    gprint (GP_ERR, "  -outMask vect : save outlier mask in vec\n");
     gprint (GP_ERR, "  -outlier-limit : fraction of average weight to reject on (default 0.1)\n");
     return (FALSE);
@@ -109,6 +81,4 @@
   // Nsubset : unmasked points
   int Ntotal = tvec->Nelements; // XXX check other lengths
-  if (dPvec) ResetVector (dPvec, OPIHI_FLT, Ntotal);
-  if (outMask) ResetVector (outMask, OPIHI_INT, Ntotal);
   
   double Rmean, Dmean, Tmean;
@@ -125,26 +95,11 @@
   PlxFitDataAlloc (&fitdata, Ntotal);
   PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
+  // fitdata only contains the points which are not masked
 
   PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
 
-  if (0) {
-    // This should be accomplished via the IRLS call
-    // determine dPsig for detections based on Noutlier attempts
-    if (Noutlier) {
-      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;
+    int maskValue = mask ? mask[n] : 1;
     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]);
   }
@@ -159,33 +114,33 @@
     return FALSE;
   }
-
-  // Set the output mask based on the input mask and the outlier limits.
-  if (outMask) {
-    double Sum_Wx = 0;
-    double Sum_Wy = 0;
-    int *outMask_V = outMask->elements.Int;
-    for (i = 0; i < fitdata.Npts; i++) {
-      Sum_Wx += fitdata.Wx[i];
-      Sum_Wy += fitdata.Wy[i];
-    }
-    for (i = 0; i < Ntotal; i++) {
-      outMask_V[i] = mask ? mask[i] : 0;
-    }
-    
-    for (i = 0; i < fitdata.Npts; i++) {
-      if ((fitdata.Wx[i] > outlier_limit * Sum_Wx / (1.0 * fitdata.Npts))||
-	  (fitdata.Wy[i] > outlier_limit * Sum_Wy / (1.0 * fitdata.Npts))) {
-	int n = fitdata.index[i];
-	outMask_V[n] = 1;
-
-	if (VERBOSE == 2) {
-	  fprintf (stderr, "%f %f : %f %d : %f %f %f : %f %f %f %f\n", R[n], D[n], T[n], outMask_V[n], fitdata.t[i], fitdata.X[i], fitdata.Y[i], fitdata.Wx[i], fitdata.Wy[i], Sum_Wx, Sum_Wy);
-	}
-      }
-    }
-  }
-
-  
+  // FitPMandPar_IRLS sets the values of Wx,Wy based on the fit distance
+
+  // update the mask based on the input mask and the outlier limits.
+  double Sum_Wx = 0;
+  double Sum_Wy = 0;
+  mask = mvec->elements.Int;
+  for (i = 0; i < fitdata.Npts; i++) {
+    Sum_Wx += fitdata.Wx[i];
+    Sum_Wy += fitdata.Wy[i];
+  }
+  for (i = 0; i < fitdata.Npts; i++) {
+    // fitdata only includes the previously unmasked points
+    if ((fitdata.Wx[i] < outlier_limit * Sum_Wx / (1.0 * fitdata.Npts))||
+	(fitdata.Wy[i] < outlier_limit * Sum_Wy / (1.0 * fitdata.Npts))) {
+      int n = fitdata.index[i];
+      mask[n] = 0;
+      
+      if (VERBOSE == 2) {
+	  fprintf (stderr, "%f %f : %f %d : %f %f %f : %f %f %f %f\n", R[n], D[n], T[n], mask[n], fitdata.t[i], fitdata.X[i], fitdata.Y[i], fitdata.Wx[i], fitdata.Wy[i], Sum_Wx, Sum_Wy);
+      }
+    }
+  }
+
   if (Nresample){
+    // now that the mask has been updated, we need to recalculate mean epoch and positions
+    // XXX make this conditional on actually masking unmasked points above
+    PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
+    PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
+
     PlxFitData sample;
     PlxFitDataAlloc (&sample, fitdata.Npts);
@@ -319,4 +274,5 @@
   }
 
+  // XXX EAM : I suggest making these static as well like the A and B matrices
   // Allocate the IRLS specific objects
   ALLOCATE (Cov, double *, 5);
@@ -327,6 +283,4 @@
   ALLOCATE(Beta, double, 5);
   ALLOCATE(Beta_prev, double, 5);
-  //  ALLOCATE(Wx, double, Npts);
-  //  ALLOCATE(Wy, double, Npts);
   ALLOCATE(rx,  double, Npts);
   ALLOCATE(ry,  double, Npts);
@@ -348,9 +302,13 @@
   sigma_ols = 0.0;
   for (i = 0; i < Npts; i++) {
+    // double Xf = B[0][0] + B[1][0]*T[i] + B[4][0]*pR[i];
+    // double Yf = B[2][0] + B[3][0]*T[i] + B[4][0]*pD[i];
+
     rx[i] = X[i] - (T[i] * B[1][0] + B[0][0] + B[4][0] * pR[i]);
     ry[i] = Y[i] - (T[i] * B[3][0] + B[2][0] + B[4][0] * pD[i]);
+
+
     //    u[i] = r[i] /
     sigma_ols += SQ(rx[i]) + SQ(ry[i]);
-
   }
   sigma_ols = sqrt(sigma_ols / (Npts - 5.0));
@@ -390,5 +348,5 @@
 
     // Solve
-    if (!weighted_LS_PLX(T,pD,pR,X,Wx,Y,Wy,Npts,A,B,VERBOSE)) {
+    if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,A,B,VERBOSE)) {
       // Handle fail case
       return(FALSE);
@@ -478,6 +436,4 @@
     }
   }
-  
-  
 
   // Finish