Index: trunk/Ohana/src/opihi/cmd.data/Makefile
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 39640)
+++ trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 39926)
@@ -124,4 +124,5 @@
 $(SRC)/resize.$(ARCH).o		\
 $(SRC)/relocate.$(ARCH).o	\
+$(SRC)/rndseed.$(ARCH).o		\
 $(SRC)/roll.$(ARCH).o		\
 $(SRC)/rotate.$(ARCH).o	\
Index: trunk/Ohana/src/opihi/cmd.data/init.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/init.c	(revision 39640)
+++ trunk/Ohana/src/opihi/cmd.data/init.c	(revision 39926)
@@ -113,4 +113,5 @@
 int reindex          PROTO((int, char **));
 int relocate         PROTO((int, char **));
+int rndseed          PROTO((int, char **));
 int roll             PROTO((int, char **));
 int rotate           PROTO((int, char **));
@@ -293,4 +294,5 @@
   {1, "relocate",     relocate,         "set graphics/image window position"},
   {1, "roll",         roll,             "roll image to new start point"},
+  {1, "rndseed",      rndseed,          "set the pseudo-random seed"},
   {1, "rotate",       rotate,           "rotate image"},
   {1, "save",         save,             "save an SAOimage style image overlay"},
Index: trunk/Ohana/src/opihi/cmd.data/match2d.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/match2d.c	(revision 39640)
+++ trunk/Ohana/src/opihi/cmd.data/match2d.c	(revision 39926)
@@ -1,6 +1,9 @@
 # include "data.h"
 
-int find_matches2d (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2);
+int find_matches2d (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch);
 int find_matches2d_closest (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index);
+
+int find_matches2d_sphere (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch);
+int find_matches2d_sphere_closest (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index);
 
 // match2d (X1) (Y1) (X2) (Y2) (Radius) [-index1 (index1)] [-index2 (index2)] [-nomatch1 nomatch1] [-nomatch2 nomatch2]
@@ -23,4 +26,14 @@
   }
 
+  int SPHERE_DISTANCE = FALSE;
+  if ((N = get_argument (argc, argv, "-sphere"))) {
+    remove_argument (N, &argc, argv);
+    SPHERE_DISTANCE = TRUE;
+  }
+  if ((N = get_argument (argc, argv, "-sky"))) {
+    remove_argument (N, &argc, argv);
+    SPHERE_DISTANCE = TRUE;
+  }
+
   if ((N = get_argument (argc, argv, "-index1"))) {
     remove_argument (N, &argc, argv);
@@ -37,4 +50,16 @@
   } else {
     if ((index2 = SelectVector ("index2", ANYVECTOR, TRUE)) == NULL) return (FALSE);    
+  }
+
+  
+  Vector *radiusMatch = NULL;
+  if ((N = get_argument (argc, argv, "-radius"))) {
+    if (CLOSEST) {
+      gprint (GP_ERR, "error: -radius and -closest are currently incompatible\n");
+      return (FALSE);
+    }
+    remove_argument (N, &argc, argv);
+    if ((radiusMatch = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);    
+    remove_argument (N, &argc, argv);
   }
 
@@ -84,11 +109,19 @@
   }
 
-  if (CLOSEST) {
+  if (SPHERE_DISTANCE) {
+    if (CLOSEST) {
+      find_matches2d_sphere_closest (X1vec, Y1vec, X2vec, Y2vec, Radius, index1);
+      find_matches2d_sphere_closest (X2vec, Y2vec, X1vec, Y1vec, Radius, index2);
+    } else {
+      find_matches2d_sphere (X1vec, Y1vec, X2vec, Y2vec, Radius, index1, index2, radiusMatch);
+    }
+  } else {
+    if (CLOSEST) {
       find_matches2d_closest (X1vec, Y1vec, X2vec, Y2vec, Radius, index1);
       find_matches2d_closest (X2vec, Y2vec, X1vec, Y1vec, Radius, index2);
-  } else {
-      find_matches2d (X1vec, Y1vec, X2vec, Y2vec, Radius, index1, index2);
-  }
-
+    } else {
+      find_matches2d (X1vec, Y1vec, X2vec, Y2vec, Radius, index1, index2, radiusMatch);
+    }
+  }
   return (TRUE);
 
@@ -108,9 +141,13 @@
   gprint (GP_ERR, "use 'reindex' to generate new vectors based on these index vectors\n");
 
+  gprint (GP_ERR, "if -sphere or -sky is supplied, (x1,y1) and (x2,y2) are treaded as (ra,dec) or (long,lat) pairs in degrees\n");
+
+  gprint (GP_ERR, "if -radius (vector) is supplied, the vector will be filled with the distance between the matched pairs\n");
+  gprint (GP_ERR, "  not valid with -closest\n");
   return FALSE;
 }
 
 // we are not defining a relative offset DX,DY for now
-int find_matches2d (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2) {
+int find_matches2d (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch) {
   
   off_t i, j, first_j, I, J, *N1, *N2, Nmatch, NMATCH, DMATCH;
@@ -122,4 +159,5 @@
   ResetVector (index1, OPIHI_INT, NMATCH);
   ResetVector (index2, OPIHI_INT, NMATCH);
+  if (radiusMatch) ResetVector (radiusMatch, OPIHI_FLT, NMATCH);
 
   ALLOCATE (N1, off_t, X1->Nelements);
@@ -154,4 +192,5 @@
 	index1->elements.Int[Nmatch] = I;
 	index2->elements.Int[Nmatch] = J;
+	if (radiusMatch) radiusMatch->elements.Flt[Nmatch] = dR;
 
 	// XXX track matches 1 and 2 with internal vector, save new nomatch index vectors
@@ -163,4 +202,5 @@
 	  REALLOCATE (index1->elements.Int, opihi_int, NMATCH);
 	  REALLOCATE (index2->elements.Int, opihi_int, NMATCH);
+	  if (radiusMatch) { REALLOCATE (radiusMatch->elements.Flt, opihi_flt, NMATCH); }
 	}
       }
@@ -171,4 +211,5 @@
   index1->Nelements = Nmatch;
   index2->Nelements = Nmatch;
+  if (radiusMatch) radiusMatch->Nelements = Nmatch;
 
   free (N1);
@@ -246,2 +287,231 @@
   return (TRUE);
 }
+
+double gcdist (double r1, double d1, double r2, double d2) {
+  double num,den;
+  r1 *= RAD_DEG;
+  d1 *= RAD_DEG;
+  r2 *= RAD_DEG;
+  d2 *= RAD_DEG;
+
+  num = sqrt(pow((cos(d2) * sin(r2 - r1)),2) +
+	     pow((cos(d1) * sin(d2) -
+		  sin(d1) * cos(d2) * cos(r2 - r1)),2));
+  den = (sin(d1) * sin(d2) + cos(d1) * cos(d2) * cos(r2 - r1));
+  return(atan2(num,den) * (180 / M_PI));
+}
+
+typedef struct {
+  double sD;
+  double cD;
+  double sR;
+  double cR;
+} Match2D_PreCalc;
+
+double gcdist_PreCalc_v0 (Match2D_PreCalc *p1, Match2D_PreCalc *p2) {
+  double num,den;
+
+  num = sqrt(pow((p2->cD * (p2->sR*p1->cR - p1->sR*p2->cR)),2) +
+	     pow((p1->cD * p2->sD -
+		  p1->sD * p2->cD * (p2->cR*p1->cR + p2->sR*p1->sR)),2));
+  den = (p1->sD * p2->sD + p1->cD * p2->cD * (p2->cR*p1->cR + p2->sR*p1->sR));
+  return(atan2(num,den) * (180 / M_PI));
+}
+
+// we are not defining a relative offset DX,DY for now
+double gcdist_PreCalc (Match2D_PreCalc *p1, Match2D_PreCalc *p2) {
+  double num,den;
+
+  // double Qd = p1->sD * p1->cD * p2->sD * p2->cD;
+  // double Qr = p1->sR * p1->cR * p2->sR * p2->cR;
+
+  //  double Xa 
+  //    = SQ(p2->cD * p1->cR * p2->sR) 
+  //    + SQ(p2->cD * p1->sR * p2->cR) -
+  //    - 2 * SQ(p2->cD) * Qr;
+
+  double cdR = (p2->cR*p1->cR + p2->sR*p1->sR);
+
+  double Xa = SQ(p2->cD * p2->sR * p1->cR - p2->cD * p1->sR * p2->cR);
+  double Xb = SQ(p1->cD * p2->sD - p1->sD * p2->cD * cdR);
+
+  num = sqrt(Xa + Xb);
+  den = (p1->sD * p2->sD + p1->cD * p2->cD * cdR);
+  return(atan2(num,den) * (180 / M_PI));
+}
+
+// we are not defining a relative offset DX,DY for now
+int find_matches2d_sphere (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch) {
+  
+  off_t i, j, first_j, I, J, *N1, *N2, Nmatch, NMATCH, DMATCH;
+  double dY, dR;
+
+  NMATCH = MAX(MAX(0.05*X1->Nelements, 0.05*X2->Nelements), 1000);
+  DMATCH = NMATCH;
+
+  ResetVector (index1, OPIHI_INT, NMATCH);
+  ResetVector (index2, OPIHI_INT, NMATCH);
+  if (radiusMatch) ResetVector (radiusMatch, OPIHI_FLT, NMATCH);
+
+  ALLOCATE (N1, off_t, X1->Nelements);
+  ALLOCATE (N2, off_t, X2->Nelements);
+
+  ALLOCATE_PTR (A1, Match2D_PreCalc, X1->Nelements);
+  ALLOCATE_PTR (A2, Match2D_PreCalc, X2->Nelements);
+
+  for (i = 0; i < X1->Nelements; i++) { 
+    A1[i].sR = sin(RAD_DEG*X1->elements.Flt[i]);
+    A1[i].cR = cos(RAD_DEG*X1->elements.Flt[i]);
+    A1[i].sD = sin(RAD_DEG*Y1->elements.Flt[i]);
+    A1[i].cD = cos(RAD_DEG*Y1->elements.Flt[i]);
+    N1[i] = i; 
+  }
+  for (i = 0; i < X2->Nelements; i++) { 
+    A2[i].sR = sin(RAD_DEG*X2->elements.Flt[i]);
+    A2[i].cR = cos(RAD_DEG*X2->elements.Flt[i]);
+    A2[i].sD = sin(RAD_DEG*Y2->elements.Flt[i]);
+    A2[i].cD = cos(RAD_DEG*Y2->elements.Flt[i]);
+    N2[i] = i; 
+  }
+
+  // sort from one pole to the other
+  sort_coords_indexonly (Y1->elements.Flt, X1->elements.Flt, N1, X1->Nelements);
+  sort_coords_indexonly (Y2->elements.Flt, X2->elements.Flt, N2, X2->Nelements);
+
+  Nmatch = 0;
+  for (i = j = 0; (i < X1->Nelements) && (j < X2->Nelements);) {
+    I = N1[i];
+    J = N2[j];
+
+    // we can use dY as minimal requirement: if dY > Radius, we are too far apart
+    dY = Y1->elements.Flt[I] - Y2->elements.Flt[J];
+
+    if (dY <= -1.02*Radius) { i++; continue; }
+    if (dY >= +1.02*Radius) { j++; continue; }
+
+    // look for all matches of list2() to list1(i)
+    first_j = j;
+    for (j = first_j; (dY > -1.02*Radius) && (j < X2->Nelements); j++) {
+      J = N2[j];
+
+      dR = gcdist_PreCalc (&A1[I], &A2[J]);
+      // dR = gcdist (X1->elements.Flt[I], Y1->elements.Flt[I], X2->elements.Flt[J], Y2->elements.Flt[J]);
+
+      if (dR < Radius) {
+	index1->elements.Int[Nmatch] = I;
+	index2->elements.Int[Nmatch] = J;
+	if (radiusMatch) radiusMatch->elements.Flt[Nmatch] = dR;
+
+	// XXX track matches 1 and 2 with internal vector, save new nomatch index vectors
+	// after this loop
+
+	Nmatch ++;
+	if (Nmatch >= NMATCH) {
+	  NMATCH += DMATCH;
+	  REALLOCATE (index1->elements.Int, opihi_int, NMATCH);
+	  REALLOCATE (index2->elements.Int, opihi_int, NMATCH);
+	  if (radiusMatch) { REALLOCATE (radiusMatch->elements.Flt, opihi_flt, NMATCH); }
+	}
+      }
+    }
+    j = first_j;
+    i++;
+  }
+  index1->Nelements = Nmatch;
+  index2->Nelements = Nmatch;
+  if (radiusMatch) radiusMatch->Nelements = Nmatch;
+
+  free (A1);
+  free (A2);
+
+  free (N1);
+  free (N2);
+
+  return (TRUE);
+}
+
+// we are not defining a relative offset DX,DY for now
+int find_matches2d_sphere_closest (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index) {
+  
+  off_t i, j, Jmin, Ji, I, J, *N1, *N2;
+  double dY, dR, Rmin;
+
+  ResetVector (index, OPIHI_INT, X1->Nelements);
+  for (i = 0; i < index->Nelements; i++) { index->elements.Int[i] = -1; }
+
+  ALLOCATE (N1, off_t, X1->Nelements);
+  ALLOCATE (N2, off_t, X2->Nelements);
+
+  ALLOCATE_PTR (A1, Match2D_PreCalc, X1->Nelements);
+  ALLOCATE_PTR (A2, Match2D_PreCalc, X2->Nelements);
+
+  for (i = 0; i < X1->Nelements; i++) { 
+    A1[i].sR = sin(RAD_DEG*X1->elements.Flt[i]);
+    A1[i].cR = cos(RAD_DEG*X1->elements.Flt[i]);
+    A1[i].sD = sin(RAD_DEG*Y1->elements.Flt[i]);
+    A1[i].cD = cos(RAD_DEG*Y1->elements.Flt[i]);
+    N1[i] = i; 
+  }
+  for (i = 0; i < X2->Nelements; i++) { 
+    A2[i].sR = sin(RAD_DEG*X2->elements.Flt[i]);
+    A2[i].cR = cos(RAD_DEG*X2->elements.Flt[i]);
+    A2[i].sD = sin(RAD_DEG*Y2->elements.Flt[i]);
+    A2[i].cD = cos(RAD_DEG*Y2->elements.Flt[i]);
+    N2[i] = i; 
+  }
+
+  // sort from one pole to the other
+  sort_coords_indexonly (Y1->elements.Flt, X1->elements.Flt, N1, X1->Nelements);
+  sort_coords_indexonly (Y2->elements.Flt, X2->elements.Flt, N2, X2->Nelements);
+
+  for (i = j = 0; (i < X1->Nelements) && (j < X2->Nelements);) {
+    I = N1[i];
+    J = N2[j];
+
+    // we can use dY as minimal requirement: if dY > Radius, we are too far apart
+    dY = Y1->elements.Flt[I] - Y2->elements.Flt[J];
+
+    if (dY <= -1.02*Radius) { 
+      // no match in list 2 to this entry
+      index->elements.Int[I] = -1; // (probably not needed --- index is init'ed above)o
+      i++; 
+      continue; 
+    }
+    if (dY >= +1.02*Radius) { j++; continue; }
+
+    // look for all matches of list2() to list1(i)
+    Jmin = -1;
+    Rmin = Radius;
+    for (Ji = j; (dY > -1.02*Radius) && (Ji < X2->Nelements); Ji++) {
+      J = N2[Ji];
+
+      dR = gcdist_PreCalc (&A1[I], &A2[J]);
+      // dR = gcdist (X1->elements.Flt[I], Y1->elements.Flt[I], X2->elements.Flt[J], Y2->elements.Flt[J]);
+
+      if (dR < Rmin) {
+	Rmin = dR;
+	Jmin = J;
+      }
+    }
+    
+
+    // no match in list 2 to this entry
+    if (Jmin == -1) {
+      index->elements.Int[I] = -1;
+      i++;
+      continue;
+    }
+    index->elements.Int[I] = Jmin;
+    i++;
+  }
+
+  free (A1);
+  free (A2);
+
+  free (N1);
+  free (N2);
+
+  return (TRUE);
+}
+
+
Index: trunk/Ohana/src/opihi/cmd.data/rndseed.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/rndseed.c	(revision 39926)
+++ trunk/Ohana/src/opihi/cmd.data/rndseed.c	(revision 39926)
@@ -0,0 +1,18 @@
+# include "data.h"
+
+int rndseed (int argc, char **argv) {
+  
+  if (argc < 2) {
+    gprint (GP_ERR, "USAGE: rndseed (value)\n");
+    return (FALSE);
+  }
+
+  /* init srand for rnd numbers elsewhere */
+  long A = atol (argv[1]);
+  srand48(A);
+
+  return (TRUE);
+}
+
+ 
+  
Index: trunk/Ohana/src/opihi/cmd.data/stats-new.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/stats-new.c	(revision 39640)
+++ trunk/Ohana/src/opihi/cmd.data/stats-new.c	(revision 39926)
@@ -129,7 +129,5 @@
   ALLOCATE (values, float, Nsample);
     
-  A = time(NULL);
-  for (B = 0; A == time(NULL); B++);
-  srand48(B);
+  // srand48() is called by startup.c
  
   *buffer = (float *) matrix[0].buffer;
Index: trunk/Ohana/src/opihi/cmd.data/test/fit1d_irls.sh
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/test/fit1d_irls.sh	(revision 39926)
+++ trunk/Ohana/src/opihi/cmd.data/test/fit1d_irls.sh	(revision 39926)
@@ -0,0 +1,347 @@
+
+list tests
+ test1
+ test2
+ test3
+ test4
+ test5
+ test6
+ test7
+ test8
+end
+
+## NOTE: this is not a TAP test
+macro mean.test
+  if ($0 != 2)
+    echo "USAGE: mean.test (errorbar)"
+    break
+  end
+
+  create x 0 100
+
+  delete C0fit dC0fit Cnfit Mfit dMfit
+  for i 0 500
+    gaussdev y x[] 0 1.0   ; # mean of 0, sigma of 1
+    
+    set dy = 1.0 + zero(x) ; # set error to be correct
+    fit1d_irls -q x y 0 -dy dy -wt wt -mask mask -use-median -Nboot 1000
+
+    concat $C0    C0fit
+    concat $dC0  dC0fit
+
+    vstat -q y
+    concat $MEAN Mfit
+    concat {$SIGMA / sqrt($NPTS)} dMfit
+  end
+
+  vstat -q C0fit;  set MeanC  = $MEAN; set StdevC = $SIGMA
+  vstat -q dC0fit; set MeanDC = $MEAN
+  vstat -q Mfit;   set MeanM  = $MEAN; set StdevM = $SIGMA
+  vstat -q dMfit;  set MeanDM = $MEAN
+
+  fprintf "%7.3f %7.3f" $MeanC $MeanM
+  fprintf "%7.3f %7.3f" $StdevC $StdevM
+end
+
+macro cliptest.plot
+  if ($0 != 3)
+    echo "USAGE: cliptest.plot (errorbar) (output)"
+    break
+  end
+
+  set dy = $1 + zero(x) ; # set error to be correct
+  fit1d_irls x y 0 -dy dy -wt wt -mask mask -use-median
+
+  subset ygood = y where not(mask)
+
+  vstat y
+  vstat ygood
+
+  histogram wt Nwt 0 1 0.005 -range dw
+  vstat wt
+
+  subset wts = wt where (wt < 0.3*$MEDIAN)
+  echo "0.3 of MEDIAN : wts[]"
+
+  subset wts = wt where (wt < 0.3*$MEAN)
+  echo "0.3 of MEAN : wts[]"
+
+  dev -n 0; resize 1200 600
+  label -fn courier 24; 
+  section a 0.4 0.0 0.6 1.0; lim x y; clear; box -ypad 0.1 -xpad 3.2 -labelpadx 2.5 -labels 1001; label -x sequence +y value +x "red: clipped, blue: unclipped"; 
+  plot x y; plot -c blue -pt 7 x y where not(mask); plot -c red -pt 7 -lw 2 x y where mask
+
+  section b 0.0 0.0 0.4 1.0; lim dw Nwt; box +ypad 0.1 -xpad 3.2 -labelpadx 2.5; label -x "weight" -y "Npoints"; plot -x 1 dw Nwt
+  png -name $2
+end
+
+## NOTE: this is not a TAP test
+macro cliptest
+
+  create x 0 100
+  gaussdev y x[] 0 1.0   ; # mean of 0, sigma of 1
+
+  if (1)
+    y[10] = 5.0
+    y[20] = 7.0
+    y[30] = 4.0
+  end
+
+  cliptest.plot 1.0 irls.sample.v0.png
+
+  cliptest.plot 0.1 irls.sample.v1.png
+
+  cliptest.plot 0.01 irls.sample.v2.png
+end
+
+## NOTE: this is not a TAP test
+macro outlier.test
+  if ($0 != 2)
+    echo "USAGE: outlier.test (errorbar)"
+    break
+  end
+
+  create x 0 100
+
+  delete C0fit dC0fit Cnfit Mfit dMfit
+  for i 0 300
+    gaussdev y x[] 0 1.0   ; # mean of 0, sigma of 1
+    
+    if (0)
+      y[10] = 5.0
+      y[20] = 7.0
+      y[30] = 4.0
+    end
+    
+    set dy = $1 + zero(x) ; # set error to be correct
+    fit1d_irls -q x y 0 -dy dy -wt wt -mask mask -use-median -Nboot 1000
+    concat $C0 C0fit
+    concat $dC0 dC0fit
+    concat $Cnfit Cnfit
+
+    vstat -q y
+    concat $MEAN Mfit
+    concat {$SIGMA / sqrt($NPTS)} dMfit
+  end
+
+  vstat C0fit
+  vstat dC0fit
+  vstat Cnfit
+  vstat Mfit
+  vstat dMfit
+end
+
+# fit a line without errors
+macro test1
+ $PASS = 1
+ break -auto off
+
+ delete -q x y
+
+ create x 0 100
+ set y = 3 + 5*x
+ fit -q x y 1
+
+ if ($Cn != 1)
+   $PASS = 0
+ end
+ if (abs($C0 - 3) > 1e-5)
+   $PASS = 0
+ end
+ if (abs($C1 - 5) > 1e-5)
+   $PASS = 0
+ end
+end
+
+# fit a line with errors
+macro test2
+ $PASS = 1
+ break -auto off
+
+ delete -q x y dy
+
+ create x 0 100
+ set dy = 0.1*rnd(x) - 0.05
+ set y = 3 + 5*x + dy
+ fit -q x y 1
+
+ if ($Cn != 1)
+   $PASS = 0
+   echo "wrong number of elements : $Cn"
+ end
+ if (abs($C0 - 3) > 0.06)
+   $PASS = 0
+   echo "wrong value for C0 : $C0"
+ end
+ if (abs($C1 - 5) > 0.06)
+   $PASS = 0
+   echo "wrong value for C1 : $C1"
+ end
+end
+
+# fit a line with errors and weights
+macro test3
+ $PASS = 1
+ break -auto off
+
+ delete -q x y dy
+
+ create x 0 100
+ set dy = 0.1*rnd(x) - 0.05
+ set y = 3 + 5*x + dy
+ set dy = 0.1 + zero(x)
+ fit -q x y 1 -dy dy
+
+ if ($Cn != 1)
+   $PASS = 0
+ end
+ if (abs($C0 - 3) > 0.06)
+   $PASS = 0
+ end
+ if (abs($C1 - 5) > 0.06)
+   $PASS = 0
+ end
+end
+
+# fit a line with errors, weights, and outliers 
+macro test4
+ $PASS = 1
+ break -auto off
+
+ delete -q x y dy
+
+ create x 0 100
+ set dy = 0.1*rnd(x) - 0.05
+ set y = 3 + 5*x + dy
+ set dy = 0.1 + zero(x)
+ y[5] = 23
+ y[20] = -10
+ y[50] = 0.0
+ fit -q x y 1 -dy dy -clip 3 3
+
+ if ($Cn != 1)
+   $PASS = 0
+ end
+ if ($Cnv != 97)
+   $PASS = 0
+ end
+ if (abs($C0 - 3) > 0.06)
+   $PASS = 0
+ end
+ if (abs($C1 - 5) > 0.06)
+   $PASS = 0
+ end
+end
+
+# fit a quadratic without errors
+macro test5
+ $PASS = 1
+ break -auto off
+
+ delete -q x y
+
+ create x 0 100
+ set y = 3 + 5*x - 4*x^2
+ fit -q x y 2
+
+ if ($Cn != 2)
+   $PASS = 0
+ end
+ if (abs($C0 - 3) > 1e-5)
+   $PASS = 0
+ end
+ if (abs($C1 - 5) > 1e-5)
+   $PASS = 0
+ end
+ if (abs($C2 + 4) > 1e-5)
+   $PASS = 0
+ end
+end
+
+# fit a quadratic with errors
+macro test6
+ $PASS = 1
+ break -auto off
+
+ delete -q x y dy
+
+ create x 0 100
+ set dy = 0.1*rnd(x) - 0.05
+ set y = 3 + 5*x - 4*x^2 + dy
+ fit -q x y 2
+
+ if ($Cn != 2)
+   $PASS = 0
+ end
+ if (abs($C0 - 3) > 0.06)
+   $PASS = 0
+ end
+ if (abs($C1 - 5) > 0.06)
+   $PASS = 0
+ end
+ if (abs($C2 + 4) > 0.06)
+   $PASS = 0
+ end
+end
+
+# fit a quadratic with errors and weights
+macro test7
+ $PASS = 1
+ break -auto off
+
+ delete -q x y dy
+
+ create x 0 100
+ set dy = 0.1*rnd(x) - 0.05
+ set y = 3 + 5*x - 4*x^2 + dy
+ set dy = 0.1 + zero(x)
+ fit -q x y 2 -dy dy
+
+ if ($Cn != 2)
+   $PASS = 0
+ end
+ if (abs($C0 - 3) > 0.06)
+   $PASS = 0
+ end
+ if (abs($C1 - 5) > 0.06)
+   $PASS = 0
+ end
+ if (abs($C2 + 4) > 0.06)
+   $PASS = 0
+ end
+end
+
+# fit a quadratic with errors, weights, and outliers 
+macro test8
+ $PASS = 1
+ break -auto off
+
+ delete -q x y dy
+
+ create x 0 100
+ set dy = 0.1*rnd(x) - 0.05
+ set y = 3 + 5*x - 4*x^2 + dy
+ set dy = 0.1 + zero(x)
+ y[5] = 23
+ y[20] = -10
+ y[50] = 0.0
+
+ # it takes 4 iterations to successfully reject the outliers above...
+ fit -q x y 2 -dy dy -clip 3 4
+
+ if ($Cn != 2)
+   $PASS = 0
+ end
+ if ($Cnv != 97)
+   $PASS = 0
+ end
+ if (abs($C0 - 3) > 0.06)
+   $PASS = 0
+ end
+ if (abs($C1 - 5) > 0.06)
+   $PASS = 0
+ end
+ if (abs($C2 + 4) > 0.06)
+   $PASS = 0
+ end
+end