Index: trunk/Ohana/src/opihi/cmd.data/fit1d_irls.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/fit1d_irls.c	(revision 39639)
+++ trunk/Ohana/src/opihi/cmd.data/fit1d_irls.c	(revision 39640)
@@ -2,9 +2,7 @@
 
 // These should probably be tunable:
-# define MAX_ITERATIONS 10
 # define FIT_TOLERANCE 1e-4
 # define FLT_TOLERANCE 1e-6
 # define WEIGHT_THRESHOLD 0.3
-# define NBOOTSTRAP 100
 
 int fit_least_squares (double *s, double **b, double **c, int nterm, int mterm, opihi_flt *x, opihi_flt *y, opihi_flt *dy, opihi_flt *wt, char *mask, int Npts);
@@ -32,7 +30,51 @@
   }
 
+  Vector *wtVec = NULL;
+  if ((N = get_argument (argc, argv, "-wt"))) {
+    remove_argument (N, &argc, argv);
+    if (!(wtVec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
+    remove_argument (N, &argc, argv);
+  }
+
+  Vector *yfitVec = NULL;
+  if ((N = get_argument (argc, argv, "-yfit"))) {
+    remove_argument (N, &argc, argv);
+    if (!(yfitVec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
+    remove_argument (N, &argc, argv);
+  }
+  Vector *mkVec = NULL;
+  if ((N = get_argument (argc, argv, "-mask"))) {
+    remove_argument (N, &argc, argv);
+    if (!(mkVec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
+    remove_argument (N, &argc, argv);
+  }
+
+  int NBOOTSTRAP = 100;
+  if ((N = get_argument (argc, argv, "-Nboot"))) {
+    remove_argument (N, &argc, argv);
+    NBOOTSTRAP = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  int MaxIterations = 10;
+  if ((N = get_argument (argc, argv, "-max-iterations"))) {
+    remove_argument (N, &argc, argv);
+    MaxIterations = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  float WeightThreshold = 0.3;
+  if ((N = get_argument (argc, argv, "-weight-threshold"))) {
+    remove_argument (N, &argc, argv);
+    WeightThreshold = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  int UseMean = TRUE;
+  if ((N = get_argument (argc, argv, "-use-median"))) {
+    remove_argument (N, &argc, argv);
+    UseMean = FALSE;
+  }
+
   if ((argc != 4) || !dyvec) {
     // require a y error value
-    gprint (GP_ERR, "USAGE: fit1d x y order -dy wt [-quiet/-q] [-clip Nsigma Niter]\n");
+    gprint (GP_ERR, "USAGE: fit1d x y order -dy wt [-quiet/-q] [-wt wtOut] [-yfit yfit] [-mask maskOut]\n");
     return (FALSE);
   }
@@ -53,4 +95,6 @@
   }
 
+  if (mkVec) ResetVector (mkVec, OPIHI_INT, xvec[0].Nelements);
+
   /* nterm is number of polynomial terms, starting at x^0 */
   order = atof (argv[3]);
@@ -87,5 +131,5 @@
 
   int converged = FALSE;
-  for (int iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations++) {
+  for (int iterations = 0; !converged && (iterations < MaxIterations); iterations++) {
 
     for (i = 0; i < xvec[0].Nelements; i++) {
@@ -116,11 +160,34 @@
   }
       
+  /* XXX TEST 
+  for (i = 0; i < nterm; i++) {
+    snprintf (name, 64, "C%d", i);
+    set_variable (name, b[i][0]);
+    if (!Quiet) gprint (GP_ERR, "%f x^%d ", b[i][0], i);
+    snprintf (name, 64, "dC%d", i);
+    set_variable (name, sqrt(c[i][i]));
+    if (!Quiet) gprint (GP_ERR, "%f x^%d ", sqrt(c[i][i]), i);
+  }
+  if (!Quiet) gprint (GP_ERR, "\n");
+  return TRUE; */
+
+  for (i = 0; i < nterm; i++) {
+    b_old[i][0] = b[i][0];
+  }
+
+
   // calculate the weight thresholds to mask the bad points:
   double Sum_W = 0.0;
+  ALLOCATE_PTR (wtlist, double, xvec[0].Nelements);
   for (i = 0; i < xvec[0].Nelements; i++) {
     wt[i] = weight_cauchy (yoff[i] / dy[i]);
+    wtlist[i] = wt[i];
     Sum_W += wt[i];
   }
-  double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * xvec[0].Nelements);
+  dsort (wtlist, xvec[0].Nelements);
+  int midpt = 0.5 * xvec[0].Nelements;
+  double WtMedian = (xvec[0].Nelements % 2) ? wtlist[midpt] : 0.5*(wtlist[midpt] + wtlist[midpt-1]);
+  double WtThreshold = UseMean ? WeightThreshold * Sum_W / (1.0 * xvec[0].Nelements) : WeightThreshold * WtMedian;
+  free (wtlist);
 
   // generate the unmasked subset
@@ -132,5 +199,9 @@
   int Npoints = 0;
   for (i = 0; i < xvec[0].Nelements; i++) {
-    if (wt[i] < WtThreshold) continue;
+    if (mkVec) mkVec->elements.Int[i] = 0;
+    if (wt[i] < WtThreshold) {
+      if (mkVec) mkVec->elements.Int[i] = 1;
+      continue;
+    }
      xunmask[Npoints] =  xvec[0].elements.Flt[i];
      yunmask[Npoints] =  yvec[0].elements.Flt[i];
@@ -184,6 +255,6 @@
   for (i = 0; i < nterm; i++) {
     snprintf (name, 64, "C%d", i);
-    set_variable (name, b[i][0]);
-    if (!Quiet) gprint (GP_ERR, "%f x^%d ", b[i][0], i);
+    set_variable (name, b_old[i][0]);
+    if (!Quiet) gprint (GP_ERR, "%f x^%d ", b_old[i][0], i);
   }
   if (!Quiet) gprint (GP_ERR, "\n");
@@ -216,6 +287,15 @@
   free (s);
   free (yoff);
-  free (yfit);
-  free (wt);
+
+  if (wtVec) {
+    SetVectorValues (wtVec, wt, OPIHI_FLT, xvec[0].Nelements);
+  } else {
+    free (wt);
+  }
+  if (yfitVec) {
+    SetVectorValues (yfitVec, yfit, OPIHI_FLT, xvec[0].Nelements);
+  } else {
+    free (yfit);
+  }
   
   free (xunmask);
Index: trunk/Ohana/src/opihi/include/dvomath.h
===================================================================
--- trunk/Ohana/src/opihi/include/dvomath.h	(revision 39639)
+++ trunk/Ohana/src/opihi/include/dvomath.h	(revision 39640)
@@ -150,4 +150,5 @@
 int           CopyVector            PROTO((Vector *out, Vector *in));
 int           ResetVector           PROTO((Vector *vec, char type, int Nelements));
+int           SetVectorValues       PROTO((Vector *vec, void *data, char type, int Nelements));
 int           SetVector             PROTO((Vector *vec, char type, int Nelements));
 int           CastVector            PROTO((Vector *vec, char type));
Index: trunk/Ohana/src/opihi/lib.shell/VectorOps.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/VectorOps.c	(revision 39639)
+++ trunk/Ohana/src/opihi/lib.shell/VectorOps.c	(revision 39640)
@@ -229,4 +229,20 @@
   } else {
     REALLOCATE (vec[0].elements.Int, opihi_int, MAX(1, Nelements));
+    vec[0].type = OPIHI_INT;
+  }
+  return TRUE;
+}
+
+int SetVectorValues (Vector *vec, void *data, char type, int Nelements) {
+
+  // a vector can only have >= 0 elements
+  vec[0].Nelements = MAX(Nelements,0);
+  if (type == OPIHI_FLT) {
+    free (vec[0].elements.Flt);
+    vec[0].elements.Flt = (opihi_flt *) data;
+    vec[0].type = OPIHI_FLT;
+  } else {
+    free (vec[0].elements.Int);
+    vec[0].elements.Int = (opihi_int *) data;
     vec[0].type = OPIHI_INT;
   }