Index: trunk/Ohana/src/opihi/cmd.astro/csystem.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/csystem.c	(revision 39640)
+++ trunk/Ohana/src/opihi/cmd.astro/csystem.c	(revision 39926)
@@ -125,6 +125,6 @@
   opihi_flt *yptr = yvec[0].elements.Flt;
 
-  opihi_flt *uxptr = uxvec ? xvec[0].elements.Flt : NULL;
-  opihi_flt *uyptr = uyvec ? yvec[0].elements.Flt : NULL;
+  opihi_flt *uxptr = uxvec ? uxvec[0].elements.Flt : NULL;
+  opihi_flt *uyptr = uyvec ? uyvec[0].elements.Flt : NULL;
 
   for (i = 0; i < xvec[0].Nelements; i++, xptr++, yptr++) {
Index: trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c	(revision 39640)
+++ trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c	(revision 39926)
@@ -158,6 +158,6 @@
     // now that the mask has been updated, we need to recalculate mean epoch and positions
     // XXX make this conditional on actually masking unmasked points above
-    PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
-    PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
+    // PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
+    // PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
 
     PlxFitData sample;
@@ -218,13 +218,37 @@
   // fprintf (stderr, "%f +/- %f | %f %f\n", fit.p, fit.dp, fit.uR, fit.uD);
 
-/*
-  FILE *f = fopen ("test.pf.dat", "w");
+  Vector *dRresPMP, *dDresPMP, *dRresPLX, *dDresPLX;
+
+  // save fit residuals (with only pm removed, and pm and plx removed)
+  if ((dRresPMP = SelectVector ("dRresPMP", ANYVECTOR, TRUE)) == NULL) ESCAPE ("cannot generate vector %s\n", "dRresPMP");
+  if ((dDresPMP = SelectVector ("dDresPMP", ANYVECTOR, TRUE)) == NULL) ESCAPE ("cannot generate vector %s\n", "dDresPMP");
+  if ((dRresPLX = SelectVector ("dRresPLX", ANYVECTOR, TRUE)) == NULL) ESCAPE ("cannot generate vector %s\n", "dRresPLX");
+  if ((dDresPLX = SelectVector ("dDresPLX", ANYVECTOR, TRUE)) == NULL) ESCAPE ("cannot generate vector %s\n", "dDresPLX");
+    
+  ResetVector (dRresPMP, OPIHI_FLT, Ntotal);
+  ResetVector (dDresPMP, OPIHI_FLT, Ntotal);
+  ResetVector (dRresPLX, OPIHI_FLT, Ntotal);
+  ResetVector (dDresPLX, OPIHI_FLT, Ntotal);
+  
   for (i = 0; i < Ntotal; i++) {
-    double Xf = fit.Ro + fit.uR*fitdata.t[i] + fit.p*fitdata.pX[i];
-    double Yf = fit.Do + fit.uD*fitdata.t[i] + fit.p*fitdata.pY[i];
-    fprintf (f, "%f : %f %f : %f %f : %f : %f %f : %f %f\n", T[i], R[i], D[i], Xf, Yf, fitdata.t[i], fitdata.X[i], fitdata.Y[i], fitdata.pX[i], fitdata.pY[i]);
-  }
-  fclose (f);
-*/
+    
+    double x0, y0;
+    RD_to_XY (&x0, &y0, R[i], D[i], &coords);
+
+    double pX0, pY0;
+    ParFactor (&pX0, &pY0, R[i], D[i], T[i]);
+
+    double t0 = (T[i] - Tmean)/365.25;
+
+    double Xpmp = fit.Ro + fit.uR*t0 + fit.p*pX0;
+    double Ypmp = fit.Do + fit.uD*t0 + fit.p*pY0;
+    double Xplx = fit.Ro + fit.uR*t0;
+    double Yplx = fit.Do + fit.uD*t0;
+
+    dRresPMP->elements.Flt[i] = x0 - Xpmp;
+    dDresPMP->elements.Flt[i] = y0 - Ypmp;
+    dRresPLX->elements.Flt[i] = x0 - Xplx;
+    dDresPLX->elements.Flt[i] = y0 - Yplx;
+  }
 
   // fprintf (stderr, "Roff, Doff: %f, %f; dRo, dDo: %f, %f\n", fit.Ro, fit.Do, fit.dRo, fit.dDo);
Index: trunk/Ohana/src/opihi/cmd.astro/spex1dgas.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/spex1dgas.c	(revision 39640)
+++ trunk/Ohana/src/opihi/cmd.astro/spex1dgas.c	(revision 39926)
@@ -108,12 +108,7 @@
   IDList = NULL;
   XoList = NULL;
-  { 
-    // init random numbers
-    long A, B;
-    A = time(NULL);
-    for (B = 0; A == time(NULL); B++);
-    srand48(B);
-  }
  
+  // srand48() is called by startup.c
+
   if (argc != 11) goto usage;
 
Index: trunk/Ohana/src/opihi/cmd.astro/spex2dgas.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/spex2dgas.c	(revision 39640)
+++ trunk/Ohana/src/opihi/cmd.astro/spex2dgas.c	(revision 39926)
@@ -120,12 +120,6 @@
   float XoMax, YoMax;
 
-  { 
-    // init random numbers
-    long A, B;
-    A = time(NULL);
-    for (B = 0; A == time(NULL); B++);
-    srand48(B);
-  }
- 
+  // srand48() is called by startup.c
+
   if (argc != 9) goto usage;
 
Index: trunk/Ohana/src/opihi/cmd.astro/spexseq.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/spexseq.c	(revision 39640)
+++ trunk/Ohana/src/opihi/cmd.astro/spexseq.c	(revision 39926)
@@ -90,13 +90,7 @@
   int Nobject = 0;
 
-  { 
-    // init random numbers
-    long A, B;
-    A = time(NULL);
-    for (B = 0; A == time(NULL); B++);
-    srand48(B);
-  }
- 
   if (argc != 6) goto usage;
+
+  // srand48() is called by startup.c
 
   // XXX enforce matching lengths on the three vectors
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
Index: trunk/Ohana/src/opihi/dvo/imbox.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/imbox.c	(revision 39640)
+++ trunk/Ohana/src/opihi/dvo/imbox.c	(revision 39926)
@@ -15,4 +15,17 @@
 
   if (!style_args (&graphmode, &argc, argv, &kapa)) return FALSE;
+
+  char *xaxis = NULL;
+  if ((N = get_argument (argc, argv, "-xaxis"))) {
+    remove_argument (N, &argc, argv);
+    xaxis = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  char *yaxis = NULL;
+  if ((N = get_argument (argc, argv, "-yaxis"))) {
+    remove_argument (N, &argc, argv);
+    yaxis = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
 
   SOLO_PHU = FALSE;
@@ -61,14 +74,29 @@
     // XXX currently, image uses an unsigned short for NX,XY. this is rather restrictive
     // and needs to be at least checked.
-    haveNx = gfits_scan (&header, "IMNAXIS1",   "%d", 1, &Nx);
-    haveNy = gfits_scan (&header, "IMNAXIS2",   "%d", 1, &Ny);
+    haveNx = FALSE;
+    if (xaxis) {
+      haveNx = gfits_scan (&header, xaxis, "%d", 1, &Nx);
+    }
+    if (!haveNx) {
+      haveNx = gfits_scan (&header, "IMNAXIS1",   "%d", 1, &Nx);
+    } 
+    if (!haveNx) {
+	haveNx = gfits_scan (&header, "ZNAXIS1",   "%d", 1, &Nx);
+    }
+    if (!haveNx) {
+	haveNx = gfits_scan (&header, "NAXIS1",   "%d", 1, &Nx);
+    }
 
-    if (!haveNx && !haveNy) {
-	haveNx = gfits_scan (&header, "ZNAXIS1",   "%d", 1, &Nx);
+    haveNy = FALSE;
+    if (yaxis) {
+      haveNy = gfits_scan (&header, yaxis, "%d", 1, &Ny);
+    }
+    if (!haveNy) {
+      haveNy = gfits_scan (&header, "IMNAXIS2",   "%d", 1, &Ny);
+    }
+    if (!haveNy) {
 	haveNy = gfits_scan (&header, "ZNAXIS2",   "%d", 1, &Ny);
     }
-
-    if (!haveNx && !haveNy) {
-	haveNx = gfits_scan (&header, "NAXIS1",   "%d", 1, &Nx);
+    if (!haveNy) {
 	haveNy = gfits_scan (&header, "NAXIS2",   "%d", 1, &Ny);
     }
@@ -137,4 +165,6 @@
   free (Xvec.elements.Flt);
   free (Yvec.elements.Flt);
+  if (xaxis) free (xaxis);
+  if (yaxis) free (yaxis);
   return (TRUE);
 
Index: trunk/Ohana/src/opihi/dvo/imdense.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/imdense.c	(revision 39640)
+++ trunk/Ohana/src/opihi/dvo/imdense.c	(revision 39926)
@@ -4,5 +4,4 @@
 int imdense (int argc, char **argv) {
   
-  long A, B;
   off_t i, Nimage;
   int kapa, N, status, NPTS;
@@ -27,7 +26,5 @@
   Rmax = graphmode.coords.crval1 + 182.0;
   
-  A = time(NULL);
-  for (B = 0; A == time(NULL); B++);
-  srand48(B);
+  // srand48() is called by startup.c
 
   N = 0;
Index: trunk/Ohana/src/opihi/dvo/skycat.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/skycat.c	(revision 39640)
+++ trunk/Ohana/src/opihi/dvo/skycat.c	(revision 39926)
@@ -78,6 +78,10 @@
       if (table) {
 	int hostID = (regions[i][0].hostFlags & DATA_USE_BCK) ? regions[i][0].backupID : regions[i][0].hostID;
-	int index = table->index[hostID];
-	snprintf (hostfile, 1024, "%s/%s.cpt", table->hosts[index].pathname, regions[i][0].name);
+	if (hostID) {
+	  int index = table->index[hostID];
+	  snprintf (hostfile, 1024, "%s/%s.cpt", table->hosts[index].pathname, regions[i][0].name);
+	} else {
+	  strcpy (hostfile, skylist[0].filename[i]);
+	}
       } else {
 	strcpy (hostfile, skylist[0].filename[i]);
Index: trunk/Ohana/src/opihi/lib.data/style_args.c
===================================================================
--- trunk/Ohana/src/opihi/lib.data/style_args.c	(revision 39640)
+++ trunk/Ohana/src/opihi/lib.data/style_args.c	(revision 39926)
@@ -29,5 +29,5 @@
   if ((N = get_argument (*argc, argv, "-lt"))) {
     remove_argument (N, argc, argv);
-    graphmode[0].ltype = atof(argv[N]);
+    graphmode[0].ltype = KapaLineTypeFromString(argv[N]);
     remove_argument (N, argc, argv);
   }
@@ -39,5 +39,5 @@
   if ((N = get_argument (*argc, argv, "-pt"))) {
     remove_argument (N, argc, argv);
-    graphmode[0].ptype = atof(argv[N]);
+    graphmode[0].ptype = KapaPointStyleFromString(argv[N]);
     remove_argument (N, argc, argv);
   }
@@ -63,5 +63,5 @@
   if ((N = get_argument (*argc, argv, "-x"))) {
     remove_argument (N, argc, argv);
-    graphmode[0].style = atof(argv[N]);
+    graphmode[0].style = KapaPlotStyleFromString(argv[N]);
     remove_argument (N, argc, argv);
   }
Index: trunk/Ohana/src/opihi/lib.shell/startup.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/startup.c	(revision 39640)
+++ trunk/Ohana/src/opihi/lib.shell/startup.c	(revision 39926)
@@ -51,4 +51,15 @@
     gfits_set_unsign_mode (FALSE);
     
+    set_variable ("M_PI", M_PI);
+    set_variable ("M_E",  M_E);
+    set_variable ("M_c",  299792459.0); // meter / second
+    set_variable ("M_c_cgs", 29979245900.0); // cm / second
+
+    set_variable ("M_h",  6.62607004e-34); // meter^2 kg / second (J s)
+    set_variable ("M_h_cgs",  6.62607004e-27); // erg s
+
+    set_variable ("M_kB",  1.38064853e-23); // J / K
+    set_variable ("M_kB_cgs",  1.38064853e-16); // erg / K
+
   /* check history file permission */
   {