Index: trunk/Ohana/src/relastro/Makefile
===================================================================
--- trunk/Ohana/src/relastro/Makefile	(revision 39237)
+++ trunk/Ohana/src/relastro/Makefile	(revision 39238)
@@ -29,6 +29,7 @@
 $(SRC)/FitMosaic.$(ARCH).o           \
 $(SRC)/FitSimple.$(ARCH).o           \
-$(SRC)/FitAstromOps.$(ARCH).o           \
+$(SRC)/FitAstromOps.$(ARCH).o        \
 $(SRC)/FitPM.$(ARCH).o               \
+$(SRC)/FitPM_IRLS.$(ARCH).o          \
 $(SRC)/FitPMandPar.$(ARCH).o         \
 $(SRC)/FitPosPMfixed.$(ARCH).o       \
@@ -103,6 +104,7 @@
 $(SRC)/ConfigInit.$(ARCH).o	     \
 $(SRC)/FitSimple.$(ARCH).o           \
-$(SRC)/FitAstromOps.$(ARCH).o           \
+$(SRC)/FitAstromOps.$(ARCH).o        \
 $(SRC)/FitPM.$(ARCH).o               \
+$(SRC)/FitPM_IRLS.$(ARCH).o          \
 $(SRC)/FitPMandPar.$(ARCH).o         \
 $(SRC)/FitPosPMfixed.$(ARCH).o       \
Index: trunk/Ohana/src/relastro/include/relastro.h
===================================================================
--- trunk/Ohana/src/relastro/include/relastro.h	(revision 39237)
+++ trunk/Ohana/src/relastro/include/relastro.h	(revision 39238)
@@ -134,7 +134,12 @@
   double **A;
   double **B;
+  double **Cov;
+  double *Beta;
+  double *Beta_prev;
   int Nterms;
 } FitAstromData;
 
+// XXX do we need doubles for all of these?  I actually only have of order 100 of these
+// allocated at a time, so size is not an issue.
 typedef struct {
   double X, dX;
@@ -143,4 +148,7 @@
   double D, dD;
   double T, dT;
+  double Wx, Wy;
+  double rx, ry;
+  double u;
   double pR;
   double pD;
@@ -148,4 +156,5 @@
   double C_red;
   int measure;
+  int mask;
 } FitAstromPoint;
 
@@ -739,2 +748,11 @@
 int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average);
 void AstromErrorSetLoop (int Nloop, int isImageMode);
+
+int FitPM_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);
+
+double weight_cauchy (double x);
+double dpsi_cauchy (double x);
Index: trunk/Ohana/src/relastro/src/FitAstromOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitAstromOps.c	(revision 39237)
+++ trunk/Ohana/src/relastro/src/FitAstromOps.c	(revision 39238)
@@ -93,6 +93,10 @@
 
   /* do I need to do this as 2 2x2 matrix equations? */
-  fit->A = array_init (Nterms, Nterms);
-  fit->B = array_init (Nterms, 1);
+  fit->B   = array_init (Nterms, 1);
+  fit->A   = array_init (Nterms, Nterms);
+  fit->Cov = array_init (Nterms, Nterms);
+
+  ALLOCATE (fit->Beta, double, Nterms);
+  ALLOCATE (fit->Beta_prev, double, Nterms);
   fit->Nterms = Nterms;
 
@@ -106,4 +110,9 @@
   array_free (fit->A, fit->Nterms);
   array_free (fit->B, fit->Nterms);
+  array_free (fit->Cov, fit->Nterms);
+
+  free (fit->Beta);
+  free (fit->Beta_prev);
+
   free (fit);
   return;
@@ -126,4 +135,13 @@
   object->C_red  = 0.0;
   object->measure= -1;
+
+  object->Wx     = 1.0;
+  object->Wy     = 1.0;
+
+  object->rx     = 0.0;
+  object->ry     = 0.0;
+  object->u      = 0.0;
+
+  object->mask   = 0;
   return;
 }
Index: trunk/Ohana/src/relastro/src/FitPM_IRLS.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPM_IRLS.c	(revision 39237)
+++ trunk/Ohana/src/relastro/src/FitPM_IRLS.c	(revision 39238)
@@ -1,285 +1,238 @@
 # include "relastro.h"
 
+// These should probably be tunable:
+# define MAX_ITERATIONS 10
+# define FIT_TOLERANCE 1e-4
+# define WEIGHT_THRESHOLD 0.3
+
 /* do we want an init function which does the alloc and a clear function to free? */
-int FitPMonly_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
+int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
 
   int i,j;
 
-  static double **A, **B;
-
-  double chisq, Xf, Yf;
-
-  double **Cov;
-  double *Beta, *Beta_prev;
-  
-  double sigma_ols, sigma_hat;
-  double *Wx, *Wy;
-  double *rx, *ry;
-  //  double *ux, *uy;
-  double *u;
   int dof = 2 * Npoints - 4;
   int p   = 4;
   int n   = 2 * Npoints;
-  double tolerance;
-  int converged;
-  int iterations;
-  
-  /* do I need to do this as 2 2x2 matrix equations? */
-  if (A == NULL) {
-    ALLOCATE (A, double *, 4);
-    ALLOCATE (B, double *, 4);
-    for (i = 0; i < 4; i++) {
-      ALLOCATE (A[i], double, 4);
-      ALLOCATE (B[i], double, 1);
-      memset (A[i], 0, 4*sizeof(double));
-      memset (B[i], 0, 1*sizeof(double));
-    }
-  }
-
-  // things we need
-  ALLOCATE (Cov, double *, 4);
-  for (i = 0; i < 4; i++) {
-    ALLOCATE ( Cov[i], double, 4);
-  }
-
-  ALLOCATE(Beta, double, 4);
-  ALLOCATE(Beta_prev, double, 4);
-  ALLOCATE(Wx, double, Npoints);
-  ALLOCATE(Wy, double, Npoints);
-  ALLOCATE(rx,  double, Npoints);
-  ALLOCATE(ry,  double, Npoints);
-  ALLOCATE(u,  double, Npoints);
+  
+  // 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++) {
-    Wx[i] = 1 / dX[i];
-    Wy[i] = 1 / dY[i];
+    points[i].Wx = 1 / points[i].dX;
+    points[i].Wy = 1 / points[i].dY;
   }
   
   // Solve OLS equation  
-  if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npoints,
-		   A,B,VERBOSE)) {
-    // Handle fail case
+  if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
+    myAbort ("handle failures, please!");
     return(FALSE);
   }
 
   // Calculate r vector of residuals and least squares sigma
-  sigma_ols = 0.0;
-  for (i = 0; i < Npoints; i++) {
-    rx[i] = X[i] - (T[i] * B[1][0] + B[0][0]);
-    ry[i] = Y[i] - (T[i] * B[3][0] + B[2][0]);
-    //    u[i] = r[i] /
-    sigma_ols += SQ(rx[i]) + SQ(ry[i]);
-
+  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 / dof);
 
-  // Save OLS covariance;
+  // 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++) {
-      Cov[i][j] = A[i][j];
-    }
-  }
-
-  // Save Beta
-  for (i = 0; i < 4; i++) {
-    Beta[i] = B[i][0];
+      data->Cov[i][j] = data->A[i][j];
+    }
+    data->Beta[i] = data->B[i][0];
   }
 
   // Iterately reweight and solve
-  converged = FALSE;
-  iterations = 0;
-  do {
+  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 < 4; i ++) {
-      Beta_prev[i] = Beta[i];
-    }
-
-    // Assign W
+      data->Beta_prev[i] = data->Beta[i];
+    }
+
+    // Assign weights based on the deviation
     for (i = 0; i < Npoints; i++) {
-      Wx[i] = weight_cauchy(rx[i] / dX[i]);
-      Wy[i] = weight_cauchy(ry[i] / dY[i]);
+      points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
+      points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     }    
 
-    // Solve
-    if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npoints,
-		     A,B,VERBOSE)) {
-      // Handle fail case
+    // Solve with the new weights
+    if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
+      myAbort ("handle failures, please!");
       return(FALSE);
     }
 
+    // store the new Beta.
     for (i = 0; i < 4; i++) {
-      Beta[i] = B[i][0];
-    }
-
-    // r
-    sigma_hat = 0.0;
+      data->Beta[i] = data->B[i][0];
+    }
+
+    // calculate the residuals:
     for (i = 0; i < Npoints; i++) {
-      rx[i] = X[i] - (T[i] * B[1][0] + B[0][0]);
-      ry[i] = Y[i] - (T[i] * B[3][0] + B[2][0]);
-      u[i] = sqrt(SQ(rx[i] / dX[i]) + SQ(ry[i] / dY[i]));
-    }
-    sigma_hat = MedianAbsDeviation(u,Npoints) / 0.6745;
+      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;
-    tolerance = 1e-4;  // This should probably be tunable.
     for (i = 0; i < 4; i++) {
-      if (fabs(Beta[i] - Beta_prev[i]) > tolerance * abs(Beta[i])) {
+      if (fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) {
 	converged = FALSE;
       }
     }
-
-    iterations++;
-    if (iterations >= 10) {
-      converged = TRUE;
-      // Throw a warning or something here.
-    }
-    
-  } while (!converged);
-
-  double ax, ay;
-  double bx, by;
-  double lambda;
-  double sigma_robust_x, sigma_robust_y;
-  double sigma_final_x,  sigma_final_y;
-  double Sum_Wx, Sum_Wy;
-  
-  Sum_Wx = 0.0;
-  Sum_Wy = 0.0;
-  ax = 0.0; ay = 0.0;
-  bx = 0.0; by = 0.0;
-  lambda = 0.0;
-  for (i = 0; i < Npoints; i++) {
-    Wx[i] = weight_cauchy(rx[i] / dX[i]);
-    Wy[i] = weight_cauchy(ry[i] / dY[i]);
-    
-    ax += dpsi_cauchy(rx[i] / dX[i]);
-    ay += dpsi_cauchy(ry[i] / dY[i]);
-
-    bx += SQ(Wx[i]);
-    by += SQ(Wy[i]);
-
-    Sum_Wx += Wx[i];
-    Sum_Wy += Wy[i];
-  }
-  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);
-  
-  sigma_robust_x = lambda * sqrt(bx) * sigma_hat * 2.385 / ax;
-  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)
-  sigma_final_x  = MAX(SQ(sigma_robust_x), (n * SQ(sigma_robust_x) + SQ(p * sigma_ols)) / (n + SQ(p)));
-  sigma_final_y  = MAX(SQ(sigma_robust_y), (n * SQ(sigma_robust_y) + SQ(p * sigma_ols)) / (n + SQ(p)));
-
-  for (i = 0; i < 4; i++) {
-    for (j = 0; j < 4; j++) {
-      // This uses the original OLS covariance.
-      if ((i < 2)&&(j < 2)) { // Upper portion
-	Cov[i][j] *= sigma_final_x;
-      }
-      else if ((i > 1)&&(j > 1)) { // Lower portion
-	Cov[i][j] *= sigma_final_y;
-      }
-      else { // Cross term
-	Cov[i][j] *= sqrt(sigma_final_x * sigma_final_y);
-      }
-    }
-  }
-
-  // Finish.
-  fit[0].Ro = Beta[0];
-  fit[0].uR = Beta[1];
-  fit[0].Do = Beta[2];
-  fit[0].uD = Beta[3];
-  
-  fit[0].dRo = sqrt(Cov[0][0]);
-  fit[0].duR = sqrt(Cov[1][1]);
-  fit[0].dDo = sqrt(Cov[2][2]);
-  fit[0].duD = sqrt(Cov[3][3]);
-
-  // Sort out the final weight threshold.
-
-  // add up the chi square for the fit
-  chisq = 0.0;
+  }
+  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[9].dRo *= sigma_final_x;
+    fit[9].duR *= sigma_final_x;
+
+    fit[9].dDo *= sigma_final_y;
+    fit[9].duD *= 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;
+    }
+  }
+
+  // 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 ((Wx[i] > 0.1 * Sum_Wx / (1.0 * Npoints))||
-	(Wy[i] > 0.1 * Sum_Wy / (1.0 * Npoints))) {
-      Xf = fit[0].Ro + fit[0].uR*T[i];
-      Yf = fit[0].Do + fit[0].uD*T[i];
-      chisq += SQ(X[i] - Xf) / SQ(dX[i]);
-      chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
-      fit[0].Nfit += 1;
-    }
-    // if (VERBOSE) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
-  }
-  //  fit[0].Nfit = Npoints;
-
-  // the reduced chisq is divided by (Ndof = 2*Npoints - 4)
-  fit[0].chisq = chisq / (2.0*Npoints - 4.0);
+    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);
+    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);
 }
 
-int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, double **A, double **B, int VERBOSE) {
-
-  int i,j;
+int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
+
+  int i;
   double 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 += WX[i];
-    Wy += WY[i];
-
-    Tx += T[i]*WX[i];
-    Ty += T[i]*WY[i];
-    
-    Tx2 += SQ(T[i])*WX[i];
-    Ty2 += SQ(T[i])*WY[i];
-    
-    Xs += X[i]*WX[i];
-    Ys += Y[i]*WY[i];
-
-    XT += X[i]*T[i]*WX[i];
-    YT += Y[i]*T[i]*WY[i];
+
+  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;
+
+    Tx += TWx;
+    Ty += TWy;
+    
+    Tx2 += points[i].T*TWx;
+    Ty2 += points[i].T*TWy;
+    
+    Xs += points[i].X*points[i].Wx;
+    Ys += points[i].Y*points[i].Wy;
+
+    XT += points[i].X*TWx;
+    YT += points[i].Y*TWy;
   }
 
   // X^T W X
-  A[0][0] = Wx;
-  A[0][1] = Tx;
-
-  A[1][0] = Tx;
-  A[1][1] = Tx2;
-
-  A[2][2] = Wy;
-  A[2][3] = Ty;
-
-  A[3][2] = Ty;
-  A[3][3] = Ty2;
+  data->A[0][0] = Wx;
+  data->A[0][1] = Tx;
+
+  data->A[1][0] = Tx;
+  data->A[1][1] = Tx2;
+  data->A[2][2] = Wy;
+  data->A[2][3] = Ty;
+  data->A[3][2] = Ty;
+  data->A[3][3] = Ty2;
 
   // X^T W Y
-  B[0][0] = Xs;
-  B[1][0] = XT;
-  B[2][0] = Ys;
-  B[3][0] = YT;
-
-  if (!dgaussjordan ((double **)A, (double **)B, 4, 1)) {
+  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");
-    if (VERBOSE == 2) {
-      for (i = 0; i < 4; i++) {
-	for (j = 0; j < 4; j++) {
-	  fprintf (stderr, "%e ", A[i][j]);
-	}
-	fprintf (stderr, " : %e\n", A[i][0]);
+    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;
   }
 
-  // A => (X^T W X)^{-1}
-  // B => beta
-  
+  fit->Ro = data->B[0][0];
+  fit->uR = data->B[1][0];
+  fit->Do = data->B[2][0];
+  fit->uD = data->B[3][0];
+  fit->p  = 0.0;
+
   return TRUE;
 }
@@ -299,34 +252,33 @@
 // median absolute deviation
 // MAD = median(abs(x - median(x)))
-double MedianAbsDeviation(double *in, int N) {
+double MedianAbsDeviation(FitAstromPoint *points, int Npoints) {
+
   double *x;
   double median = 0.0;
   int i;
   
-  ALLOCATE(x,double,N);
-  for (i = 0; i < N; i++) {
-    x[i] = in[i];
-  }
-
-  dsort(x,N);
-
-  if (N % 2) {
-    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
+  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*N)];
-  }
-
-  for (i = 0; i < N; i++ ) {
+    median = x[(int)(0.5*Npoints)];
+  }
+
+  for (i = 0; i < Npoints; i++ ) {
     x[i] = fabs(x[i] - median);
   }
-
-  dsort(x,N);
-
-  if (N % 2) {
-    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
+  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*N)];
-  }
-
-  return(median);
-}
+    median = x[(int)(0.5*Npoints)];
+  }
+
+  return median;
+}
Index: trunk/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 39237)
+++ trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 39238)
@@ -153,11 +153,14 @@
     if (fitStats->NfitAlloc == 1) {
       // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling:
-      FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
-      // FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
+      if (1) {
+	FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
+      } else {
+	FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, VERBOSE);
+      }
     } else {
       fitStats->Nfit = 0;
       for (k = 0; k < fitStats->NfitAlloc; k++) {
 	BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
-	// if (!FitPM_IRLS (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
+	// if (!FitPM_IRLS (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints, VERBOSE)) continue;
 	if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
 	fitStats->Nfit ++;
@@ -207,5 +210,5 @@
 	BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
 	// FitPMandPar_IRLS (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints);
-	FitPMandPar_IRLS (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints);
+	FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints);
 	fitStats->Nfit ++;
       }