Index: trunk/Ohana/src/opihi/cmd.astro/astrom_ops.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/astrom_ops.c	(revision 39596)
+++ 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 39596)
+++ 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 39596)
+++ 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);