Index: trunk/Ohana/src/relastro/Makefile
===================================================================
--- trunk/Ohana/src/relastro/Makefile	(revision 39239)
+++ trunk/Ohana/src/relastro/Makefile	(revision 39241)
@@ -170,4 +170,6 @@
 $(SRC)/FitPM.$(ARCH).o               \
 $(SRC)/FitPM_IRLS.$(ARCH).o               \
+$(SRC)/FitPMandPar.$(ARCH).o               \
+$(SRC)/FitPMandPar_IRLS.$(ARCH).o               \
 $(SRC)/mkpolyterm.$(ARCH).o            \
 $(SRC)/fitpoly.$(ARCH).o
Index: trunk/Ohana/src/relastro/include/relastro.h
===================================================================
--- trunk/Ohana/src/relastro/include/relastro.h	(revision 39239)
+++ trunk/Ohana/src/relastro/include/relastro.h	(revision 39241)
@@ -750,8 +750,10 @@
 
 int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
+int FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
 
 double MedianAbsDeviation(FitAstromPoint *points, int Npoints);
 
 int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
+int weighted_LS_PLX (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
 
 double weight_cauchy (double x);
Index: trunk/Ohana/src/relastro/src/FitAstromOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitAstromOps.c	(revision 39239)
+++ trunk/Ohana/src/relastro/src/FitAstromOps.c	(revision 39241)
@@ -225,2 +225,48 @@
   return;
 }
+
+double weight_cauchy (double x) {
+  double r = x / 2.385;
+  return (1.0 / (1.0 + SQ(r)));
+}
+
+// dpsi = (d/dx) (x * weight(x))
+double dpsi_cauchy (double x) {
+  double r2 = SQ(x / 2.385);
+  return ((1.0 - r2) / (SQ(1 + r2)));
+}
+
+
+// median absolute deviation
+// MAD = median(abs(x - median(x)))
+double MedianAbsDeviation(FitAstromPoint *points, int Npoints) {
+
+  double *x;
+  double median = 0.0;
+  int i;
+  
+  ALLOCATE(x, double, Npoints);
+  for (i = 0; i < Npoints; i++) {
+    x[i] = points[i].u;
+  }
+  dsort(x, Npoints);
+
+  if (Npoints % 2) {
+    median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]);
+  } else {
+    median = x[(int)(0.5*Npoints)];
+  }
+
+  for (i = 0; i < Npoints; i++ ) {
+    x[i] = fabs(x[i] - median);
+  }
+  dsort(x, Npoints);
+
+  if (Npoints % 2) {
+    median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]);
+  } else {
+    median = x[(int)(0.5*Npoints)];
+  }
+
+  return median;
+}
Index: trunk/Ohana/src/relastro/src/FitPM_IRLS.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPM_IRLS.c	(revision 39239)
+++ trunk/Ohana/src/relastro/src/FitPM_IRLS.c	(revision 39241)
@@ -4,14 +4,14 @@
 # define MAX_ITERATIONS 10
 # define FIT_TOLERANCE 1e-4
+# define FLT_TOLERANCE 1e-6
 # define WEIGHT_THRESHOLD 0.3
 
-/* do we want an init function which does the alloc and a clear function to free? */
 int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
 
   int i,j;
 
-  int dof = 2 * Npoints - 4;
   int p   = 4;
   int n   = 2 * Npoints;
+  int dof = n - p;
   
   // data->A,B,Cov,Beta,Beta_prev are allocated outside by FitAstromDataInit()
@@ -20,6 +20,8 @@
   // Convert the measurement errors into initial weights.
   for (i = 0; i < Npoints; i++) {
-    points[i].Wx = 1 / points[i].dX;
-    points[i].Wy = 1 / points[i].dY;
+    // points[i].Wx = 1 / points[i].dX;
+    // points[i].Wy = 1 / points[i].dY;
+    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);
   }
   
@@ -40,6 +42,6 @@
 
   // Save OLS covariance and Beta (solution vector, which is actually also saved in fit)
-  for (i = 0; i < 4; i++) {
-    for (j = 0; j < 4; j++) {
+  for (i = 0; i < data->Nterms; i++) {
+    for (j = 0; j < data->Nterms; j++) {
       data->Cov[i][j] = data->A[i][j];
     }
@@ -47,5 +49,5 @@
   }
 
-  // Iterately reweight and solve
+  // Iteratively reweight and solve
   double sigma_hat = 0.0; // save for the error model
   int converged = FALSE;
@@ -56,5 +58,5 @@
   for (iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations ++) {
     // Save Beta.
-    for (i = 0; i < 4; i ++) {
+    for (i = 0; i < data->Nterms; i ++) {
       data->Beta_prev[i] = data->Beta[i];
     }
@@ -73,5 +75,5 @@
 
     // store the new Beta.
-    for (i = 0; i < 4; i++) {
+    for (i = 0; i < data->Nterms; i++) {
       data->Beta[i] = data->B[i][0];
     }
@@ -87,6 +89,8 @@
     // Check convergence
     converged = TRUE;
-    for (i = 0; i < 4; i++) {
-      if (fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) {
+    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;
       }
@@ -134,9 +138,8 @@
     fit[0].duD = sqrt(data->Cov[3][3]);
 
-    fit[9].dRo *= sigma_final_x;
-    fit[9].duR *= sigma_final_x;
-
-    fit[9].dDo *= sigma_final_y;
-    fit[9].duD *= sigma_final_y;
+    fit[0].dRo *= sigma_final_x;
+    fit[0].duR *= sigma_final_x;
+    fit[0].dDo *= sigma_final_y;
+    fit[0].duD *= sigma_final_y;
   }
 
@@ -154,19 +157,22 @@
   }
 
-  // add up the chi square for the fit, only counting the unmasked points
+  // (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;
-
+      
     double Xf = fit[0].Ro + fit[0].uR*points[i].T;
     double Yf = fit[0].Do + fit[0].uD*points[i].T;
-    chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
-    chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
+    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;
     fit[0].Nfit ++;
   }
-
+    
   // the reduced chisq is divided by (Ndof = 2*Nfit - 4)
   fit[0].chisq = chisq / (2.0*fit[0].Nfit - 4.0);
+
   return (TRUE);
 }
@@ -175,13 +181,16 @@
 
   int i;
-  double Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
+  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;
-
-    double TWx = points[i].T*points[i].Wx;
-    double TWy = points[i].T*points[i].Wy;
+    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;
@@ -191,6 +200,6 @@
     Ty2 += points[i].T*TWy;
     
-    Xs += points[i].X*points[i].Wx;
-    Ys += points[i].Y*points[i].Wy;
+    Xs += points[i].X*wx;
+    Ys += points[i].Y*wy;
 
     XT += points[i].X*TWx;
@@ -237,48 +246,2 @@
   return TRUE;
 }
-
-double weight_cauchy (double x) {
-  double r = x / 2.385;
-  return (1.0 / (1.0 + SQ(r)));
-}
-
-// dpsi = (d/dx) (x * weight(x))
-double dpsi_cauchy (double x) {
-  double r2 = SQ(x / 2.385);
-  return ((1.0 - r2) / (SQ(1 + r2)));
-}
-
-
-// median absolute deviation
-// MAD = median(abs(x - median(x)))
-double MedianAbsDeviation(FitAstromPoint *points, int Npoints) {
-
-  double *x;
-  double median = 0.0;
-  int i;
-  
-  ALLOCATE(x, double, Npoints);
-  for (i = 0; i < Npoints; i++) {
-    x[i] = points[i].u;
-  }
-  dsort(x, Npoints);
-
-  if (Npoints % 2) {
-    median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]);
-  } else {
-    median = x[(int)(0.5*Npoints)];
-  }
-
-  for (i = 0; i < Npoints; i++ ) {
-    x[i] = fabs(x[i] - median);
-  }
-  dsort(x, Npoints);
-
-  if (Npoints % 2) {
-    median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]);
-  } else {
-    median = x[(int)(0.5*Npoints)];
-  }
-
-  return median;
-}
Index: trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c	(revision 39241)
+++ trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c	(revision 39241)
@@ -0,0 +1,276 @@
+# 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 FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
+
+  int i,j;
+
+  int p   = 5;
+  int n   = 2 * Npoints;
+  int dof = n - p;
+  
+  // data->A,B,Cov,Beta,Beta_prev are allocated outside by FitAstromDataInit()
+  // points->Wx,Wy,rx,ry,u are elements of FitAstromPoint
+  
+  // Convert the measurement errors into initial weights.
+  for (i = 0; i < Npoints; i++) {
+    // points[i].Wx = 1 / points[i].dX;
+    // points[i].Wy = 1 / points[i].dY;
+    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_PLX(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].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 / dof);
+  // sigma_ols = sqrt(sigma_ols / (Npts - 5.0));
+  
+  // 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_PLX(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].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;
+      }
+    }
+  }
+  if (!converged) {
+    myAbort ("raise a warning on non-convergence");
+  }
+
+  // this section calculates the formal error on the weighted fit using the covariance values
+  double Sum_Wx = 0.0;
+  double Sum_Wy = 0.0;
+  if (1) {
+    double ax = 0.0, ay = 0.0;
+    double bx = 0.0, by = 0.0;
+    double lambda = 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);
+    
+      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);
+
+      Sum_Wx += points[i].Wx;
+      Sum_Wy += points[i].Wy;
+    }
+    ax /= 1.0 * Npoints;  // mean(psi_dot(r))
+    ay /= 1.0 * Npoints; 
+    bx /= 1.0 * (Npoints - p); // mean(psi^2(r)) * (N / (N-p))
+    by /= 1.0 * (Npoints - p);
+  
+    double sigma_robust_x = lambda * sqrt(bx) * sigma_hat * 2.385 / ax;
+    double sigma_robust_y = lambda * 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), (n * SQ(sigma_robust_x) + SQ(p * sigma_ols)) / (n + SQ(p)));
+    double sigma_final_y  = MAX(SQ(sigma_robust_y), (n * SQ(sigma_robust_y) + SQ(p * sigma_ols)) / (n + SQ(p)));
+
+    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);
+  }
+
+  // set a mask (which can be used by the bootstrap resampling analysis)
+  double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints);
+  double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints);
+
+  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;
+    }
+  }
+
+  // (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;
+      
+    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 = (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;
+    fit[0].Nfit ++;
+  }
+    
+  // the reduced chisq is divided by (Ndof = 2*Nfit - 4)
+  fit[0].chisq = chisq / (2.0*fit[0].Nfit - 5.0);
+  
+  return (TRUE);
+}
+
+int weighted_LS_PLX (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;
+  double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2;
+
+  PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0;
+  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
+
+  for (i = 0; i < Npoints; i++) {
+
+    // if (VERBOSE == 2) fprintf (stderr, "%f %f : %f %f : %f : %f %f\n", X[i], WX[i], Y[i], WY[i], T[i], pR[i], pD[i]);
+
+    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;
+    
+    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;
+    
+    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;
+  }
+
+  data->A[0][0] = Wx;
+  data->A[0][1] = Tx;
+  data->A[0][4] = PR;
+
+  data->A[1][0] = Tx;
+  data->A[1][1] = Tx2;
+  data->A[1][4] = PRT;
+
+  data->A[2][2] = Wy;
+  data->A[2][3] = Ty;
+  data->A[2][4] = PD;
+
+  data->A[3][2] = Ty;
+  data->A[3][3] = Ty2;
+  data->A[3][4] = PDT;
+
+  data->A[4][0] = PR;
+  data->A[4][1] = PRT;
+  data->A[4][2] = PD;
+  data->A[4][3] = PDT;
+  data->A[4][4] = PR2 + PD2;
+
+  data->B[0][0] = Xs;
+  data->B[1][0] = XT;
+  data->B[2][0] = Ys;
+  data->B[3][0] = YT;
+  data->B[4][0] = PRX + PDY;
+
+  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;
+  }
+
+  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  = data->B[4][0];
+
+  return TRUE;
+}
Index: trunk/Ohana/src/relastro/src/fitpm.c
===================================================================
--- trunk/Ohana/src/relastro/src/fitpm.c	(revision 39239)
+++ trunk/Ohana/src/relastro/src/fitpm.c	(revision 39241)
@@ -1,7 +1,4 @@
 # include "relastro.h"
 
-double rnd_gauss (double mean, double sigma);
-void gauss_init (int Nbin);
-double gaussian (double x, double mean, double sigma);
 int mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad);
 
@@ -15,4 +12,5 @@
 # define OUT_ERROR 0.500 /* arcsec */
 
+# define N_STARS 20000
 # define N_POINTS 100
 # define N_OUTLIERS 10
@@ -21,6 +19,33 @@
 # define dsin(THETA) sin(RAD_DEG*THETA)
 
+enum {
+  FIT_PM_NONE,
+  FIT_PM_IRLS,
+  FIT_PM_NOCLIP,
+  FIT_PLX_IRLS,
+  FIT_PLX_NOCLIP,
+};
+
 int main (int argc, char **argv) {
   
+  if (argc != 2) {
+    fprintf (stderr, "USAGE: %s (mode)\n", argv[0]);
+    exit (2);
+  }
+
+  int mode = FIT_NONE;
+  if (!strcasecmp(argv[1], "pm")) {
+    mode = FIT_PM_NOCLIP;
+  }
+  if (!strcasecmp(argv[1], "pm-irls")) {
+    mode = FIT_PM_IRLS;
+  }
+  if (!strcasecmp(argv[1], "plx")) {
+    mode = FIT_PLX_NOCLIP;
+  }
+  if (!strcasecmp(argv[1], "plx-irls")) {
+    mode = FIT_PLX_IRLS;
+  }
+
   // plan_tests (14);
 
@@ -34,5 +59,6 @@
     srand48(A);
   }
-  gauss_init (2048);
+
+  gaussdev_init ();
 
   FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, 0);
@@ -40,14 +66,5 @@
   int i;
 
-  // XXX try a single star to see if it looks good:
-  if (0) {
-    mkstar (fitStats, 3.0, 2.0, 0.2, 0.4, 0.1, N_POINTS, N_OUTLIERS);
-    for (i = 0; i < fitStats->Npoints; i++) {
-      fprintf (stdout, "%f %f\n", 3600*(fitStats->points[i].R - 3.0), 3600*(fitStats->points[i].D - 2.0));
-    }
-    exit (0);
-  }
-
-  int Nstars = 1000;
+  int Nstars = N_STARS;
   for (i = 0; i < Nstars; i++) {
 
@@ -58,5 +75,19 @@
     double uR  = 0.5*(drand48() - 0.5);
     double uD  = 0.5*(drand48() - 0.5);
-    double plx = 0.0*(drand48() + 0.0);
+
+    double plx;
+    switch (mode) {
+      case FIT_PM_IRLS:
+      case FIT_PM_NOCLIP:
+	plx = 0.0*(drand48() + 0.0);
+	break;
+      case FIT_PLX_IRLS:
+      case FIT_PLX_NOCLIP:
+	plx = 0.5*(drand48() + 0.0);
+	break;
+      default:
+	myAbort("oops");
+    }
+
 
     mkstar (fitStats, Ro, Do, uR, uD, plx, N_POINTS, N_OUTLIERS);
@@ -68,8 +99,19 @@
     FitAstromResultInit (&fitPM);
 
-    if (0) {
-      FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
-    } else {
-      FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, 0);
+    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;
+      default:
+	myAbort("oops");
     }
 
@@ -83,4 +125,6 @@
   }
   // return exit_status();
+  gaussdev_free();
+
   exit (0);
 }
@@ -103,6 +147,6 @@
     ParFactor (&pR, &pD, Ro, Do, points[i].T);
 
-    double dR = rnd_gauss(0.0, POS_ERROR);
-    double dD = rnd_gauss(0.0, POS_ERROR);
+    double dR = gaussdev_rnd(0.0, POS_ERROR);
+    double dD = gaussdev_rnd(0.0, POS_ERROR);
 
     points[i].R = Ro + (dR + plx*pR + uR*(points[i].T - MJD_REF_YRS))/3600/dcos(Do);
@@ -135,71 +179,2 @@
   return TRUE;
 }
-
-static int Ngaussint = 0;
-static double *gaussint;
-
-extern double drand48();
-
-double gaussian (double x, double mean, double sigma) {
-
-  double f;
-
-  f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
-
-  return (f);
-
-}
-
-/* integrate a gaussian from -5 sigma to +5 sigma */
-void gauss_init (int Nbin) {
- 
-  int i;
-  double val, x, dx, dx1, dx2, dx3, df;
-  double mean, sigma;
- 
-  /* no need to generate this if it already exists */
-  if (Ngaussint == Nbin) return;
-
-  // A = time(NULL);
-  // // XXX this is expensive if called a lot (1 sec min)
-  // // for (B = 0; A == time(NULL); B++);
-  // B = A + 10000;
-  // srand48(B);
- 
-  Ngaussint = Nbin;
-  ALLOCATE (gaussint, double, Ngaussint + 1);
-
-  val = 0;
-  dx = 1.0 / Ngaussint;
-  dx1 = dx / 3.0;
-  dx2 = 2.0*dx/3.0;
-  dx3 = dx;
-  mean = 0.0;
-  sigma = 1.0;
- 
-  for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx)  {
-    df = (3.0*gaussian(x    , mean, sigma) + 
-          9.0*gaussian(x+dx1, mean, sigma) +
-          9.0*gaussian(x+dx2, mean, sigma) + 
-          3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0);
-    val += df;
-    if (val > (i + 0.5) / (double) Ngaussint) {
-      gaussint[i] = x + dx / 2.0;
-      i++;
-    }
-  }
-}
-
-double rnd_gauss (double mean, double sigma) {
- 
-  int i;
-  double y;
- 
-  y = drand48();
-  i = Ngaussint*y;
-  y = gaussint[i]*sigma + mean;
- 
-  return (y);
- 
-}
-