Index: trunk/Ohana/src/relastro/src/FitPM.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPM.c	(revision 39286)
+++ trunk/Ohana/src/relastro/src/FitPM.c	(revision 39287)
@@ -13,4 +13,7 @@
 
   for (i = 0; i < Npoints; i++) {
+    if (points[i].mask) continue; // respect the mask if set
+    fit->Nfit ++;
+
     /* handle case where dX or dY = 0.0 */
     wx = 1.0 / SQ(points[i].dX);
@@ -64,6 +67,4 @@
   fit->dp  = 0.0;
   
-  fit->Nfit = Npoints;
-
   return (TRUE);
 }
Index: trunk/Ohana/src/relastro/src/FitPosPMfixed_IRLS.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPosPMfixed_IRLS.c	(revision 39287)
+++ trunk/Ohana/src/relastro/src/FitPosPMfixed_IRLS.c	(revision 39287)
@@ -0,0 +1,248 @@
+# include "relastro.h"
+
+// These should probably be tunable:
+# define MAX_ITERATIONS 10
+# define FIT_TOLERANCE 1e-4
+# define FLT_TOLERANCE 1e-6
+# define WEIGHT_THRESHOLD 0.3
+
+int FitPosPMfixed_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
+
+  int i,j;
+
+  int Ndof = 2 * Npoints - data->Nterms;
+  
+  // Convert the measurement errors into initial weights.
+  for (i = 0; i < Npoints; i++) {
+    points[i].Wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX);
+    points[i].Wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY);
+  }
+  
+  // Solve OLS equation  
+  if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
+    myAbort ("handle failures, please!");
+    return(FALSE);
+  }
+
+  // Calculate r vector of residuals and least squares sigma
+  double sigma_ols = 0.0;
+  for (i = 0; i < Npoints; i++) {
+    points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro);
+    points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do);
+    sigma_ols += SQ(points[i].rx) + SQ(points[i].ry);
+  }
+  sigma_ols = sqrt(sigma_ols / (float)Ndof);
+
+  // Save OLS covariance and Beta (solution vector, which is actually also saved in fit)
+  for (i = 0; i < data->Nterms; i++) {
+    for (j = 0; j < data->Nterms; j++) {
+      data->Cov[i][j] = data->A[i][j];
+    }
+    data->Beta[i] = data->B[i][0];
+  }
+
+  // Iteratively reweight and solve
+  double sigma_hat = 0.0; // save for the error model
+  int converged = FALSE;
+  int iterations = 0;
+
+  // modify the weight based on the distance from the previous fit.  try up to MAX_ITERATIONS.
+  // at the end "fit", has the last fit parameters
+  for (iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations ++) {
+    // Save Beta.
+    for (i = 0; i < data->Nterms; i ++) {
+      data->Beta_prev[i] = data->Beta[i];
+    }
+
+    // Assign weights based on the deviation
+    for (i = 0; i < Npoints; i++) {
+      points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
+      points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
+    }    
+
+    // Solve with the new weights
+    if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
+      myAbort ("handle failures, please!");
+      return(FALSE);
+    }
+
+    // store the new Beta.
+    for (i = 0; i < data->Nterms; i++) {
+      data->Beta[i] = data->B[i][0];
+    }
+
+    // calculate the residuals:
+    for (i = 0; i < Npoints; i++) {
+      points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro);
+      points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do);
+      points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY));
+    }
+    sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745;
+    
+    // Check convergence
+    converged = TRUE;
+    for (i = 0; i < data->Nterms; i++) {
+      // if we are within FIT_TOLERANCE as a fractional error or FLT_TOLERANCE as an absolute error, we are good
+      if ((fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) && 
+	  (fabs(data->Beta[i] - data->Beta_prev[i]) > FLT_TOLERANCE)) {
+	converged = FALSE;
+      }
+    }
+  }
+  if (!converged) {
+    myAbort ("raise a warning on non-convergence");
+  }
+
+  // calculate the weight thresholds to mask the bad points:
+  double Sum_Wx = 0.0;
+  double Sum_Wy = 0.0;
+  for (i = 0; i < Npoints; i++) {
+    points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
+    points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
+    
+    Sum_Wx += points[i].Wx;
+    Sum_Wy += points[i].Wy;
+  }
+  double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints);
+  double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints);
+
+  // set a mask (which can be used by the bootstrap resampling analysis)
+  for (i = 0; i < Npoints; i++) {
+    // keep if either is above threshold?
+    // drop if either is below threshold?
+    // points are marked as keep by default
+    if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) {
+      points[i].mask = 1; // keep point if mask == 0
+    }
+  }
+
+  // this section calculates the formal error on the weighted fit using the covariance values
+  // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37
+  if (data->getError) {
+    double ax = 0.0, ay = 0.0;
+    double bx = 0.0, by = 0.0;
+
+    for (i = 0; i < Npoints; i++) {
+      ax += dpsi_cauchy(points[i].rx / points[i].dX);
+      ay += dpsi_cauchy(points[i].ry / points[i].dY);
+
+      bx += SQ(points[i].Wx);
+      by += SQ(points[i].Wy);
+    }
+    ax /= 1.0 * Npoints;  // mean(psi_dot(r))
+    ay /= 1.0 * Npoints; 
+    bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p))
+    by /= 1.0 * (Npoints - data->Nterms);
+  
+    double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax;
+    double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay;
+  
+    double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax;
+    double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay;
+
+    // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
+    double sigma_final_x  = MAX(SQ(sigma_robust_x), (2 * Npoints * SQ(sigma_robust_x) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
+    double sigma_final_y  = MAX(SQ(sigma_robust_y), (2 * Npoints * SQ(sigma_robust_y) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
+
+    fit[0].dRo = sqrt(data->Cov[0][0]);
+    fit[0].duR = sqrt(data->Cov[1][1]);
+    fit[0].dDo = sqrt(data->Cov[2][2]);
+    fit[0].duD = sqrt(data->Cov[3][3]);
+
+    fit[0].dRo *= sigma_final_x;
+    fit[0].duR *= sigma_final_x;
+    fit[0].dDo *= sigma_final_y;
+    fit[0].duD *= sigma_final_y;
+  }
+
+  // (optionally) add up the chi square for the fit, only counting the unmasked points
+  double chisq = 0.0;
+  fit[0].Nfit = 0;
+  for (i = 0; i < Npoints; i++) {
+    if (points[i].mask) continue;
+    fit[0].Nfit ++;
+      
+    if (data->getChisq) {
+      double Xf = fit[0].Ro + fit[0].uR*points[i].T;
+      double Yf = fit[0].Do + fit[0].uD*points[i].T;
+      double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX);
+      double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY);
+      chisq += SQ(points[i].X - Xf) * wx;
+      chisq += SQ(points[i].Y - Yf) * wy;
+    }
+  }
+    
+  // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms)
+  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);
+
+  return (TRUE);
+}
+
+int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
+
+  int i;
+  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
+  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
+
+  for (i = 0; i < Npoints; i++) {
+    wx = points[i].Wx;
+    wy = points[i].Wy;
+
+    Wx += wx;
+    Wy += wy;
+
+    double TWx = points[i].T*wx;
+    double TWy = points[i].T*wy;
+
+    Tx += TWx;
+    Ty += TWy;
+    
+    Tx2 += points[i].T*TWx;
+    Ty2 += points[i].T*TWy;
+    
+    Xs += points[i].X*wx;
+    Ys += points[i].Y*wy;
+
+    XT += points[i].X*TWx;
+    YT += points[i].Y*TWy;
+  }
+
+  // X^T W X
+  data->A[0][0] = Wx;
+  data->A[0][1] = Tx;
+
+  data->A[1][0] = Tx;
+  data->A[1][1] = Tx2;
+  data->A[2][2] = Wy;
+  data->A[2][3] = Ty;
+  data->A[3][2] = Ty;
+  data->A[3][3] = Ty2;
+
+  // X^T W Y
+  data->B[0][0] = Xs;
+  data->B[1][0] = XT;
+  data->B[2][0] = Ys;
+  data->B[3][0] = YT;
+
+  if (!dgaussjordan (data->A, data->B, 4, 1)) {
+# if (DEBUG)
+    if (VERBOSE) fprintf (stderr, "error in fit\n");
+    int j;
+    for (i = 0; i < 4; i++) {
+      for (j = 0; j < 4; j++) {
+	fprintf (stderr, "%e ", data->A[i][j]);
+      }
+      fprintf (stderr, " : %e\n", data->B[i][0]);
+    }
+# endif
+    return FALSE;
+  }
+
+  fit->Ro = data->B[0][0];
+  fit->uR = data->B[1][0];
+  fit->Do = data->B[2][0];
+  fit->uD = data->B[3][0];
+  fit->p  = 0.0;
+
+  return TRUE;
+}
Index: trunk/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/ImageOps.c	(revision 39286)
+++ trunk/Ohana/src/relastro/src/ImageOps.c	(revision 39287)
@@ -1195,4 +1195,6 @@
 }
 
+# define SUPER_VERBOSE 0
+
 /** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ 
 // we only optionally apply the sigma limit: for object averages, this should not be used (should it?)
@@ -1203,7 +1205,7 @@
   float mag;
 
-  if (!finite(measure[0].R) || !finite(measure[0].D)) return FALSE;
-  if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks?
-  if (!finite(measure[0].dM)) return FALSE; //XXX is this necessary for all relastro tasks?
+  if (!finite(measure[0].R) || !finite(measure[0].D)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 1\n"); return FALSE; };
+  if (!finite(measure[0].M)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 2\n"); return FALSE; }; //XXX is this necessary for all relastro tasks?
+  if (!finite(measure[0].dM)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 3\n"); return FALSE; }; //XXX is this necessary for all relastro tasks?
   
   /* select measurements by photcode, or equiv photcode, if specified */
@@ -1214,5 +1216,5 @@
       if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
     }
-    if (!found) return FALSE;
+    if (!found) { if (SUPER_VERBOSE) fprintf (stderr, "filter 4\n"); return FALSE; };
   }
   
@@ -1223,16 +1225,16 @@
       if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
     }
-    if (found) return FALSE;
+    if (found) { if (SUPER_VERBOSE) fprintf (stderr, "filter 5\n"); return FALSE; };
   }  
   
   if ((MinBadQF > 0.0) && isGPC1chip(measure[0].photcode)) {
-    if (!isfinite(measure[0].psfQF)) return FALSE;
-    if (measure[0].psfQF < MinBadQF) return FALSE;
+    if (!isfinite(measure[0].psfQF)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 6\n"); return FALSE; };
+    if (measure[0].psfQF < MinBadQF) { if (SUPER_VERBOSE) fprintf (stderr, "filter 7\n"); return FALSE; };
   }
 
   /* select measurements by time */
   if (TimeSelect) {
-    if (measure[0].t < TSTART) return FALSE;
-    if (measure[0].t > TSTOP) return FALSE;
+    if (measure[0].t < TSTART) { if (SUPER_VERBOSE) fprintf (stderr, "filter 8\n"); return FALSE; };
+    if (measure[0].t > TSTOP) { if (SUPER_VERBOSE) fprintf (stderr, "filter 9\n"); return FALSE; };
   }
   
@@ -1243,13 +1245,13 @@
     } else {
       code = GetPhotcodebyCode (measure[0].photcode);
-      if (!code) return FALSE;
+      if (!code) { if (SUPER_VERBOSE) fprintf (stderr, "filter 10\n"); return FALSE; };
       mask = code[0].astromBadMask;
     }
-    if (mask & measure[0].photFlags) return FALSE;
+    if (mask & measure[0].photFlags) { if (SUPER_VERBOSE) fprintf (stderr, "filter 11\n"); return FALSE; };
   }
 
   /* select measurements by measurement error */
   if (applySigmaLim && (SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) {
-    return FALSE;
+    { if (SUPER_VERBOSE) fprintf (stderr, "filter 12\n"); return FALSE; };
   }
   
@@ -1257,5 +1259,5 @@
   if (ImagSelect) {
     mag = PhotInst (measure, MAG_CLASS_PSF);
-    if (mag < ImagMin || mag > ImagMax) return FALSE;
+    if (mag < ImagMin || mag > ImagMax) { if (SUPER_VERBOSE) fprintf (stderr, "filter 13\n"); return FALSE; };
   }
   
Index: trunk/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 39286)
+++ trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 39287)
@@ -211,4 +211,5 @@
 
     if (average[0].flags & ID_STAR_USE_PAR) {
+      // get the best pm fit and chisq given the set of points (mask is respected)
       if (!FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
 	myAbort ("oops");
Index: trunk/Ohana/src/relastro/src/bcatalog.c
===================================================================
--- trunk/Ohana/src/relastro/src/bcatalog.c	(revision 39286)
+++ trunk/Ohana/src/relastro/src/bcatalog.c	(revision 39287)
@@ -61,7 +61,10 @@
   }
 
+  int myNskip1 = 0, myNskip2 = 0, myNskip3 = 0, myNskip4 = 0, myNskip5 = 0, myNskip6 = 0;
+
   /* exclude stars not in range or with too few measurements */
   for (i = 0; i < catalog[0].Naverage; i++) {
     if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) {
+      myNskip1 ++;
       continue;
     }
@@ -108,4 +111,5 @@
 	REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
       }
+      myNskip2 ++;
       continue;
     }
@@ -146,4 +150,5 @@
 	  Nskip1 ++;
 	}
+	myNskip3 ++;
 	continue;
       }
@@ -159,4 +164,5 @@
 	  Nskip2 ++;
 	}
+	myNskip4 ++;
 	continue;
       }
@@ -180,4 +186,5 @@
 	    free (date);
 	  }
+	  myNskip5 ++;
 	  continue;
 	}
@@ -220,4 +227,5 @@
     if (Nm <= SRC_MEAS_TOOFEW) {
       Nmeasure -= Nm;
+      myNskip6 ++;
       continue;
     }
@@ -230,4 +238,5 @@
     }
   }
+  fprintf (stderr, "skips: %d %d %d %d %d %d\n", myNskip1, myNskip2, myNskip3, myNskip4, myNskip5, myNskip6);
   REALLOCATE (subcatalog[0].average,  Average,     MAX (Naverage, 1));
   REALLOCATE (subcatalog[0].measureT, MeasureTiny, MAX (Nmeasure, 1));
@@ -244,5 +253,5 @@
     LimitDensityCatalog_ByNmeasureGrid (subcatalog, catalog);
   } else {
-    if (VERBOSE2) {
+    if (VERBOSE) {
       char *basename = filebasename (catalog[0].filename);
       fprintf (stderr, "subset of "OFF_T_FMT" ("OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" total) measures for catalog %s\n", 
Index: trunk/Ohana/src/relastro/src/fitpm.c
===================================================================
--- trunk/Ohana/src/relastro/src/fitpm.c	(revision 39286)
+++ trunk/Ohana/src/relastro/src/fitpm.c	(revision 39287)
@@ -14,5 +14,5 @@
 # define N_STARS 3000
 # define N_POINTS 100
-# define N_OUTLIERS 100
+# define N_OUTLIERS 5
 # define N_BOOTSTRAP 100