Index: trunk/Ohana/src/opihi/cmd.astro/astrom_ops.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/astrom_ops.c	(revision 39609)
+++ trunk/Ohana/src/opihi/cmd.astro/astrom_ops.c	(revision 39610)
@@ -54,4 +54,8 @@
   ALLOCATE (data->Wx, double, N);
   ALLOCATE (data->Wy, double, N);
+  ALLOCATE (data->Qx, double, N);
+  ALLOCATE (data->Qy, double, N);
+  ALLOCATE (data->qx, double, N);
+  ALLOCATE (data->qy, double, N);
   ALLOCATE (data->index, int, N);
   return TRUE;
@@ -68,4 +72,8 @@
   FREE (data->Wx);
   FREE (data->Wy);
+  FREE (data->Qx);
+  FREE (data->Qy);
+  FREE (data->qx);
+  FREE (data->qy);
   FREE (data->index);
 }
Index: trunk/Ohana/src/opihi/cmd.astro/fitplx.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/fitplx.c	(revision 39609)
+++ trunk/Ohana/src/opihi/cmd.astro/fitplx.c	(revision 39610)
@@ -597,17 +597,20 @@
   dsort (values, Npts);
 
-  if (Npts % 2) {
-    int Ncenter = Npts / 2;
-    *median = values[Ncenter];
-  } else {
-    int Ncenter = Npts / 2 - 1;
-    *median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
-  }
-
-  double Slo = VectorFractionInterpolate (values, 0.158655, Npts);
-  double Shi = VectorFractionInterpolate (values, 0.841345, Npts);
-
-  *sigma = (Shi - Slo) / 2.0;
-
+  if (median) {
+    if (Npts % 2) {
+      int Ncenter = Npts / 2;
+      *median = values[Ncenter];
+    } else {
+      int Ncenter = Npts / 2 - 1;
+      *median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
+    }
+  }
+    
+  if (sigma) {
+    double Slo = VectorFractionInterpolate (values, 0.158655, Npts);
+    double Shi = VectorFractionInterpolate (values, 0.841345, Npts);
+    *sigma = (Shi - Slo) / 2.0;
+  }
+  
   return TRUE;
 }
Index: trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c	(revision 39609)
+++ trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c	(revision 39610)
@@ -119,4 +119,6 @@
 			 fitdata.t, fitdata.pX, fitdata.pY,
 			 fitdata.Wx, fitdata.Wy,
+			 fitdata.Qx, fitdata.Qy,
+			 fitdata.qx, fitdata.qy,
 			 fitdata.Npts, max_iterations, outlier_limit, binning_step, VERBOSE)) {
     return FALSE;
@@ -204,10 +206,12 @@
     }
 
-    // now calculate median and sigma for each vector
-    VectorRobustStats (pvec,  &fit.p,  &fit.dp);
-    VectorRobustStats (uRvec, &fit.uR, &fit.duR);
-    VectorRobustStats (uDvec, &fit.uD, &fit.duD);
-    VectorRobustStats (Rvec,  &fit.Ro, &fit.dRo);
-    VectorRobustStats (Dvec,  &fit.Do, &fit.dDo);
+    // now calculate robust sigma for each vector
+    VectorRobustStats (pvec,  NULL, &fit.dp);
+    VectorRobustStats (uRvec, NULL, &fit.duR);
+    VectorRobustStats (uDvec, NULL, &fit.duD);
+    VectorRobustStats (Rvec,  NULL, &fit.dRo);
+    VectorRobustStats (Dvec,  NULL, &fit.dDo);
+
+    PlxFitDataFree (&sample);    
   }
 
@@ -250,5 +254,5 @@
   set_variable ("Nfit",  fit.Nfit);
 
-
+  PlxFitDataFree (&fitdata);
   
   return (TRUE);
@@ -257,5 +261,5 @@
 /* do we want an init function which does the alloc and a clear function to free? */
 int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD,
-		      double *Wx, double *Wy,
+		      double *Wx, double *Wy, double *Qx, double *Qy, double *qx, double *qy,
 		      int Npts, int max_iterations, double outlier_limit, double binning_step, int VERBOSE) {
 
@@ -270,5 +274,4 @@
 
   double sigma_ols, sigma_hat; // Sigma estimates.
-  //  double *Wx, *Wy;             // Weight vectors.  Not errors.
   double *rx, *ry;             // Deviation from model
   double *u;                   // Deviation magnitude
@@ -305,6 +308,6 @@
   // Convert the measurement errors into the initial weights.
   for (i = 0; i < Npts; i++) {
-    Wx[i] = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
-    Wy[i] = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
+    Qx[i] = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
+    Qy[i] = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
   }
   
@@ -313,6 +316,5 @@
   // Solve OLS equation
   if (binning_step == 0.0) {
-    if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,
-			 A,B,VERBOSE)) {
+    if (!weighted_LS_PLX(T,pR,pD,X,Qx,Y,Qy,Npts,A,B,VERBOSE)) {
       // Handle fail case
       return(FALSE);
@@ -329,5 +331,5 @@
     double *WYbin = NULL;
 
-    if (!bin_points_PLX(T, pR, pD, X, Wx, Y, Wy, Npts,
+    if (!bin_points_PLX(T, pR, pD, X, Qx, Y, Qy, Npts,
 			&Tbin, &pRbin, &pDbin, &Xbin, &WXbin, &Ybin, &WYbin, &Nbins,  binning_step) ) {
       return(FALSE);
@@ -396,8 +398,10 @@
       Wx[i] = weight_cauchy(rx[i] / dX[i]);
       Wy[i] = weight_cauchy(ry[i] / dY[i]);
+      qx[i] = Qx[i] * Wx[i];
+      qy[i] = Qy[i] * Wy[i];
     }
 
     // Solve
-    if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,A,B,VERBOSE)) {
+    if (!weighted_LS_PLX(T,pR,pD,X,qx,Y,qy,Npts,A,B,VERBOSE)) {
       // Handle fail case
       return(FALSE);
Index: trunk/Ohana/src/opihi/include/astro.h
===================================================================
--- trunk/Ohana/src/opihi/include/astro.h	(revision 39609)
+++ trunk/Ohana/src/opihi/include/astro.h	(revision 39610)
@@ -20,6 +20,10 @@
   double *dX;
   double *dY;
-  double *Wx;
-  double *Wy;
+  double *Wx; // modified weight based on distance from fit
+  double *Wy; // modified weight based on distance from fit
+  double *Qx; // raw error-based weight (1/dX^2)
+  double *Qy;
+  double *qx; // modified error-based weight (1/dX^2) * Wx
+  double *qy;
   int *index;
   int Npts;
@@ -59,5 +63,5 @@
 
 int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, double binning_step, int VERBOSE);
-int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *Wx, double *Wy, int Npts, int max_iterations, double outlier_limit, double binning_step, int VERBOSE);
+int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *Wx, double *Wy, double *Qx, double *Qy, double *qx, double *qy, int Npts, int max_iterations, double outlier_limit, double binning_step, int VERBOSE);
 
 int weighted_LS_PLX (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE);
Index: trunk/Ohana/src/relastro/include/relastro.h
===================================================================
--- trunk/Ohana/src/relastro/include/relastro.h	(revision 39609)
+++ trunk/Ohana/src/relastro/include/relastro.h	(revision 39610)
@@ -163,5 +163,7 @@
   double D, dD;
   double T, dT;
-  double Wx, Wy;
+  double Qx, Qy; // unmodified error-based weight (1/dX^2)
+  double qx, qy; // IRLS-modified error-based weight (1/dX^2)
+  double Wx, Wy; // IRLS weight factor
   double rx, ry;
   double u;
Index: trunk/Ohana/src/relastro/src/FitAstromOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitAstromOps.c	(revision 39609)
+++ trunk/Ohana/src/relastro/src/FitAstromOps.c	(revision 39610)
@@ -227,4 +227,9 @@
   object->Wy     = 1.0;
 
+  object->Qx     = 0.0;
+  object->Qy     = 0.0;
+  object->qx     = 0.0;
+  object->qy     = 0.0;
+
   object->rx     = 0.0;
   object->ry     = 0.0;
Index: trunk/Ohana/src/relastro/src/FitPMandPar.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPMandPar.c	(revision 39609)
+++ trunk/Ohana/src/relastro/src/FitPMandPar.c	(revision 39610)
@@ -1,3 +1,5 @@
 # include "relastro.h"
+# define OLD_METHOD 0
+
 int FitPMandPar_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
 int FitPMandPar_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
@@ -16,6 +18,15 @@
   // Convert the measurement errors into initial weights.
   for (i = 0; i < Npoints; i++) {
+# if (OLD_METHOD)
     points[i].Wx = 1.0 / SQ(points[i].dX);
     points[i].Wy = 1.0 / SQ(points[i].dY);
+# else
+    points[i].Wx = 1.0;
+    points[i].Wy = 1.0;
+    points[i].Qx = 1.0 / SQ(points[i].dX);
+    points[i].Qy = 1.0 / SQ(points[i].dY);
+    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
+    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
+# endif
   }
 
@@ -34,6 +45,15 @@
   // Convert the measurement errors into initial weights.
   for (i = 0; i < Npoints; i++) {
+# if (OLD_METHOD)
     points[i].Wx = 1.0 / SQ(points[i].dX);
     points[i].Wy = 1.0 / SQ(points[i].dY);
+    points[i].qx = points[i].Wx; // Wx, Wy start out at 1.0
+    points[i].qy = points[i].Wy; // Wx, Wy start out at 1.0
+# else
+    points[i].Qx = 1.0 / SQ(points[i].dX);
+    points[i].Qy = 1.0 / SQ(points[i].dY);
+    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
+    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
+# endif
   }
   
@@ -77,4 +97,6 @@
       points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
       points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
+      points[i].qx = points[i].Wx * points[i].Qx;
+      points[i].qy = points[i].Wy * points[i].Qy;
     }    
 
@@ -209,6 +231,11 @@
     Nfit ++;
 
+# if (OLD_METHOD)    
     wx = points[i].Wx;
     wy = points[i].Wy;
+# else
+    wx = points[i].qx;
+    wy = points[i].qy;
+# endif
 
     Wx += wx;
Index: trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c	(revision 39609)
+++ trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c	(revision 39610)
@@ -1,3 +1,4 @@
 # include "relastro.h"
+# define OLD_METHOD 1
 
 // These should probably be tunable:
@@ -15,6 +16,13 @@
   // Convert the measurement errors into initial weights.
   for (i = 0; i < Npoints; i++) {
+# if (OLD_METHOD)
     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);
+# else
+    points[i].Qx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX);
+    points[i].Qy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY);
+# endif
+    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
+    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
   }
   
@@ -58,4 +66,6 @@
       points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
       points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
+      points[i].qx = points[i].Wx * points[i].Qx;
+      points[i].qy = points[i].Wy * points[i].Qy;
     }    
 
@@ -138,6 +148,6 @@
       ay += dpsi_cauchy(points[i].ry / points[i].dY);
 
-      bx += SQ(points[i].Wx);
-      by += SQ(points[i].Wy);
+      bx += SQ(points[i].qx);
+      by += SQ(points[i].qy);
     }
     ax /= 1.0 * Npoints;  // mean(psi_dot(r))
@@ -207,6 +217,11 @@
     // 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]);
 
+# if (OLD_METHOD)
     wx = points[i].Wx;
     wy = points[i].Wy;
+# else 
+    wx = points[i].Wx;
+    wy = points[i].Wy;
+# endif
 
     Wx += wx;
Index: trunk/Ohana/src/relastro/src/fitobj.c
===================================================================
--- trunk/Ohana/src/relastro/src/fitobj.c	(revision 39609)
+++ trunk/Ohana/src/relastro/src/fitobj.c	(revision 39610)
@@ -121,5 +121,5 @@
   for (i = 0; i < Nstars; i++) {
 
-    double Ro  = 360.0*(drand48() - 0.5);
+    double Ro  = 360.0*drand48();
     double Phi = 2.0*(drand48() - 0.5);
     double Do  = DEG_RAD * asin(Phi);
Index: trunk/Ohana/src/relastro/src/fitpm.c
===================================================================
--- trunk/Ohana/src/relastro/src/fitpm.c	(revision 39609)
+++ trunk/Ohana/src/relastro/src/fitpm.c	(revision 39610)
@@ -155,5 +155,5 @@
   for (i = 0; i < Nstars; i++) {
 
-    double Ro  = 360.0*(drand48_cnt() - 0.5);
+    double Ro  = 360.0*drand48_cnt();
     double Phi = 2.0*(drand48_cnt() - 0.5);
     double Do  = DEG_RAD * asin(Phi);
@@ -181,7 +181,10 @@
     }
 
-    int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 points
-    int Noutliers_max = Npoints * F_OUTLIERS;
-    int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points
+    // int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 points
+    // int Noutliers_max = Npoints * F_OUTLIERS;
+    // int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points
+
+    int Npoints = N_POINTS; // minimum of 0 points
+    int Noutliers = Npoints * F_OUTLIERS;
 
     // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points
@@ -287,5 +290,6 @@
   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);
+	     Ro, Do, uR, uD, plx, 
+	     NAN, NAN, NAN, NAN, NAN, 0.0, NAN, NAN, NAN, NAN, NAN, 0, 0.0);
   }
   // return exit_status();
Index: trunk/Ohana/src/relastro/test/mana.sh
===================================================================
--- trunk/Ohana/src/relastro/test/mana.sh	(revision 39609)
+++ trunk/Ohana/src/relastro/test/mana.sh	(revision 39610)
@@ -91,6 +91,7 @@
   # foreach outliers 0 3 10 30 100 
 
-  foreach outliers 10
-    exec ../bin/fitpm.lin64 -Nstars 3000 -Noutliers $outliers -Npoints 100 -Nbootstrap 100 $1 > test.$1.dat
+  foreach outliers 0
+    echo ../bin/fitpm.lin64 -Nstars 3000 -Npoints 100 -Nbootstrap 100 $1 
+    exec ../bin/fitpm.lin64 -Nstars 3000 -Npoints 100 -Nbootstrap 100 $1 > test.$1.dat
     reload test.$1.dat 0
 
@@ -107,6 +108,6 @@
   data $1
   if ($2)
-    read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26
-    # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24
+    # read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26
+    read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24
   else
     read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22