Index: trunk/Ohana/src/relastro/src/BootstrapOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/BootstrapOps.c	(revision 39363)
+++ trunk/Ohana/src/relastro/src/BootstrapOps.c	(revision 39364)
@@ -26,5 +26,5 @@
 }
 
-// calculate mean and sigma points for the 5 fit parameter
+// calculate sigma values for the 5 fit parameters
 int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) {
 
Index: trunk/Ohana/src/relastro/src/FitAstromOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitAstromOps.c	(revision 39363)
+++ trunk/Ohana/src/relastro/src/FitAstromOps.c	(revision 39364)
@@ -228,4 +228,5 @@
   fit->chisq = NAN;
   fit->Nfit = 0;
+  fit->converged = FALSE;
   
   return;
Index: trunk/Ohana/src/relastro/src/FitPM.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPM.c	(revision 39363)
+++ trunk/Ohana/src/relastro/src/FitPM.c	(revision 39364)
@@ -1,51 +1,248 @@
 # include "relastro.h"
+int FitPM_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
+int FitPM_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
+
+// These should probably be tunable:
+# define MAX_ITERATIONS 10
+# define FIT_TOLERANCE 1e-4
+# define FLT_TOLERANCE 1e-6
+# define WEIGHT_THRESHOLD 0.3
 
 // initial values of *fit are ignored
-int FitPM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+int FitPM_Basic (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
 
   int i;
 
+  // Convert the measurement errors into initial weights.
+  for (i = 0; i < Npoints; i++) {
+    points[i].Wx = 1.0 / SQ(points[i].dX);
+    points[i].Wy = 1.0 / SQ(points[i].dY);
+  }
+
+  int status = FitPM_MinChisq (fit, data, points, Npoints);
+  if (!status) return FALSE;
+
+  FitPM_SetChisq (fit, data, points, Npoints);
+  return TRUE;
+}
+
+int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  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 = 1.0 / SQ(points[i].dX);
+    points[i].Wy = 1.0 / SQ(points[i].dY);
+  }
+  
+  // Solve OLS equation  
+  if (!FitPM_MinChisq(fit, data, points, Npoints)) {
+    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 (!FitPM_MinChisq(fit, data, points, Npoints)) {
+
+      // restore the last solution and break
+      fit->Ro = data->Beta_prev[0];
+      fit->uR = data->Beta_prev[1];
+      fit->Do = data->Beta_prev[2];
+      fit->uD = data->Beta_prev[3];
+      
+      // 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;
+      break;
+    }
+
+    // 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;
+      }
+    }
+  }
+  fit->converged = converged;
+
+  // 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;
+  }
+
+  FitPM_SetChisq (fit, data, points, Npoints);
+  return (TRUE);
+}
+
+int FitPM_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  myAssert (data->Nterms == 4, "invalid fit arrays");
+
+  int i;
   double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
-
-  myAssert (data->Nterms == 4, "invalid fit arrays");
-
   Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
 
   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);
-    wy = 1.0 / SQ(points[i].dY);
+
+    wx = points[i].Wx;
+    wy = points[i].Wy;
 
     Wx += wx;
     Wy += wy;
 
-    Tx += points[i].T*wx;
-    Ty += points[i].T*wy;
-    
-    Tx2 += SQ(points[i].T)*wx;
-    Ty2 += SQ(points[i].T)*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*points[i].T*wx;
-    YT += points[i].Y*points[i].T*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[0][2] = 0.0;
+  data->A[0][3] = 0.0;
 
   data->A[1][0] = Tx;
   data->A[1][1] = Tx2;
-
+  data->A[1][2] = 0.0;
+  data->A[1][3] = 0.0;
+
+  data->A[2][0] = 0.0;
+  data->A[2][1] = 0.0;
   data->A[2][2] = Wy;
   data->A[2][3] = Ty;
 
+  data->A[3][0] = 0.0;
+  data->A[3][1] = 0.0;
   data->A[3][2] = Ty;
   data->A[3][3] = Ty2;
 
+  // X^T W Y
   data->B[0][0] = Xs;
   data->B[1][0] = XT;
@@ -53,5 +250,7 @@
   data->B[3][0] = YT;
 
-  dgaussjordan (data->A, data->B, 4, 1);
+  if (!dgaussjordan (data->A, data->B, 4, 1)) {
+    return FALSE;
+  }
 
   fit->Ro = data->B[0][0];
@@ -60,5 +259,5 @@
   fit->uD = data->B[3][0];
   fit->p  = 0.0;
-  
+
   fit->dRo = sqrt(data->A[0][0]);
   fit->duR = sqrt(data->A[1][1]);
@@ -66,5 +265,32 @@
   fit->duD = sqrt(data->A[3][3]);
   fit->dp  = 0.0;
-  
-  return (TRUE);
-}
+
+  return TRUE;
+}
+
+int FitPM_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  int i;
+
+  // count unmasked points and (optionally) add up the chi square for the fit
+  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 = 1.0 / SQ(points[i].dX);
+      double wy = 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;
+}
Index: trunk/Ohana/src/relastro/src/FitPM_IRLS.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPM_IRLS.c	(revision 39363)
+++ trunk/Ohana/src/relastro/src/FitPM_IRLS.c	(revision 39364)
@@ -20,6 +20,5 @@
   
   // Solve OLS equation  
-  if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
-    myAbort ("handle failures, please!");
+  if (!FitPM_MinChisq(fit, data, points, Npoints)) {
     return(FALSE);
   }
@@ -62,7 +61,20 @@
 
     // Solve with the new weights
-    if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
-      myAbort ("handle failures, please!");
-      return(FALSE);
+    if (!FitPM_MinChisq(fit, data, points, Npoints)) {
+
+      // restore the last solution and break
+      fit->Ro = data->Beta_prev[0];
+      fit->uR = data->Beta_prev[1];
+      fit->Do = data->Beta_prev[2];
+      fit->uD = data->Beta_prev[3];
+      
+      // 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;
+      break;
     }
 
@@ -90,7 +102,5 @@
     }
   }
-  if (!converged) {
-    myAbort ("raise a warning on non-convergence");
-  }
+  fit->converged = converged;
 
   // calculate the weight thresholds to mask the bad points:
@@ -156,5 +166,5 @@
   }
 
-  // (optionally) add up the chi square for the fit, only counting the unmasked points
+  // count unmasked points and (optionally) add up the chi square for the fit
   double chisq = 0.0;
   fit[0].Nfit = 0;
@@ -179,5 +189,7 @@
 }
 
-int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
+int FitPM_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  myAssert (data->Nterms == 4, "invalid fit arrays");
 
   int i;
@@ -186,4 +198,6 @@
 
   for (i = 0; i < Npoints; i++) {
+    if (points[i].mask) continue; // respect the mask if set
+
     wx = points[i].Wx;
     wy = points[i].Wy;
@@ -211,9 +225,19 @@
   data->A[0][0] = Wx;
   data->A[0][1] = Tx;
+  data->A[0][2] = 0.0;
+  data->A[0][3] = 0.0;
 
   data->A[1][0] = Tx;
   data->A[1][1] = Tx2;
+  data->A[1][2] = 0.0;
+  data->A[1][3] = 0.0;
+
+  data->A[2][0] = 0.0;
+  data->A[2][1] = 0.0;
   data->A[2][2] = Wy;
   data->A[2][3] = Ty;
+
+  data->A[3][0] = 0.0;
+  data->A[3][1] = 0.0;
   data->A[3][2] = Ty;
   data->A[3][3] = Ty2;
@@ -226,14 +250,4 @@
 
   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;
   }
@@ -245,4 +259,11 @@
   fit->p  = 0.0;
 
+  fit->dRo = sqrt(data->A[0][0]);
+  fit->duR = sqrt(data->A[1][1]);
+  fit->dDo = sqrt(data->A[2][2]);
+  fit->duD = sqrt(data->A[3][3]);
+  fit->dp  = 0.0;
+
   return TRUE;
 }
+
Index: trunk/Ohana/src/relastro/src/FitPMandPar.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPMandPar.c	(revision 39363)
+++ trunk/Ohana/src/relastro/src/FitPMandPar.c	(revision 39364)
@@ -1,6 +1,199 @@
 # include "relastro.h"
+int FitPMandPar_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
+int FitPMandPar_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
+
+// These should probably be tunable:
+# define MAX_ITERATIONS 10
+# define FIT_TOLERANCE 1e-4
+# define FLT_TOLERANCE 1e-6
+# define WEIGHT_THRESHOLD 0.3
 
 // initial values of *fit are ignored
-int FitPMandPar (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+int FitPMandPar_Basic (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  int i;
+
+  // Convert the measurement errors into initial weights.
+  for (i = 0; i < Npoints; i++) {
+    points[i].Wx = 1.0 / SQ(points[i].dX);
+    points[i].Wy = 1.0 / SQ(points[i].dY);
+  }
+
+  int status = FitPMandPar_MinChisq (fit, data, points, Npoints);
+  if (!status) return FALSE;
+
+  FitPMandPar_SetChisq (fit, data, points, Npoints);
+  return TRUE;
+}
+
+int FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  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 = 1.0 / SQ(points[i].dX);
+    points[i].Wy = 1.0 / SQ(points[i].dY);
+  }
+  
+  // Solve OLS equation: failure here means the chisq matrix is degenerate, give up entirely
+  if (!FitPMandPar_MinChisq(fit, data, points, Npoints)) {
+    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].pR * fit->p);
+    points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p);
+    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 (!FitPMandPar_MinChisq(fit, data, points, Npoints)) {
+
+      // restore the last solution and break
+      fit->Ro = data->Beta_prev[0];
+      fit->uR = data->Beta_prev[1];
+      fit->Do = data->Beta_prev[2];
+      fit->uD = data->Beta_prev[3];
+      fit->p  = data->Beta_prev[4];
+      
+      // calculate the residuals:
+      for (i = 0; i < Npoints; i++) {
+	points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro + points[i].pR * fit->p);
+	points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p);
+	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;
+      break;
+    }
+
+    // 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].pR * fit->p);
+      points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p);
+      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;
+      }
+    }
+  }
+  fit->converged = converged;
+
+  // 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].dp  = sqrt(data->Cov[4][4]);
+
+    fit[0].dRo *= sigma_final_x;
+    fit[0].duR *= sigma_final_x;
+    fit[0].dDo *= sigma_final_y;
+    fit[0].duD *= sigma_final_y;
+    fit[0].dp  *= sqrt(sigma_final_x * sigma_final_y);
+  }
+
+  FitPMandPar_SetChisq (fit, data, points, Npoints);
+  return (TRUE);
+}
+
+int FitPMandPar_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  myAssert (data->Nterms == 5, "invalid fit arrays");
 
   int i;
@@ -9,6 +202,4 @@
   double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2;
 
-  myAssert (data->Nterms == 5, "invalid fit arrays");
-
   PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0;
   Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
@@ -16,50 +207,62 @@
   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);
-    wy = 1.0 / SQ(points[i].dY);
+
+    wx = points[i].Wx;
+    wy = points[i].Wy;
 
     Wx += wx;
     Wy += wy;
 
-    Tx += points[i].T*wx;
-    Ty += points[i].T*wy;
-    
-    Tx2 += SQ(points[i].T)*wx;
-    Ty2 += SQ(points[i].T)*wy;
+    double TWx = points[i].T*wx;
+    double TWy = points[i].T*wy;
+
+    double XWx = points[i].X*wx;
+    double YWy = points[i].Y*wy;
+
+    Tx += TWx;
+    Ty += TWy;
+    
+    Tx2 += points[i].T*TWx;
+    Ty2 += points[i].T*TWy;
     
     PR += points[i].pR*wx;
     PD += points[i].pD*wy;
     
-    PRT += points[i].pR*points[i].T*wx;
-    PDT += points[i].pD*points[i].T*wy;
-    
-    PRX += points[i].pR*points[i].X*wx;
-    PDY += points[i].pD*points[i].Y*wy;
+    PRT += points[i].pR*TWx;
+    PDT += points[i].pD*TWy;
+    
+    PRX += points[i].pR*XWx;
+    PDY += points[i].pD*YWy;
     
     PR2 += SQ(points[i].pR)*wx;
     PD2 += SQ(points[i].pD)*wy;
 
-    Xs += points[i].X*wx;
-    Ys += points[i].Y*wy;
-
-    XT += points[i].X*points[i].T*wx;
-    YT += points[i].Y*points[i].T*wy;
+    Xs += XWx;
+    Ys += YWy;
+
+    XT += points[i].X*TWx;
+    YT += points[i].Y*TWy;
   }
 
   data->A[0][0] = Wx;
   data->A[0][1] = Tx;
+  data->A[0][2] = 0.0;
+  data->A[0][3] = 0.0;
   data->A[0][4] = PR;
 
   data->A[1][0] = Tx;
   data->A[1][1] = Tx2;
+  data->A[1][2] = 0.0;
+  data->A[1][3] = 0.0;
   data->A[1][4] = PRT;
 
+  data->A[2][0] = 0.0;
+  data->A[2][1] = 0.0;
   data->A[2][2] = Wy;
   data->A[2][3] = Ty;
   data->A[2][4] = PD;
 
+  data->A[3][0] = 0.0;
+  data->A[3][1] = 0.0;
   data->A[3][2] = Ty;
   data->A[3][3] = Ty2;
@@ -78,5 +281,7 @@
   data->B[4][0] = PRX + PDY;
 
-  dgaussjordan (data->A, data->B, 5, 1);
+  if (!dgaussjordan (data->A, data->B, 5, 1)) {
+    return FALSE;
+  }
 
   fit->Ro = data->B[0][0];
@@ -92,4 +297,31 @@
   fit->dp  = sqrt(data->A[4][4]);
   
-  return (TRUE);
+  return TRUE;
 }
+
+int FitPMandPar_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  int i;
+
+  // count unmasked points and (optionally) add up the chi square for the fit
+  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 + fit[0].p*points[i].pR;
+      double Yf = fit[0].Do + fit[0].uD*points[i].T + fit[0].p*points[i].pD;
+      double wx = 1.0 / SQ(points[i].dX);
+      double wy = 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;
+}
Index: trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c	(revision 39363)
+++ trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c	(revision 39364)
@@ -19,7 +19,6 @@
   }
   
-  // Solve OLS equation
+  // Solve OLS equation: failure here means the chisq matrix is degenerate, give up entirely
   if (!weighted_LS_PLX(fit, data, points, Npoints, VERBOSE)) {
-    myAbort ("handle failures, please!");
     return(FALSE);
   }
@@ -63,6 +62,20 @@
     // Solve with the new weights
     if (!weighted_LS_PLX(fit, data, points, Npoints, VERBOSE)) {
-      myAbort ("handle failures, please!");
-      return(FALSE);
+
+      // restore the last solution and break
+      fit->Ro = data->Beta_prev[0];
+      fit->uR = data->Beta_prev[1];
+      fit->Do = data->Beta_prev[2];
+      fit->uD = data->Beta_prev[3];
+      fit->p  = data->Beta_prev[4];
+      
+      // calculate the residuals:
+      for (i = 0; i < Npoints; i++) {
+	points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro + points[i].pR * fit->p);
+	points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p);
+	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;
+      break;
     }
 
@@ -90,7 +103,5 @@
     }
   }
-  if (!converged) {
-    fprintf (stderr, "raise a warning on non-convergence\n");
-  }
+  fit->converged = converged;
 
   // calculate the weight thresholds to mask the bad points:
@@ -158,5 +169,5 @@
   }
 
-  // (optionally) add up the chi square for the fit, only counting the unmasked points
+  // count unmasked points and (optionally) add up the chi square for the fit
   double chisq = 0.0;
   fit[0].Nfit = 0;
@@ -183,4 +194,6 @@
 int weighted_LS_PLX (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
 
+  myAssert (data->Nterms == 5, "invalid fit arrays");
+
   int i;
   double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
@@ -227,14 +240,22 @@
   data->A[0][0] = Wx;
   data->A[0][1] = Tx;
+  data->A[0][2] = 0.0;
+  data->A[0][3] = 0.0;
   data->A[0][4] = PR;
 
   data->A[1][0] = Tx;
   data->A[1][1] = Tx2;
+  data->A[1][2] = 0.0;
+  data->A[1][3] = 0.0;
   data->A[1][4] = PRT;
 
+  data->A[2][0] = 0.0;
+  data->A[2][1] = 0.0;
   data->A[2][2] = Wy;
   data->A[2][3] = Ty;
   data->A[2][4] = PD;
 
+  data->A[3][0] = 0.0;
+  data->A[3][1] = 0.0;
   data->A[3][2] = Ty;
   data->A[3][3] = Ty2;
@@ -254,14 +275,4 @@
 
   if (!dgaussjordan (data->A, data->B, 5, 1)) {
-# if (DEBUG)
-    if (VERBOSE) fprintf (stderr, "error in fit\n");
-    int j;
-    for (i = 0; i < 5; i++) {
-      for (j = 0; j < 5; j++) {
-	fprintf (stderr, "%e ", data->A[i][j]);
-      }
-      fprintf (stderr, " : %e\n", data->B[i][0]);
-    }
-# endif
     return FALSE;
   }
Index: trunk/Ohana/src/relastro/src/FitPosPMfixed.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPosPMfixed.c	(revision 39363)
+++ trunk/Ohana/src/relastro/src/FitPosPMfixed.c	(revision 39364)
@@ -1,18 +1,200 @@
 # include "relastro.h"
-
-int FitPosPMfixed (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+int FitPosPMfixed_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
+int FitPosPMfixed_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
+
+// These should probably be tunable:
+# define MAX_ITERATIONS 10
+# define FIT_TOLERANCE 1e-4
+# define FLT_TOLERANCE 1e-6
+# define WEIGHT_THRESHOLD 0.3
+
+// fit->uR,uD are used in the fit
+int FitPosPMfixed_Basic (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
 
   int i;
 
+  // Convert the measurement errors into initial weights.
+  for (i = 0; i < Npoints; i++) {
+    points[i].Wx = 1.0 / SQ(points[i].dX);
+    points[i].Wy = 1.0 / SQ(points[i].dY);
+  }
+
+  int status = FitPosPMfixed_MinChisq (fit, data, points, Npoints);
+  if (!status) return FALSE;
+
+  FitPosPMfixed_SetChisq (fit, data, points, Npoints);
+  return TRUE;
+}
+
+int FitPosPMfixed_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  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 = 1.0 / SQ(points[i].dX);
+    points[i].Wy = 1.0 / SQ(points[i].dY);
+  }
+  
+  // Solve OLS equation  
+  if (!FitPosPMfixed_MinChisq(fit, data, points, Npoints)) {
+    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 (!FitPosPMfixed_MinChisq(fit, data, points, Npoints)) {
+
+      // restore the last solution and break
+      fit->Ro = data->Beta_prev[0];
+      fit->Do = data->Beta_prev[1];
+      
+      // 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;
+      break;
+    }
+
+    // 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;
+      }
+    }
+  }
+  fit->converged = converged;
+
+  // 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].dDo = sqrt(data->Cov[1][1]);
+
+    fit[0].dRo *= sigma_final_x;
+    fit[0].dDo *= sigma_final_y;
+  }
+
+  FitPosPMfixed_SetChisq (fit, data, points, Npoints);
+  return (TRUE);
+}
+
+int FitPosPMfixed_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  myAssert (data->Nterms == 2, "invalid fit arrays");
+
+  int i;
   double wx, wy, Wx, Wy, Tx, Ty, Xs, Ys;
-
-  myAssert (data->Nterms == 2, "invalid fit arrays");
-
   Wx = Wy = Tx = Ty = Xs = Ys = 0.0;
 
   for (i = 0; i < Npoints; i++) {
-    /* handle case where dX or dY = 0.0 */
-    wx = 1.0 / SQ(points[i].dX);
-    wy = 1.0 / SQ(points[i].dY);
+    if (points[i].mask) continue; // respect the mask if set
+
+    wx = points[i].Wx;
+    wy = points[i].Wy;
 
     Wx += wx;
@@ -26,15 +208,56 @@
   }
 
-  fit->Ro = (Xs - fit->uR*Tx) / Wx;
-  fit->Do = (Ys - fit->uD*Ty) / Wy;
-
+  // X^T W X
+  data->A[0][0] = Wx;
+  data->A[0][1] = 0.0;
+
+  data->A[1][0] = 0.0;
+  data->A[1][1] = Wy;
+
+  // X^T W Y
+  data->B[0][0] = Xs - fit->uR*Tx;
+  data->B[1][0] = Ys - fit->uD*Ty;
+
+  if (!dgaussjordan (data->A, data->B, 2, 1)) {
+    return FALSE;
+  }
+
+  fit->Ro = data->B[0][0];
+  fit->Do = data->B[1][0];
   fit->p  = 0.0;
-  
-  // should be ~ 1.0 / Wx
-  fit->dRo = sqrt(1.0 / Wx);
-  fit->dDo = sqrt(1.0 / Wy);
-  
-  fit->Nfit = Npoints;
-
-  return (TRUE);
-}
+
+  fit->dRo = sqrt(data->A[0][0]);
+  fit->dDo = sqrt(data->A[1][1]);
+  fit->duR = NAN;
+  fit->duD = NAN;
+  fit->dp  = 0.0;
+
+  return TRUE;
+}
+
+int FitPosPMfixed_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  int i;
+
+  // count unmasked points and (optionally) add up the chi square for the fit
+  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 = 1.0 / SQ(points[i].dX);
+      double wy = 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;
+}
Index: trunk/Ohana/src/relastro/src/FitPosPMfixed_IRLS.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPosPMfixed_IRLS.c	(revision 39363)
+++ trunk/Ohana/src/relastro/src/FitPosPMfixed_IRLS.c	(revision 39364)
@@ -20,6 +20,5 @@
   
   // Solve OLS equation  
-  if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
-    myAbort ("handle failures, please!");
+  if (!weighted_LS_POS(fit, data, points, Npoints, VERBOSE)) {
     return(FALSE);
   }
@@ -62,7 +61,18 @@
 
     // Solve with the new weights
-    if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
-      myAbort ("handle failures, please!");
-      return(FALSE);
+    if (!weighted_LS_POS(fit, data, points, Npoints, VERBOSE)) {
+
+      // restore the last solution and break
+      fit->Ro = data->Beta_prev[0];
+      fit->Do = data->Beta_prev[1];
+      
+      // 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;
+      break;
     }
 
@@ -90,7 +100,5 @@
     }
   }
-  if (!converged) {
-    myAbort ("raise a warning on non-convergence");
-  }
+  fit->converged = converged;
 
   // calculate the weight thresholds to mask the bad points:
@@ -146,15 +154,11 @@
 
     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].dDo = sqrt(data->Cov[1][1]);
 
     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
+  }
+
+  // count unmasked points and (optionally) add up the chi square for the fit
   double chisq = 0.0;
   fit[0].Nfit = 0;
@@ -179,9 +183,12 @@
 }
 
-int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
+int weighted_LS_POS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
+
+  myAssert (data->Nterms == 2, "invalid fit arrays");
 
   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;
+  double wx, wy, Wx, Wy, Tx, Ty, Xs, Ys;
+
+  Wx = Wy = Tx = Ty = Xs = Ys = 0.0;
 
   for (i = 0; i < Npoints; i++) {
@@ -192,55 +199,28 @@
     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;
+    Tx += points[i].T*wx;
+    Ty += points[i].T*wy;
     
     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;
+  data->A[0][1] = 0.0;
+
+  data->A[1][0] = 0.0;
+  data->A[1][1] = Wy;
 
   // 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
+  data->B[0][0] = Xs - fit->uR*Tx;
+  data->B[1][0] = Ys - fit->uD*Ty;
+
+  if (!dgaussjordan (data->A, data->B, 2, 1)) {
     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->Do = data->B[1][0];
   fit->p  = 0.0;
 
Index: trunk/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 39363)
+++ trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 39364)
@@ -1,4 +1,10 @@
 # include "relastro.h"
 # define PAR_TOOFEW 5
+
+typedef enum {
+  SELECT_MEAS_HAS_DATA  = 1,
+  SELECT_MEAS_HAS_STACK = 2,
+  SELECT_MEAS_HAS_2MASS = 4,
+} SelectMeasureStatus;
 
 int DumpObjectsWith2MASS (Catalog *catalog, int Ncatalog);
@@ -112,4 +118,5 @@
 
   int mode = FIT_MODE; // start with the globally-defined fit mode
+  // mode may be one of FIT_AVERAGE (position only), FIT_PM_ONLY, FIT_PM_AND_PAR
 
   int XVERB = FALSE;
@@ -135,5 +142,6 @@
   FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
 
-  // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
+  // To judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq.
+  // we start with the parallax fits, and step down to fewer and fewer parameters
 
   // if we have too few good detections for the desired fit, or too limited a baseline,
@@ -141,16 +149,17 @@
 
   // fit the parallax + proper-motion model
-  // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail
-  // to fit PM, we do not attempt PAR.  thus failure to fit PAR falls back to PM-only
   if (mode == FIT_PM_AND_PAR) {
     if (Trange < PM_DT_MIN) {
+      // not enough baseline for proper motion, only set mean position
       mode = FIT_AVERAGE;
       goto justPosition;
     }
     if (parRange < PAR_FACTOR_MIN) {
+      // not enough parallax factor range, skip parallax
       mode = FIT_PM_ONLY;
       goto skipParallax;
     }
     if (fitStats->Npoints <= PAR_TOOFEW) {
+      // not enough data, skip parallax
       mode = FIT_PM_ONLY;
       goto skipParallax;
@@ -159,16 +168,21 @@
     // we are going to use the IRLS analysis to calculate the mean solution and the masking
     // then run N_BOOTSTRAP_SAMPLES to measure the errors
-    fitStats->fitdataPar->getError = (N_BOOTSTRAP_SAMPLES < 3);
-    FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, VERBOSE);
-
-    if (N_BOOTSTRAP_SAMPLES >= 3) {
+    fitStats->fitdataPar->getError = !N_BOOTSTRAP_SAMPLES;
+    if (!FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) {
+      mode = FIT_PM_ONLY;
+      goto skipParallax;
+    }
+
+    if (N_BOOTSTRAP_SAMPLES) {
       fitStats->Nfit = 0;
       int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
       if (Nnomask < 5) {
-	myAbort ("handle this case, please");
+	// if we do not have enough points to fit parallax, we cannot do the bootstrap analysis.
+	mode = FIT_PM_ONLY;
+	goto skipParallax;
       }
       for (k = 0; k < fitStats->NfitAlloc; k++) {
 	BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
-	FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask);
+	if (!FitPMandPar_Basic (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask)) continue;
 	fitStats->Nfit ++;
       }
@@ -190,5 +204,15 @@
     // XXX a hard-wired hack...
     // unless there is a clear problems (below) with the parallax fit, we will use it
-    if ((fabs(fitPar.uR) > 4.0) || (fabs(fitPar.uD) > 4.0)) {
+    int valid = TRUE;
+    valid = valid && isfinite(fitPar.uR);
+    valid = valid && isfinite(fitPar.uD);
+    valid = valid && isfinite(fitPar.p);
+    valid = valid && isfinite(fitPar.duR);
+    valid = valid && isfinite(fitPar.duD);
+    valid = valid && isfinite(fitPar.dp);
+    valid = valid && (fabs(fitPar.uR) < 4.0);
+    valid = valid && (fabs(fitPar.uD) < 4.0);
+    valid = valid && (fabs(fitPar.p) < 2.0);
+    if (!valid) {
       mode = FIT_PM_ONLY;
     } else {
@@ -212,20 +236,24 @@
     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");
+      if (!FitPM_Basic (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
+	average[0].flags |= ID_STAR_BAD_PM;
+	goto justPosition;
       }
     } else {
-      fitStats->fitdataPM->getError = (N_BOOTSTRAP_SAMPLES < 3);
-      FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, VERBOSE);
-
-      if (N_BOOTSTRAP_SAMPLES >= 3) {
+      fitStats->fitdataPM->getError = !N_BOOTSTRAP_SAMPLES;
+      if (!FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
+	mode = FIT_AVERAGE;
+	goto justPosition;
+      }
+      if (N_BOOTSTRAP_SAMPLES) {
 	fitStats->Nfit = 0;
 	int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
 	if (Nnomask < 5) {
-	  myAbort ("handle this case, please");
+	  mode = FIT_AVERAGE;
+	  goto justPosition;
 	}
 	for (k = 0; k < fitStats->NfitAlloc; k++) {
 	  BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
-	  if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
+	  if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
 	  fitStats->Nfit ++;
 	}
@@ -244,7 +272,13 @@
     fitStats->Npm ++;
 
-    // XXX a hard-wired hack...
     // unless there is a clear problems (below) with the proper-motion fit or we have a parallax fit, we will use pm fit
-    if ((fabs(fitPM.uR) > 4.0) || (fabs(fitPM.uD) > 4.0)) {
+    int valid = TRUE;
+    valid = valid && isfinite(fitPM.uR);
+    valid = valid && isfinite(fitPM.uD);
+    valid = valid && isfinite(fitPM.duR);
+    valid = valid && isfinite(fitPM.duD);
+    valid = valid && (fabs(fitPM.uR) < 4.0);
+    valid = valid && (fabs(fitPM.uD) < 4.0);
+    if (!valid) {
       mode = FIT_AVERAGE;
       average[0].flags |= ID_STAR_BAD_PM;
@@ -260,29 +294,42 @@
     // use bootstrap resampling to check the error distribution
     // if we only have one point, this is silly...
-    
-    // XXX I need to rethink here : use a median for the position or weighted average? use IRLS anyway?
-
-    if (fitStats->NfitAlloc == 1) {
+
+    FitAstromResultSetPM (&fitPos, 1, average);
+    if ((average[0].flags & ID_STAR_USE_PAR) || (average[0].flags & ID_STAR_USE_PM)) {
+      if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 
+	goto doneWithFit;
+      }
+    } else {
       FitAstromResultSetPM (&fitPos, 1, average);
-      FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
-    } else {
-      fitStats->Nfit = 0;
-      FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
-      for (k = 0; k < fitStats->NfitAlloc; k++) {
-	BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
-	FitPosPMfixed (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, fitStats->Npoints);
-	fitStats->Nfit ++;
-      }
-      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
-      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
-    }
-    FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
+      if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
+	goto doneWithFit;
+      }
+      if (N_BOOTSTRAP_SAMPLES) {
+	fitStats->Nfit = 0;
+	int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
+	if (Nnomask < 5) {
+	  goto doneWithFit;
+	}
+	FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
+	for (k = 0; k < fitStats->NfitAlloc; k++) {
+	  BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
+	  if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, fitStats->Npoints)) continue;
+	  fitStats->Nfit ++;
+	}
+	BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
+	BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
+      }
+      average[0].flags |= ID_STAR_USE_AVE;
+    }
 
     // project Ro, Do back to RA,DEC
     XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords);
+    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
+
     average[0].flags |= ID_STAR_FIT_AVE;
     fitStats->Nave ++;
   }
 
+ doneWithFit:
   // update the bit flags of which points were used
   for (k = 0; k < fitStats->Npoints; k++) {
@@ -333,22 +380,20 @@
   /* choose the result based on the chisq values */
   // XXXX for now, just use the mode as the result:
-  int result = mode;
   FitAstromResult fit;
   FitAstromResultInit (&fit);
 
-  switch (result) {
-    case FIT_AVERAGE:
-      average[0].flags |= ID_STAR_USE_AVE;
-      fit = fitPos;
-      break;
-    case FIT_PM_ONLY:
-      average[0].flags |= ID_STAR_USE_PM;
-      fit = fitPM;
-      break;
-    case FIT_PM_AND_PAR:
-      average[0].flags |= ID_STAR_USE_PAR;
-      fit = fitPar;
-      break;
-  }
+  if (average[0].flags & ID_STAR_USE_PAR) {
+    myAssert ((average[0].flags & (ID_STAR_USE_PM | ID_STAR_USE_AVE)) == 0, "programming error");
+    fit = fitPar;
+  }
+  if (average[0].flags & ID_STAR_USE_PM) {
+    myAssert ((average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_AVE)) == 0, "programming error");
+    fit = fitPM;
+  }
+  if (average[0].flags & ID_STAR_USE_AVE) {
+    myAssert ((average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM)) == 0, "programming error");
+    fit = fitPos;
+  }
+
   if (XVERB) {
     fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) plx=(%f +/- %f)\n",
@@ -376,5 +421,5 @@
   if (!status) {
     fitStats->Nskip ++;
-    return FALSE; // XXX ??
+    return FALSE;
   }
 
@@ -441,4 +486,6 @@
 int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats) {
 
+  int status;
+
   // set the default values
   average[0].Rstk  = NAN; // RA in degrees
@@ -458,13 +505,24 @@
 
   // select the measurements to be used in this analysis
-  UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, TRUE);
-
-  // too few measurements for average position (require 2 values)
-  if (fitStats->Npoints < 1) return FALSE; // XXX ?? 
+  status = UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, TRUE);
+  if (status & SELECT_MEAS_HAS_STACK) {
+    // if we have any stack measurements, set the default value to the mean position
+    // this is set first, so the position may only be the original stack position
+    // if we exit this function without updating the stack position, it means we 
+    // have a stack measurement, but not a good stack measuremetn
+    average[0].Rstk  = average[0].R; // RA in degrees
+    average[0].Dstk  = average[0].D; // DEC in degrees
+    average[0].dRstk = average[0].dR; // RA scatter in arcsec
+    average[0].dDstk = average[0].dD; // DEC scatter in arcsec
+  }
+
+  // no stack measurements of acceptable quality (defaults to mean position)
+  if (fitStats->Npoints < 1) return FALSE;
   
   double Tmean, Trange, parRange;
   FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
 
-  FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
+  // the stack positions are not statistically independent...
+  FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
   FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
 
@@ -482,6 +540,6 @@
 		      3600*(average[0].D - fitPos.Do));
 
-  // make sure that the fit succeeded
-  int status = TRUE;
+  // check if the fit succeeded
+  status = TRUE;
   status &= finite(fitPos.Ro);
   status &= finite(fitPos.Do);
@@ -522,4 +580,5 @@
 
   int has2MASS = FALSE;
+  int hasStack = FALSE;
 
   int Npoints = fit->Npoints = 0;
@@ -550,4 +609,5 @@
     if (isStack) {
       if (!isGPC1stack(measure[k].photcode)) continue;
+      hasStack = TRUE;
     } else {
       if ( isGPC1stack(measure[k].photcode)) continue;
@@ -616,4 +676,8 @@
     if (isnan(points[Npoints].dY)) continue;
 
+    // avoid dX,dY == 0.0 so we do not have to constantly test for it
+    points[Npoints].dX = MAX(points[Npoints].dX, 0.001);
+    points[Npoints].dY = MAX(points[Npoints].dY, 0.001);
+
     points[Npoints].dT = measure[k].dt;
 
@@ -646,5 +710,9 @@
 
   fit->Npoints = Npoints;
-  return TRUE;
+
+  int status = SELECT_MEAS_HAS_DATA;
+  if (hasStack) status |= SELECT_MEAS_HAS_STACK; 
+  if (has2MASS) status |= SELECT_MEAS_HAS_2MASS; 
+  return status;
 }
 
Index: trunk/Ohana/src/relastro/src/args.c
===================================================================
--- trunk/Ohana/src/relastro/src/args.c	(revision 39363)
+++ trunk/Ohana/src/relastro/src/args.c	(revision 39364)
@@ -121,12 +121,14 @@
     FIT_MODE = FIT_PM_ONLY;
   }
-  if ((N = get_argument (argc, argv, "-par"))) {
-    remove_argument (N, &argc, argv);
-    FIT_MODE = FIT_PAR_ONLY;
-  }
   if ((N = get_argument (argc, argv, "-pmpar"))) {
     remove_argument (N, &argc, argv);
     FIT_MODE = FIT_PM_AND_PAR;
   }
+
+  // NOTE : this is no longer a valid mode
+  // if ((N = get_argument (argc, argv, "-par"))) {
+  //   remove_argument (N, &argc, argv);
+  //   FIT_MODE = FIT_PAR_ONLY;
+  // }
 
   if ((N = get_argument (argc, argv, "-high-speed"))) {
Index: trunk/Ohana/src/relastro/src/fitobj.c
===================================================================
--- trunk/Ohana/src/relastro/src/fitobj.c	(revision 39364)
+++ trunk/Ohana/src/relastro/src/fitobj.c	(revision 39364)
@@ -0,0 +1,194 @@
+# include "relastro.h"
+
+Catalog *mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad);
+
+# define MJD_MIN_MJD 55197   /* 2010/01/01,00:00:00 */
+# define MJD_MAX_MJD 57023   /* 2015/01/01,00:00:00 */
+# define MJD_MIN_YRS 10.00   /* 2010/01/01,00:00:00 */
+# define MJD_MAX_YRS 14.9993 /* 2015/01/01,00:00:00 */
+# define MJD_REF_YRS 12.4148 /* 2012/06/01,00:00:00 */
+
+# define MJD_J2000   51544.0 /* 2010/01/01,00:00:00 */
+
+# define POS_ERROR 0.010 /* arcsec */
+# define OUT_ERROR 0.500 /* arcsec */
+
+// # define N_STARS 3000
+// # define N_POINTS 100
+// # define N_OUTLIERS 5
+// # define N_BOOTSTRAP 100
+
+# define dcos(THETA) cos(RAD_DEG*THETA)
+# define dsin(THETA) sin(RAD_DEG*THETA)
+
+enum {
+  FIT_POS_NONE,
+  FIT_POS_IRLS,
+  FIT_POS_BOOT,
+  FIT_POS_NOCLIP,
+  FIT_PM_IRLS,
+  FIT_PM_BOOT,
+  FIT_PM_NOCLIP,
+  FIT_PLX_IRLS,
+  FIT_PLX_BOOT,
+  FIT_PLX_NOCLIP,
+};
+
+static int Nsecfilt = 0;
+
+int main (int argc, char **argv) {
+  
+  int N_STARS     = 3000;
+  int N_POINTS    =  100;
+  int N_OUTLIERS  =   10;
+  int N_BOOTSTRAP =  100;
+
+  int N;
+  if ((N = get_argument (argc, argv, "-Nstars"))) {
+    remove_argument (N, &argc, argv);
+    N_STARS = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Npoints"))) {
+    remove_argument (N, &argc, argv);
+    N_POINTS = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Noutliers"))) {
+    remove_argument (N, &argc, argv);
+    N_OUTLIERS = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Nbootstrap"))) {
+    remove_argument (N, &argc, argv);
+    N_BOOTSTRAP = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 2) {
+    fprintf (stderr, "USAGE: %s (mode)\n", argv[0]);
+    fprintf (stderr, "  mode: pos, pos-irls, pos-boot, pm, pm-irls, pm-boot, plx, plx-irls, plx-boot\n");
+    fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -Noutliers [10], -Nbootstrap [100]\n");
+    exit (2);
+  }
+
+  // load reference table here
+  if (!LoadPhotcodesText ("dvo.photcodes")) {
+    fprintf (stderr, "failed to load photcodes\n");
+    exit (2);
+  }
+  Nsecfilt = GetPhotcodeNsecfilt();
+
+  // init the random seed
+  { 
+    struct timeval now;
+    gettimeofday (&now, NULL);
+    long A = now.tv_sec + now.tv_usec * 1000000;
+    srand48(A);
+  }
+
+  ohana_gaussdev_init ();
+
+  int i;
+
+  FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, N_BOOTSTRAP);
+
+  int Nstars = N_STARS;
+  for (i = 0; i < Nstars; i++) {
+
+    double Ro  = 360.0*(drand48() - 0.5);
+    double Phi = 2.0*(drand48() - 0.5);
+    double Do  = DEG_RAD * asin(Phi);
+
+    double uR  = 0.5*(drand48() - 0.5);
+    double uD  = 0.5*(drand48() - 0.5);
+
+    double plx = 0.5*(drand48() + 0.0);
+
+    // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points
+    Catalog *catalog = mkstar (fitStats, Ro, Do, uR, uD, plx, N_POINTS, N_OUTLIERS);
+
+    if (!UpdateObjects_Chips (catalog->average, catalog->secfilt, catalog->measureT, NULL, Nsecfilt, fitStats, 0, 0)) goto escape;
+
+    double Tmean = ohana_sec_to_mjd(catalog->average->Tmean);
+    fprintf (stdout, "%12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f |  %8.6f %8.6f %8.6f %8.6f %8.6f  %d | %f %f %f\n",
+	     Ro, Do, uR, uD, plx, 
+	     catalog->average->R,  catalog->average->D,  catalog->average->uR,  catalog->average->uD,  catalog->average->P, Tmean, 
+	     catalog->average->dR, catalog->average->dD, catalog->average->duR, catalog->average->duD, catalog->average->dP, catalog->average->Npos,
+	     catalog->average->ChiSqAve, catalog->average->ChiSqPM, catalog->average->ChiSqPar);
+
+    continue;
+
+  escape:
+    fprintf (stdout, "%12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f |  %8.6f %8.6f %8.6f %8.6f %8.6f  %d | %f\n",
+	     NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, 0.0, NAN, NAN, NAN, NAN, NAN, 0, 0.0);
+  }
+  // return exit_status();
+  ohana_gaussdev_free();
+
+  exit (0);
+}
+
+Catalog *mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad) {
+  
+  int i;
+
+  int Ntotal = Npoints + Nbad;
+
+  ALLOCATE_PTR (catalog, Catalog, 1);
+  ALLOCATE (catalog->average, Average, 1);
+  ALLOCATE (catalog->secfilt, SecFilt, Nsecfilt);
+  ALLOCATE (catalog->measureT, MeasureTiny, Npoints + Nbad);
+  
+  dvo_average_init (catalog->average);
+  for (i = 0; i < Ntotal; i++) {
+    dvo_measureT_init (&catalog->measureT[i]);
+  }
+
+  for (i = 0; i < Npoints; i++) {
+    // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined)
+    double Tmjd = MJD_MIN_MJD + drand48()*(MJD_MAX_MJD - MJD_MIN_MJD);
+
+    catalog->measureT[i].t = ohana_mjd_to_sec (Tmjd);
+    double Tyears = (Tmjd - MJD_J2000) / 365.25;
+    
+    double pR, pD;
+    ParFactor (&pR, &pD, Ro, Do, Tyears);
+
+    double dR = ohana_gaussdev_rnd(0.0, POS_ERROR);
+    double dD = ohana_gaussdev_rnd(0.0, POS_ERROR);
+
+    catalog->measureT[i].R = Ro + (dR + plx*pR + uR*(Tyears - MJD_REF_YRS))/3600/dcos(Do);
+    catalog->measureT[i].D = Do + (dD + plx*pD + uD*(Tyears - MJD_REF_YRS))/3600;
+
+    catalog->measureT[i].dXccd = POS_ERROR;
+    catalog->measureT[i].dYccd = POS_ERROR;
+
+    catalog->measureT[i].photcode = 10233;
+  }
+
+  for (i = Npoints; i < Ntotal; i++) {
+    // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined)
+    double Tmjd = MJD_MIN_MJD + drand48()*(MJD_MAX_MJD - MJD_MIN_MJD);
+
+    catalog->measureT[i].t = ohana_mjd_to_sec (Tmjd);
+    double Tyears = (Tmjd - MJD_J2000) / 365.25;
+    
+    double pR, pD;
+    ParFactor (&pR, &pD, Ro, Do, Tyears);
+
+    double dR = OUT_ERROR*(drand48() - 0.5);
+    double dD = OUT_ERROR*(drand48() - 0.5);
+
+    catalog->measureT[i].R = Ro + (dR + plx*pR + uR*(Tyears - MJD_REF_YRS))/3600/dcos(Do);
+    catalog->measureT[i].D = Do + (dD + plx*pD + uD*(Tyears - MJD_REF_YRS))/3600;
+
+    catalog->measureT[i].dXccd = ToShortPixels(POS_ERROR);
+    catalog->measureT[i].dYccd = ToShortPixels(POS_ERROR);
+
+    catalog->measureT[i].photcode = 10233;
+  }
+  catalog->average->Nmeasure = Ntotal;
+
+  return catalog;
+}
Index: trunk/Ohana/src/relastro/src/fitpm.c
===================================================================
--- trunk/Ohana/src/relastro/src/fitpm.c	(revision 39363)
+++ trunk/Ohana/src/relastro/src/fitpm.c	(revision 39364)
@@ -12,8 +12,8 @@
 # define OUT_ERROR 0.500 /* arcsec */
 
-# define N_STARS 3000
-# define N_POINTS 100
-# define N_OUTLIERS 5
-# define N_BOOTSTRAP 100
+// # define N_STARS 3000
+// # define N_POINTS 100
+// # define N_OUTLIERS 5
+// # define N_BOOTSTRAP 100
 
 # define dcos(THETA) cos(RAD_DEG*THETA)
@@ -21,6 +21,10 @@
 
 enum {
-  FIT_PM_NONE,
+  FIT_POS_NONE,
+  FIT_POS_IRLS,
+  FIT_POS_BOOT,
+  FIT_POS_NOCLIP,
   FIT_PM_IRLS,
+  FIT_PM_BOOT,
   FIT_PM_NOCLIP,
   FIT_PLX_IRLS,
@@ -31,10 +35,48 @@
 int main (int argc, char **argv) {
   
+  int N_STARS     = 3000;
+  int N_POINTS    =  100;
+  int N_OUTLIERS  =   10;
+  int N_BOOTSTRAP =  100;
+
+  int N;
+  if ((N = get_argument (argc, argv, "-Nstars"))) {
+    remove_argument (N, &argc, argv);
+    N_STARS = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Npoints"))) {
+    remove_argument (N, &argc, argv);
+    N_POINTS = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Noutliers"))) {
+    remove_argument (N, &argc, argv);
+    N_OUTLIERS = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Nbootstrap"))) {
+    remove_argument (N, &argc, argv);
+    N_BOOTSTRAP = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   if (argc != 2) {
     fprintf (stderr, "USAGE: %s (mode)\n", argv[0]);
+    fprintf (stderr, "  mode: pos, pos-irls, pos-boot, pm, pm-irls, pm-boot, plx, plx-irls, plx-boot\n");
+    fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -Noutliers [10], -Nbootstrap [100]\n");
     exit (2);
   }
 
-  int mode = FIT_NONE;
+  int mode = FIT_POS_NONE;
+  if (!strcasecmp(argv[1], "pos")) {
+    mode = FIT_POS_NOCLIP;
+  }
+  if (!strcasecmp(argv[1], "pos-irls")) {
+    mode = FIT_POS_IRLS;
+  }
+  if (!strcasecmp(argv[1], "pos-boot")) {
+    mode = FIT_POS_BOOT;
+  }
   if (!strcasecmp(argv[1], "pm")) {
     mode = FIT_PM_NOCLIP;
@@ -42,4 +84,7 @@
   if (!strcasecmp(argv[1], "pm-irls")) {
     mode = FIT_PM_IRLS;
+  }
+  if (!strcasecmp(argv[1], "pm-boot")) {
+    mode = FIT_PM_BOOT;
   }
   if (!strcasecmp(argv[1], "plx")) {
@@ -72,4 +117,6 @@
   fitStats->fitdataPar->getChisq = TRUE;
   fitStats->fitdataPar->getError = TRUE;
+  fitStats->fitdataPos->getChisq = TRUE;
+  fitStats->fitdataPos->getError = TRUE;
   // NOTE: surprisingly, the chisq and error calculations only add ~2-3% to the exec time
 
@@ -88,7 +135,11 @@
     double plx;
     switch (mode) {
+      case FIT_POS_IRLS:
+      case FIT_POS_BOOT:
+      case FIT_POS_NOCLIP:
       case FIT_PM_IRLS:
+      case FIT_PM_BOOT:
       case FIT_PM_NOCLIP:
-	plx = 0.0*(drand48() + 0.0);
+	plx = 0.0;
 	break;
       case FIT_PLX_IRLS:
@@ -101,5 +152,5 @@
     }
 
-
+    // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points
     mkstar (fitStats, Ro, Do, uR, uD, plx, N_POINTS, N_OUTLIERS);
 
@@ -107,36 +158,81 @@
     FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
 
-    FitAstromResult fitPM;    
-    FitAstromResultInit (&fitPM);
+    FitAstromResult fitResult;    
+    FitAstromResultInit (&fitResult);
+
+    int Nnomask;
 
     switch (mode) {
-      case FIT_PM_NOCLIP:
-	FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
-	break;
-      case FIT_PM_IRLS:
-	FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, 0);
-	break;
-      case FIT_PLX_NOCLIP:
-	FitPMandPar (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints);
-	break;
-      case FIT_PLX_IRLS:
-	FitPMandPar_IRLS (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, 0);
-	break;
-
-      case FIT_PLX_BOOT:
-	FitPMandPar_IRLS (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, 0);
+      case FIT_POS_NOCLIP:
+	fitResult.uR = uR;
+	fitResult.uD = uD;
+	if (!FitPosPMfixed_Basic (&fitResult, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) goto escape;
+	break;
+      case FIT_POS_IRLS:
+	fitResult.uR = uR;
+	fitResult.uD = uD;
+	if (!FitPosPMfixed_IRLS (&fitResult, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) goto escape;
+	break;
+      case FIT_POS_BOOT:
+	fitResult.uR = uR;
+	fitResult.uD = uD;
+	if (!FitPosPMfixed_IRLS (&fitResult, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) goto escape;
 	fitStats->Nfit = 0;
-	int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
+	Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
 	for (k = 0; k < fitStats->NfitAlloc; k++) {
 	  BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
-	  FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask);
+	  fitStats->fit[k].uR = uR;
+	  fitStats->fit[k].uD = uD;
+	  if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue;
 	  fitStats->Nfit ++;
 	}
 	// these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
-	BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
-	BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
-	BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
-	BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
-	BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
+	BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
+	BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
+	break;
+
+      case FIT_PM_NOCLIP:
+	if (!FitPM_Basic (&fitResult, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) goto escape;
+	break;
+      case FIT_PM_IRLS:
+	if (!FitPM_IRLS (&fitResult, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) goto escape;
+	break;
+      case FIT_PM_BOOT:
+	if (!FitPM_IRLS (&fitResult, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) goto escape;
+	fitStats->Nfit = 0;
+	Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
+	for (k = 0; k < fitStats->NfitAlloc; k++) {
+	  BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
+	  if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, Nnomask)) continue;
+	  fitStats->Nfit ++;
+	}
+	// these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
+	BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
+	BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
+	BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
+	BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
+	break;
+
+      case FIT_PLX_NOCLIP:
+	if (!FitPMandPar_Basic (&fitResult, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) goto escape;
+	break;
+      case FIT_PLX_IRLS:
+	if (!FitPMandPar_IRLS (&fitResult, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) goto escape;
+	break;
+      case FIT_PLX_BOOT:
+	if (!FitPMandPar_IRLS (&fitResult, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) goto escape;
+	fitStats->Nfit = 0;
+	Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
+	for (k = 0; k < fitStats->NfitAlloc; k++) {
+	  BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
+	  if (!FitPMandPar_Basic (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask)) continue;
+	  fitStats->Nfit ++;
+	}
+	// these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
+	BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
+	BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
+	BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
+	BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
+	BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
 	break;
 
@@ -145,11 +241,15 @@
     }
 
-    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords);
+    XY_to_RD (&fitResult.Ro, &fitResult.Do, fitResult.Ro, fitResult.Do, &fitStats->coords);
 
     fprintf (stdout, "%12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f |  %8.6f %8.6f %8.6f %8.6f %8.6f  %d | %f\n",
 	     Ro, Do, uR, uD, plx, 
-	     fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p, Tmean, fitPM.dRo, fitPM.dDo, fitPM.duR, fitPM.duD, fitPM.dp, fitPM.Nfit, fitPM.chisq); 
-
-    // ok (TRUE, "fitted star");
+	     fitResult.Ro, fitResult.Do, fitResult.uR, fitResult.uD, fitResult.p, Tmean, fitResult.dRo, fitResult.dDo, fitResult.duR, fitResult.duD, fitResult.dp, fitResult.Nfit, fitResult.chisq); 
+
+    continue;
+
+  escape:
+    fprintf (stdout, "%12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f |  %8.6f %8.6f %8.6f %8.6f %8.6f  %d | %f\n",
+	     NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, 0.0, NAN, NAN, NAN, NAN, NAN, 0, 0.0);
   }
   // return exit_status();