Index: trunk/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 32632)
+++ trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 33662)
@@ -32,7 +32,8 @@
 $(SRC)/fixwrap.$(ARCH).o	   \
 $(SRC)/fixcols.$(ARCH).o	   \
+$(SRC)/fiximage.$(ARCH).o	   \
 $(SRC)/gauss.$(ARCH).o		   \
-$(SRC)/getvel.$(ARCH).o	   \
-$(SRC)/getlst.$(ARCH).o	   \
+$(SRC)/getvel.$(ARCH).o	           \
+$(SRC)/getlst.$(ARCH).o	           \
 $(SRC)/medianmap.$(ARCH).o	   \
 $(SRC)/mkgauss.$(ARCH).o	   \
@@ -51,12 +52,17 @@
 $(SRC)/precess.$(ARCH).o	   \
 $(SRC)/profile.$(ARCH).o	   \
-$(SRC)/radec.$(ARCH).o	   \
-$(SRC)/region.$(ARCH).o	   \
+$(SRC)/radec.$(ARCH).o	           \
+$(SRC)/region.$(ARCH).o	           \
 $(SRC)/rotcurve.$(ARCH).o	   \
 $(SRC)/scale.$(ARCH).o		   \
 $(SRC)/sexigesimal.$(ARCH).o	   \
 $(SRC)/spec.$(ARCH).o		   \
+$(SRC)/specpairfit.$(ARCH).o	   \
+$(SRC)/spexseq.$(ARCH).o	   \
+$(SRC)/spex1dgas.$(ARCH).o	   \
+$(SRC)/spex2dgas.$(ARCH).o	   \
+$(SRC)/mkclusters.$(ARCH).o	   \
 $(SRC)/star.$(ARCH).o		   \
-$(SRC)/transform.$(ARCH).o        \
+$(SRC)/transform.$(ARCH).o         \
 $(SRC)/imsub.$(ARCH).o		   \
 $(SRC)/imfit.$(ARCH).o		   \
Index: trunk/Ohana/src/opihi/cmd.astro/coords.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/coords.c	(revision 32632)
+++ trunk/Ohana/src/opihi/cmd.astro/coords.c	(revision 33662)
@@ -83,4 +83,5 @@
     }
     if (mode == PIXEL) {
+      Xin = ohana_normalize_angle_to_midpoint (Xin, coords.crval1);
       RD_to_XY (&Xout, &Yout, Xin, Yin, &coords);
       if (!Quiet) gprint (GP_LOG, "%7.2f %7.2f\n", Xout, Yout);
@@ -92,5 +93,9 @@
   if (mode == SKY) {
     for (i = 0; i < xvec[0].Nelements; i++) {
-      XY_to_RD (&xvec[0].elements.Flt[i], &yvec[0].elements.Flt[i], xvec[0].elements.Flt[i], yvec[0].elements.Flt[i], &coords);
+      double Xin = xvec[0].elements.Flt[i];
+      double Yin = yvec[0].elements.Flt[i];
+      xvec[0].elements.Flt[i] = NAN;
+      yvec[0].elements.Flt[i] = NAN;
+      XY_to_RD (&xvec[0].elements.Flt[i], &yvec[0].elements.Flt[i], Xin, Yin, &coords);
     }
     return (TRUE);
@@ -98,5 +103,10 @@
   if (mode == PIXEL) {
     for (i = 0; i < xvec[0].Nelements; i++) {
-      RD_to_XY (&xvec[0].elements.Flt[i], &yvec[0].elements.Flt[i], xvec[0].elements.Flt[i], yvec[0].elements.Flt[i], &coords);
+      double Xin = xvec[0].elements.Flt[i];
+      double Yin = yvec[0].elements.Flt[i];
+      xvec[0].elements.Flt[i] = NAN;
+      yvec[0].elements.Flt[i] = NAN;
+      Xin = ohana_normalize_angle_to_midpoint (Xin, coords.crval1);
+      RD_to_XY (&xvec[0].elements.Flt[i], &yvec[0].elements.Flt[i], Xin, Yin, &coords);
     }
     return (TRUE);
Index: trunk/Ohana/src/opihi/cmd.astro/fiximage.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/fiximage.c	(revision 33662)
+++ trunk/Ohana/src/opihi/cmd.astro/fiximage.c	(revision 33662)
@@ -0,0 +1,91 @@
+# include "astro.h"
+
+int fiximage (int argc, char **argv) {
+
+  int ix, iy, N;
+  Buffer *in, *ct, *mask;
+
+  int VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    VERBOSE = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: fiximage <data> <count> (mask)\n");
+    return (FALSE);
+  }
+
+  if ((in   = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((ct   = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((mask = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+
+  int Nx = in[0].matrix.Naxis[0];
+  int Ny = in[0].matrix.Naxis[1];
+
+  gfits_free_matrix (&mask[0].matrix);
+  gfits_free_header (&mask[0].header);
+  CreateBuffer (mask, Nx, Ny, -32, 0.0, 1.0);
+  strcpy (mask[0].file, "(empty)");
+  memset (mask[0].matrix.buffer, 0, Nx*Ny*sizeof(float));
+  // 0 -- init, cannot be repaired
+  // 1 -- does not need to be repaired
+  // 2 -- repair x only
+  // 4 -- repair y only
+  // 6 -- repair x & y only
+  
+  float *inB = (float *) in[0].matrix.buffer;
+  float *ctB = (float *) ct[0].matrix.buffer;
+  float *mkB = (float *) mask[0].matrix.buffer;
+
+  // loop over input ct and find pixels that can be repaired. 
+  for (iy = 1; iy < Ny - 1; iy++) {
+    for (ix = 1; ix < Nx - 1; ix++) {
+      
+      int Npix = ix + Nx*iy;
+
+      // does not need to be repaired
+      if (ctB[Npix]) { mkB[Npix] = 1.0; continue; }
+      
+      // if we have all 4 valid neighbors, we can fix
+      if (ctB[Npix - 1 ] && ctB[Npix + 1 ]) mkB[Npix] += 2.0;
+      if (ctB[Npix - Nx] && ctB[Npix + Nx]) mkB[Npix] += 4.0;
+    }
+  }
+	
+  // loop over input ct and repair the reparable pixels
+  for (iy = 1; iy < Ny - 1; iy++) {
+    for (ix = 1; ix < Nx - 1; ix++) {
+      
+      int Npix = ix + Nx*iy;
+
+      if (mkB[Npix] == 0.0) continue; // cannot be repaired
+      if (mkB[Npix] == 1.0) continue; // does not need to be repaired
+
+      float Vxm, Vxp, Vym, Vyp, Value;
+      if ((mkB[Npix] == 2) || (mkB[Npix] == 6)) {
+	Vxm = inB[Npix - 1 ] / ctB[Npix - 1 ];
+	Vxp = inB[Npix + 1 ] / ctB[Npix + 1 ];
+      }
+      if ((mkB[Npix] == 4) || (mkB[Npix] == 6)) {
+	Vym = inB[Npix - Nx] / ctB[Npix - Nx];
+	Vyp = inB[Npix + Nx] / ctB[Npix + Nx];
+      }
+
+      if (mkB[Npix] == 2) {
+	Value = 0.5*(Vxm + Vxp);
+      }
+      if (mkB[Npix] == 4) {
+	Value = 0.5*(Vym + Vyp);
+      }
+      if (mkB[Npix] == 6) {
+	Value = 0.25*(Vxm + Vxp + Vym + Vyp);
+      }
+
+      ctB[Npix] = 1;
+      inB[Npix] = Value;
+    }
+  }
+
+  return TRUE;
+}
Index: trunk/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 32632)
+++ trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 33662)
@@ -15,4 +15,5 @@
 int flux                    PROTO((int, char **));
 int fixwrap                 PROTO((int, char **));
+int fiximage                PROTO((int, char **));
 int fixcols                 PROTO((int, char **));
 int fixrows                 PROTO((int, char **));
@@ -45,4 +46,9 @@
 int sexigesimal             PROTO((int, char **));
 int spec                    PROTO((int, char **));
+int specpairfit             PROTO((int, char **));
+int spexseq                 PROTO((int, char **));
+int spex1dgas               PROTO((int, char **));
+int spex2dgas               PROTO((int, char **));
+int mkclusters              PROTO((int, char **));
 int star                    PROTO((int, char **));
 int times                   PROTO((int, char **));
@@ -64,4 +70,5 @@
   {1, "flux",        flux,         "flux in a convex contour"},
   {1, "fixwrap",     fixwrap,      "fix megacam over-wrapped pixels"},
+  {1, "fiximage",    fiximage,     "fix pixels in an image by interpolation"},
   {1, "fixcols",     fixcols,      "fix bad columns by comparing with others"},
   {1, "fixrows",     fixrows,      "fix bad rows by comparing with others"},
@@ -93,4 +100,9 @@
   {1, "sexigesimal", sexigesimal,  "convert to/from sexigesimal/decimal"},
   {1, "spec",        spec,         "extract a spectrum"},
+  {1, "specpairfit", specpairfit,  "fit spectrum to another spectrum"},
+  {1, "spexseq",     spexseq,      "generate the spectral sequence"},
+  {1, "spex1dgas",   spex1dgas,    "minimize distances in 1D"},
+  {1, "spex2dgas",   spex2dgas,    "minimize distances in 2D"},
+  {1, "mkclusters",  mkclusters,   "group spectra by distance"},
   {1, "star",        star,         "star stats at rough coords"},
   {1, "transform",   transform,    "geometric transformation of image"},
Index: trunk/Ohana/src/opihi/cmd.astro/mkclusters.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/mkclusters.c	(revision 33662)
+++ trunk/Ohana/src/opihi/cmd.astro/mkclusters.c	(revision 33662)
@@ -0,0 +1,233 @@
+# include "astro.h"
+
+typedef struct {
+  int *friend; // indices for all tested relationships
+  float *dist; // distance for this relationship
+  int Ndist;   // number of tested relationships
+  int NDIST;   // number of allocated relationships
+  int group;   // assigned group
+  int Nfriend; // number of friends brought into group
+} Object;
+
+typedef struct {
+  int *entry;
+  int Nentry;
+  int NENTRY;
+} Group;
+
+static Group  *group = NULL;
+static int Ngroup = 0;
+static int NGROUP = 0;
+
+static Object *object = NULL;
+static int Nobject = 0;
+// static int NOBJECT = 0;
+
+void add_to_group (int Nobj, float scale, int Nfriends) {
+
+  int j;
+
+  // fprintf (stderr, "add object %d to group %d\n", Nobj, Ngroup);
+
+  object[Nobj].group = Ngroup;
+
+  // add to this group
+  int N = group[Ngroup].Nentry;
+  group[Ngroup].entry[N] = Nobj;
+  group[Ngroup].Nentry ++;
+  CHECK_REALLOCATE (group[Ngroup].entry, int, group[Ngroup].NENTRY, group[Ngroup].Nentry, 100);
+
+  // add all friends of this object (up to Nfriensd)
+  // friends are already sorted by distance, so we add closest friends first
+  for (j = 0; (object[Nobj].Nfriend < Nfriends) && (j < object[Nobj].Ndist); j++) {
+    if (object[Nobj].dist[j] > scale) continue;
+
+    int Nnew = object[Nobj].friend[j];
+
+    if (object[Nnew].group != -1) continue;
+
+    // found a friend : 
+    object[Nobj].Nfriend ++;
+    add_to_group (Nnew, scale, Nfriends);
+  }
+  return;
+}
+
+void sortvecset (opihi_flt *X, opihi_int *IDX1, opihi_int *IDX2, int N) {
+
+# define SWAPFUNC(A,B){ opihi_flt tmp; opihi_int itmp; \
+  tmp  = X[A];    X[A]    = X[B];    X[B]    = tmp; \
+  itmp = IDX1[A]; IDX1[A] = IDX1[B]; IDX1[B] = itmp; \
+  itmp = IDX2[A]; IDX2[A] = IDX2[B]; IDX2[B] = itmp; \
+}
+
+# define COMPARE(A,B)(X[A] < X[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+void sortfriends (float *X, int *IDX1, int N) {
+
+# define SWAPFUNC(A,B){ float tmp; int itmp; \
+  tmp  = X[A];    X[A]    = X[B];    X[B]    = tmp; \
+  itmp = IDX1[A]; IDX1[A] = IDX1[B]; IDX1[B] = itmp; \
+}
+
+# define COMPARE(A,B)(X[A] < X[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+int mkclusters (int argc, char **argv) {
+  
+  int i, j;
+  Vector *index1, *index2, *distance;
+
+  if (argc != 6) goto usage;
+
+  if ((index1   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((index2   = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((distance = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
+  float scale = atof (argv[4]);
+  int Nfriends = atoi (argv[5]);
+
+  // XXX enforce matching lengths on the three vectors
+
+  CastVector (index1, OPIHI_INT);
+  CastVector (index2, OPIHI_INT);
+
+  Nobject = 0;
+  for (i = 0; i < index1->Nelements; i++) {
+    Nobject = MAX (Nobject, index1->elements.Int[i]);
+    Nobject = MAX (Nobject, index2->elements.Int[i]);
+  }
+  Nobject ++;  // after the loop, Nobject has the value of the highest index, not the count
+
+  // allocate the list of Object -- these list the possible friends of the given object[i] 
+  ALLOCATE (object, Object, Nobject);
+  for (i = 0; i < Nobject; i++) {
+    ALLOCATE (object[i].friend, int, Nobject);
+    memset (object[i].friend, 0, Nobject*sizeof(int));
+    ALLOCATE (object[i].dist, float, Nobject);
+    memset (object[i].dist, 0, Nobject*sizeof(float));
+    object[i].NDIST = Nobject;
+    object[i].Ndist = 0;
+    object[i].group = -1; // unassigned
+    object[i].Nfriend = 0; // no friends yet
+  }
+
+  // generate the set of vectors of all distances for each entry
+  // the object seq number (object[i]) will be used to match the index1 and index2 values
+  for (i = 0; i < index1->Nelements; i++) {
+    int N, m;
+    N = index1->elements.Int[i];
+    if (N >= Nobject) abort();
+    m = object[N].Ndist;
+    object[N].friend[m] = index2->elements.Int[i];
+    object[N].dist[m] = distance->elements.Flt[i];
+    object[N].Ndist ++;
+    if (object[N].Ndist == object[N].NDIST) {
+      object[N].NDIST += 100;
+      REALLOCATE (object[N].dist, float, object[N].NDIST);
+      REALLOCATE (object[N].friend, int, object[N].NDIST);
+    }
+
+    N = index2->elements.Int[i];
+    if (N >= Nobject) abort();
+    m = object[N].Ndist;
+    object[N].friend[m] = index1->elements.Int[i];
+    object[N].dist[m] = distance->elements.Flt[i];
+    object[N].Ndist ++;
+    if (object[N].Ndist == object[N].NDIST) {
+      object[N].NDIST += 100;
+      REALLOCATE (object[N].dist, float, object[N].NDIST);
+      REALLOCATE (object[N].friend, int, object[N].NDIST);
+    }
+  }
+
+  for (i = 0; i < Nobject; i++) {
+    sortfriends (object[i].dist, object[i].friend, object[i].Ndist);
+  }
+
+  // generate the set of groups. each group is a list of friends and their friends
+  Ngroup = 0;
+  NGROUP = Nobject;
+  ALLOCATE (group, Group, NGROUP);
+  for (i = 0; i < NGROUP; i++) {
+    group[i].Nentry = 0;
+    group[i].NENTRY = 100;
+    ALLOCATE (group[i].entry, int, group[i].NENTRY);
+  }
+
+  sortvecset (distance->elements.Flt, index1->elements.Int, index2->elements.Int, index2->Nelements);
+  
+  for (i = 0; i < distance->Nelements; i++) {
+    int Nnew;
+    int found;
+
+    if (distance->elements.Flt[i] > scale) continue;
+
+    found = FALSE;
+
+    Nnew = index1->elements.Int[i];
+    if (object[Nnew].group == -1) {
+      add_to_group (Nnew, scale, Nfriends);
+      // fprintf (stderr, "new group %d part 1 with %d elements\n", Ngroup, group[Ngroup].Nentry);
+      found = TRUE;
+    }
+
+    Nnew = index2->elements.Int[i];
+    if (object[Nnew].group == -1) {
+      add_to_group (Nnew, scale, Nfriends);
+      // fprintf (stderr, "new group %d part 1 with %d elements\n", Ngroup, group[Ngroup].Nentry);
+      found = TRUE;
+    }
+
+    if (found) {
+      fprintf (stderr, "group %d with %d elements\n", Ngroup, group[Ngroup].Nentry);
+      Ngroup ++;
+      if (Ngroup >= NGROUP) abort();
+    }
+  }
+
+  if (1) {
+    int Nmiss = 0;
+    for (i = 0; i < Nobject; i++) {
+      if (object[i].group > -1) continue;
+      Nmiss ++;
+      // fprintf (stderr, "object %d not assigned\n", i);
+      for (j = 0; FALSE && j < object[i].Ndist; j++) {
+	fprintf (stderr, "friend %d, dist %f\n", object[i].friend[j], object[i].dist[j]);
+      }
+    }
+    fprintf (stderr, "%d objects not assigned\n", Nmiss);
+  } 
+
+  return TRUE;
+  
+ escape: 
+  gprint (GP_ERR, "invalid vector\n");
+  return FALSE;
+  
+usage:
+  gprint (GP_ERR, "USAGE: mkclusters (index1) (index2) (distance) scale\n");
+  return FALSE;
+}
+
+/* this function takes 3 vectors: a set of distances and a pair of indicies identifying the end points.  
+   
+   the goal is to generate groups of the end points that are relatively closer.
+
+   note that the distances may not be Euclidean: nothing can be assumed about relationships between d(1,2), d(2,3), and d(1,3)
+
+   the indicies should be sequential.  the algorithm does not require it, but storage is much cleaner if it is
+
+*/
Index: trunk/Ohana/src/opihi/cmd.astro/specpairfit.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/specpairfit.c	(revision 33662)
+++ trunk/Ohana/src/opihi/cmd.astro/specpairfit.c	(revision 33662)
@@ -0,0 +1,61 @@
+# include "astro.h"
+
+int specpairfit (int argc, char **argv) {
+  
+  int i;
+  Vector *flux1, *flux2, *dflux1, *dflux2, *window;
+
+  if (argc != 7) goto usage;
+
+  if ((flux1  = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((dflux1 = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((flux2  = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((dflux2 = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((window = SelectVector (argv[5], OLDVECTOR, TRUE)) == NULL) goto escape;
+  
+  // XXX enforce matching lengths on the 6 vectors
+
+  int mask = atoi (argv[6]);
+  CastVector (window, OPIHI_INT);
+
+  // minimize (flux1 - flux2*alpha) in window defined by mask
+  // note that the mask is a SELECTION mask not an EXCLUSION mask
+
+  double F12 = 0.0; 
+  double F22 = 0.0;
+  for (i = 0; i < flux1->Nelements; i++) {
+    if ((mask & window->elements.Int[i]) == 0) continue;
+    double weight = 1.0 / (SQ(dflux1->elements.Flt[i]) + SQ(dflux2->elements.Flt[i]));
+    F12 += flux1->elements.Flt[i] * flux2->elements.Flt[i] * weight;
+    F22 += flux2->elements.Flt[i] * flux2->elements.Flt[i] * weight;
+  }
+
+  double Ao = F12 / F22;
+  double dA = sqrt(1.0 / F22);
+
+  int Ndof = -1; // 1 parameter fit
+  double chisq = 0.0;
+  for (i = 0; i < flux1->Nelements; i++) {
+    if ((mask & window->elements.Int[i]) == 0) continue;
+    double weight = 1.0 / (SQ(dflux1->elements.Flt[i]) + SQ(dflux2->elements.Flt[i]));
+    chisq += SQ(flux1->elements.Flt[i] - Ao * flux2->elements.Flt[i]) * weight;
+    Ndof ++;
+  }
+
+  double chisqNu = chisq / Ndof;
+
+  // fprintf (stderr, "Ao: %f +/- %f, chisq: %f, chisq_nu : %f for %d dof\n", Ao, dA, chisq, chisqNu, Ndof);
+  set_variable ("Ao", Ao);
+  set_variable ("dA", dA);
+  set_variable ("Xv", chisqNu);
+  set_variable ("Nd", Ndof);
+  return (TRUE);
+  
+ escape: 
+  gprint (GP_ERR, "invalid vector\n");
+  return (FALSE);
+  
+usage:
+  gprint (GP_ERR, "USAGE: specpairfit (flux1) (dflux1) (flux2) (dflux2) (window) (mask) [options]\n");
+  return FALSE;
+}
Index: trunk/Ohana/src/opihi/cmd.astro/spex1dgas.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/spex1dgas.c	(revision 33662)
+++ trunk/Ohana/src/opihi/cmd.astro/spex1dgas.c	(revision 33662)
@@ -0,0 +1,438 @@
+# include "astro.h"
+
+/* We have N objects with N(N-1)/2 paired distances.  We want to find the 1D distribution that
+   best describes the observed distances.  Assume the distances are in a 1D space, though in
+   principle it could be N-D or even non-euclidean.
+
+   Model the points as a gas under pressure.  
+
+   Start with the a guess for positions of the objects at some locations in 2D
+
+   P = (d_now - d_tru)
+
+   Iterate over the points and move in X,Y
+
+   Find (dP/dX,dP/dY) for modest moves in X,Y
+
+   Move based on dP
+
+   USAGE: spex1dgas incdex1 index2 distance
+
+   outline:
+
+   * load data
+   * generate the unique objects
+   * determine the max distance needed
+   * place each object in the 2D space
+
+   * iterate:
+   ** calculate dP/dX,dPdY for each object
+   ** move each object proportionally to the pressure gradient
+
+   */
+
+typedef struct {
+  int Nindex;   // number of tested relationships
+  int NINDEX;   // number of allocated relationships
+  int *index;  // indices for all other objects
+  float *Dtgt; // target distance for this relationship
+  float *Dcur; // current distance for this relationship
+  float Xo,Yo; // current X,Y position of this object
+  float dPdX;  // pressure in X
+  float dPdY;  // pressure in Y
+} Object;
+
+static Object *object = NULL;
+static int Nobject = 0;
+// static int NOBJECT = 0;
+
+void sortfriends (float *X, int *IDX1, int N);
+
+static void get_pressure_gradient (int iObj, int nCloseMax, float farFrac, int nearNeighbors, float maxPressure) {
+
+  int i;
+  
+  float dPdX = 0.0;
+
+  // only use the first N friends
+  for (i = 0; (i < object[iObj].Nindex); i++) {
+    int jObj = object[iObj].index[i];
+
+    // some options:
+    // if the iterations are small, we should only worry about getting the near neighbors right
+    // if the iterations are large, we should add in more distant objects
+    if (nearNeighbors && (i >= nCloseMax)) break;
+    if (!nearNeighbors && (i >= nCloseMax)) {
+      if (drand48() > farFrac) continue;
+    }
+
+    float Dtgt = object[iObj].Dtgt[i];
+    float dX = object[jObj].Xo - object[iObj].Xo;
+    float Dcur = dX;
+
+    // the force law as a function of (Dcur - Dtgt) : if Dcur is too large, dF is negative
+    // float dF = (Dcur < 0.01*Dtgt) ? -100.0 : (Dcur - Dtgt) / Dcur; XXX modified spring constant : too crazy
+    float dF = (Dcur - Dtgt);
+    dF = MIN (maxPressure, MAX (-maxPressure, dF));
+
+    float dPdXi;
+    if (fabs(Dcur) < 1e-6) {
+      dPdXi = 0.0;
+    } else {
+      dPdXi = dF;
+    }
+
+    if (isnan(Dtgt) || isnan(dX) || isnan(Dcur) || isnan(dF) || isnan(dPdXi)) abort();
+    dPdX += dPdXi;
+  }
+  object[iObj].dPdX = dPdX;
+  object[iObj].dPdY = 0.0;
+
+  return;
+}
+
+static void move_object (int iObj) {
+
+  object[iObj].Xo += 0.25*object[iObj].dPdX;
+  return;
+}
+
+int spex1dgas (int argc, char **argv) {
+  
+  int i, j, iter;
+  Vector *index1, *index2, *distance;
+  float *XoList;
+  int *IDList, *MidObj;
+
+  MidObj = NULL;
+  IDList = NULL;
+  XoList = NULL;
+  { 
+    // init random numbers
+    long A, B;
+    A = time(NULL);
+    for (B = 0; A == time(NULL); B++);
+    srand48(B);
+  }
+ 
+  if (argc != 11) goto usage;
+
+  // XXX enforce matching lengths on the three vectors
+  if ((index1   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((index2   = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((distance = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
+
+  int Niter = atoi (argv[4]);
+  int nCloseMax = atoi (argv[5]);
+  int nCloseIter = atoi (argv[6]);
+  float farFrac = atof (argv[7]);
+  float maxPressure = atof (argv[8]);
+
+  int idx1 = atoi(argv[9]);
+  int idx2 = atoi(argv[10]);
+
+  CastVector (index1, OPIHI_INT);
+  CastVector (index2, OPIHI_INT);
+
+  // how many objects do we have?
+  Nobject = 0;
+  for (i = 0; i < index1->Nelements; i++) {
+    Nobject = MAX (Nobject, index1->elements.Int[i]);
+    Nobject = MAX (Nobject, index2->elements.Int[i]);
+  }
+  Nobject ++;  // after the loop, Nobject has the value of the highest index, not the count
+
+  // allocate the list of Object -- these list the possible friends of the given object[i] 
+  ALLOCATE (object, Object, Nobject);
+  for (i = 0; i < Nobject; i++) {
+    ALLOCATE (object[i].index, int, Nobject);
+    memset (object[i].index, 0, Nobject*sizeof(int));
+    ALLOCATE (object[i].Dtgt, float, Nobject);
+    memset (object[i].Dtgt, 0, Nobject*sizeof(float));
+    ALLOCATE (object[i].Dcur, float, Nobject);
+    memset (object[i].Dcur, 0, Nobject*sizeof(float));
+    object[i].NINDEX = Nobject;
+    object[i].Nindex = 0;
+  }
+
+  // generate the set of vectors of all distances for each entry
+  // the object seq number (object[i]) will be used to match the index1 and index2 values
+  for (i = 0; i < index1->Nelements; i++) {
+    int N, m;
+
+    N = index1->elements.Int[i];
+    if (N >= Nobject) abort();
+    m = object[N].Nindex;
+    object[N].index[m] = index2->elements.Int[i];
+    object[N].Dtgt[m] = distance->elements.Flt[i];
+    object[N].Dcur[m] = 0.0;
+    object[N].Nindex ++;
+    if (object[N].Nindex == object[N].NINDEX) {
+      object[N].NINDEX += 100;
+      REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
+      REALLOCATE (object[N].Dcur, float, object[N].NINDEX);
+      REALLOCATE (object[N].index, int, object[N].NINDEX);
+    }
+
+    // Above we used the set (N,M,D_N,M) to assigned the distance of obj M to obj N.
+    // Below we use that set to assign the distance of obj N to obj M.  IF we supply all
+    // pairwise distances, then this duplicates all entries, so we need to skip that step
+    // in such a case.
+
+    // XX N = index2->elements.Int[i];
+    // XX if (N >= Nobject) abort();
+    // XX m = object[N].Nindex;
+    // XX object[N].index[m] = index1->elements.Int[i];
+    // XX object[N].Dtgt[m] = distance->elements.Flt[i];
+    // XX object[N].Dcur[m] = 0.0;
+    // XX object[N].Nindex ++;
+    // XX if (object[N].Nindex == object[N].NINDEX) {
+    // XX   object[N].NINDEX += 100;
+    // XX   REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
+    // XX   REALLOCATE (object[N].Dcur, float, object[N].NINDEX);
+    // XX   REALLOCATE (object[N].index, int, object[N].NINDEX);
+    // XX }
+  }
+
+  for (i = 0; i < Nobject; i++) {
+    // sort so closest friends are first
+    sortfriends (object[i].Dtgt, object[i].index, object[i].Nindex);
+  }
+
+  // find and save the max distance
+  float Dmax = 0;
+  for (i = 0; i < distance->Nelements; i++) {
+    Dmax = MAX (Dmax, distance->elements.Flt[i]);
+  }
+
+  // place the objects at the initial guess locations
+  // XXX let's try with just a simple grid dividing up the max range
+  // int Ngrid = sqrt(Nobject);
+  // float dgrid = 1.5 * Dmax / Ngrid; // XXX remove the fudge factor
+  // XXX for (i = 0; i < Nobject; i++) {
+  // XXX   // object[i].Xo = dgrid * (int) (i % Ngrid);
+  // XXX   // object[i].Yo = dgrid * (int) (i / Ngrid);
+  // XXX   object[i].Xo = Dmax*drand48();
+  // XXX   object[i].Yo = Dmax*drand48();
+  // XXX   object[i].dPdX = 0.0;
+  // XXX   object[i].dPdY = 0.0;
+  // XXX }
+
+  // place the objects at the initial guess locations
+  // use object 0 and its most distant friend to constrain:
+
+  // we are going to use the projected X distance as the location in X, and ignore the Y
+  // coordinate
+  
+  // need to get pretty close to a good starting point.  start with a guess using the first entry & its extremum
+  // int idx1 = object[0].index[(int)(0.75*object[0].Nindex)];
+  // int idx2 = object[idx1].index[(int)(0.75*object[idx1].Nindex)];
+  float A = NAN;
+  for (i = 0; isnan (A) && (i < object[idx2].Nindex); i++) {
+      if (object[idx2].index[i] == idx1) {
+	  A = object[idx2].Dtgt[i];
+      }
+  }
+  fprintf (stderr, "idx1: %d, idx2: %d, A: %f\n", idx1, idx2, A);
+  set_variable ("SP1D_A", A);
+  object[idx1].Xo = 0.0;
+  object[idx1].Yo = 0.0;
+  object[idx1].dPdX = 0.0;
+  object[idx1].dPdY = 0.0;
+
+  object[idx2].Xo = A;
+  object[idx2].Yo = 0.0;
+  object[idx2].dPdX = 0.0;
+  object[idx2].dPdY = 0.0;
+		 
+  Vector *tmpB = SelectVector ("sp1d_B", ANYVECTOR, TRUE); 
+  Vector *tmpC = SelectVector ("sp1d_C", ANYVECTOR, TRUE); 
+  ResetVector (tmpB, OPIHI_FLT, Nobject);
+  ResetVector (tmpC, OPIHI_FLT, Nobject);
+
+  for (i = 0; i < Nobject; i++) {
+    if (i == idx1) continue;
+    if (i == idx2) continue;
+    float B = NAN;
+    float C = NAN;
+    for (j = 0; (isnan(B) || isnan(C)) && (j < object[i].Nindex); j++) {
+      if (object[i].index[j] == idx1) { B = object[i].Dtgt[j]; }
+      if (object[i].index[j] == idx2) { C = object[i].Dtgt[j]; }
+    }
+    if (isnan(B) || isnan(C)) abort();
+
+    // C^2 = A^2 + B^2 - 2AB cos(t)
+    // Xo = B cos(t)
+    // Xo = A^2 + B^2 - C^2 / 2 A
+    float Xo = (SQ(A) + SQ(B) - SQ(C)) / (2*A);
+
+    object[i].Xo = Xo;
+    object[i].dPdX = 0.0;
+    object[i].Yo = 0.0;
+    object[i].dPdY = 0.0;
+    tmpB->elements.Flt[i] = B;
+    tmpC->elements.Flt[i] = C;
+  }
+  tmpB->elements.Flt[idx1] = 0;
+  tmpB->elements.Flt[idx2] = A;
+  tmpC->elements.Flt[idx1] = A;
+  tmpC->elements.Flt[idx2] = 0;
+
+  if (Niter >= -1) {
+    // now choose several in the middle of the range, and find the mean distance for each relative to that group
+    ALLOCATE (XoList, float, Nobject);
+    ALLOCATE (IDList, int, Nobject);
+    for (i = 0; i < Nobject; i++) {
+      XoList[i] = object[i].Xo;
+      IDList[i] = i;
+    }
+    sortfriends (XoList, IDList, Nobject);
+    float XoMid = XoList[(int)(Nobject/2)];
+
+# define NMID 5
+    ALLOCATE (MidObj, int, NMID);
+    MidObj[0] = IDList[(int)(Nobject/2) + 0];
+    MidObj[1] = IDList[(int)(Nobject/2) + 1];
+    MidObj[2] = IDList[(int)(Nobject/2) - 1];
+    MidObj[3] = IDList[(int)(Nobject/2) + 2];
+    MidObj[4] = IDList[(int)(Nobject/2) - 2];
+
+    for (i = 0; i < Nobject; i++) {
+      for (j = 0; j < NMID; j++) {
+	if (MidObj[j] == i) goto skip_object;
+      }
+      float Dmean = 0.0;
+      for (j = 0; j < NMID; j++) {
+	Dmean += object[i].Dtgt[MidObj[j]]; // desired distance between object (i) and object in middle section
+      }
+      Dmean = Dmean / NMID;
+      if (object[i].Xo < XoMid) {
+	object[i].Xo = XoMid - Dmean;
+      } else {
+	object[i].Xo = XoMid + Dmean;
+      }
+    skip_object:
+      continue;
+    }
+  }
+
+  if (Niter >= 0) {
+    float XoMean = 0.0;
+    // now choose one end of the range, and find the mean distance for each relative to THAT group
+    for (i = 0; i < Nobject; i++) {
+      XoList[i] = object[i].Xo;
+      IDList[i] = i;
+      XoMean += object[i].Xo;
+    }
+    XoMean = XoMean / Nobject;
+    sortfriends (XoList, IDList, Nobject);
+    float XoMin = XoList[0];
+
+    MidObj[0] = IDList[0];
+    MidObj[1] = IDList[1];
+    MidObj[2] = IDList[2];
+    MidObj[3] = IDList[3];
+    MidObj[4] = IDList[4];
+
+    for (i = 0; i < Nobject; i++) {
+      for (j = 0; j < NMID; j++) {
+	if (MidObj[j] == i) goto skip_object_p2;
+      }
+      float Dmean = 0.0;
+      for (j = 0; j < NMID; j++) {
+	Dmean += object[i].Dtgt[MidObj[j]]; // desired distance between object (i) and object in middle section
+      }
+      Dmean = Dmean / NMID;
+      if (object[i].Xo < XoMin) {
+	object[i].Xo = XoMin - Dmean;
+      } else {
+	object[i].Xo = XoMin + Dmean;
+      }
+    skip_object_p2:
+      continue;
+    }
+  }
+
+  for (iter = 0; iter < Niter; iter ++) {
+    fprintf (stderr, "iter %d\n", iter);
+
+    int nearNeighbors = (iter < nCloseIter);
+
+    // save the result
+    if (0) {
+      char name[64];
+      snprintf (name, 64, "output.%02d.dat", iter);
+      FILE *output = fopen (name, "w");
+      for (i = 0; i < Nobject; i++) {
+	fprintf (output, "%f %f : %f %f\n", object[i].Xo, object[i].Yo, object[i].dPdX, object[i].dPdY);
+      }
+      fclose (output);
+    }
+
+    // measure (dP/dX) & move object
+    for (i = 0; i < Nobject; i++) {
+      get_pressure_gradient (i, nCloseMax, farFrac, nearNeighbors, maxPressure);
+      move_object (i);
+    }
+  }
+
+  // find the mean Xo, Yo positions for all objects
+  float XoSum = 0.0;
+  for (i = 0; i < Nobject; i++) {
+    float dY2sum = 0.0;
+    int Npts = 0;
+    for (j = 0; (j < nCloseMax) && (j < object[i].Nindex); j++) {
+      int idx = object[i].index[j];
+      float dR2 = SQ(object[i].Dtgt[j]);
+      float dX2 = SQ(object[i].Xo - object[idx].Xo);
+      float dY2 = dR2 - dX2;
+      dY2sum += dY2;
+      Npts ++;
+    }
+    XoSum += object[i].Xo;
+    object[i].Yo = (dY2sum > 0.0) ? sqrt(fabs(dY2sum / Npts)) : -sqrt(fabs(dY2sum / Npts));
+  }
+  float XoMean = XoSum / Nobject;
+
+  // remove the mean Xo (uninteresting)
+  for (i = 0; i < Nobject; i++) {
+    object[i].Xo -= XoMean;
+  }
+
+  // save the result
+  {
+    Vector *outindex = SelectVector ("sp1d_idx", ANYVECTOR, TRUE); if (!outindex) goto escape;
+    Vector *outXo    = SelectVector ("sp1d_Xo",  ANYVECTOR, TRUE); if (!outXo) goto escape;
+    Vector *outYo    = SelectVector ("sp1d_Yo",  ANYVECTOR, TRUE); if (!outYo) goto escape;
+
+    ResetVector (outindex, OPIHI_INT, Nobject);
+    ResetVector (outXo, OPIHI_FLT, Nobject);
+    ResetVector (outYo, OPIHI_FLT, Nobject);
+
+    for (i = 0; i < Nobject; i++) {
+      outindex->elements.Int[i] = i;
+      outXo->elements.Flt[i] = object[i].Xo;
+      outYo->elements.Flt[i] = object[i].Yo;
+    }
+  }
+
+  return TRUE;
+  
+escape: 
+  gprint (GP_ERR, "invalid vector\n");
+  return FALSE;
+  
+usage:
+  gprint (GP_ERR, "USAGE: spex2dgas (index1) (index2) (distance) (Niter) (nCloseMax) (nCloseIter) (farFrac) (maxPressure) (idx1) (idx2)\n");
+  return FALSE;
+}
+
+/* this function takes 3 vectors: a set of distances and a pair of indicies identifying the end points.  
+   
+   the goal is to generate groups of the end points that are relatively closer.
+
+   note that the distances may not be Euclidean: nothing can be assumed about relationships between d(1,2), d(2,3), and d(1,3)
+
+   the indicies should be sequential.  the algorithm does not require it, but storage is much cleaner if it is
+
+*/
Index: trunk/Ohana/src/opihi/cmd.astro/spex2dgas.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/spex2dgas.c	(revision 33662)
+++ trunk/Ohana/src/opihi/cmd.astro/spex2dgas.c	(revision 33662)
@@ -0,0 +1,344 @@
+# include "astro.h"
+
+/* We have N objects with N(N-1)/2 paired distances.  We want to find the 2D distribution that
+   best describes the observed distances.  Assume the distances are in a 2D space, though in
+   principle it could be N-D or even non-euclidean.
+
+   Model the points as a gas under pressure.  
+
+   Start with the a guess for positions of the objects at some locations in 2D
+
+   P = (d_now - d_tru)
+
+   Iterate over the points and move in X,Y
+
+   Find (dP/dX,dP/dY) for modest moves in X,Y
+
+   Move based on dP
+
+   USAGE: spex2dgas incdex1 index2 distance
+
+   outline:
+
+   * load data
+   * generate the unique objects
+   * determine the max distance needed
+   * place each object in the 2D space
+
+   * iterate:
+   ** calculate dP/dX,dPdY for each object
+   ** move each object proportionally to the pressure gradient
+
+   */
+
+typedef struct {
+  int Nindex;   // number of tested relationships
+  int NINDEX;   // number of allocated relationships
+  int *index;  // indices for all other objects
+  float *Dtgt; // target distance for this relationship
+  float *Dcur; // current distance for this relationship
+  float Xo,Yo; // current X,Y position of this object
+  float dPdX;  // pressure in X
+  float dPdY;  // pressure in Y
+} Object;
+
+static Object *object = NULL;
+static int Nobject = 0;
+// static int NOBJECT = 0;
+
+void sortfriends (float *X, int *IDX1, int N);
+
+static void get_pressure_gradient (int iObj, int nCloseMax, float farFrac, int nearNeighbors, float maxPressure) {
+
+  int i;
+  
+  float dPdX = 0.0;
+  float dPdY = 0.0;
+
+  // only use the first N friends
+  for (i = 0; (i < object[iObj].Nindex); i++) {
+    int jObj = object[iObj].index[i];
+
+    // some options:
+    // if the iterations are small, we should only worry about getting the near neighbors right
+    // if the iterations are large, we should add in more distant objects
+    if (nearNeighbors && (i >= nCloseMax)) break;
+    if (!nearNeighbors && (i >= nCloseMax)) {
+      if (drand48() > farFrac) continue;
+    }
+
+    float Dtgt = object[iObj].Dtgt[i];
+    float dX = object[jObj].Xo - object[iObj].Xo;
+    float dY = object[jObj].Yo - object[iObj].Yo;
+    float Dcur = hypot(dX,dY);
+
+    // the force law as a function of (Dcur - Dtgt) : if Dcur is too large, dF is negative
+    // float dF = (Dcur < 0.01*Dtgt) ? -100.0 : (Dcur - Dtgt) / Dcur; XXX modified spring constant : too crazy
+    float dF = (Dcur - Dtgt);
+    dF = MIN (maxPressure, MAX (-maxPressure, dF));
+
+    float dPdXi, dPdYi;
+    if (fabs(Dcur) < 1e-6) {
+      dPdXi = 0.0;
+      dPdYi = 0.0;
+    } else {
+      dPdXi = dF * dX / Dcur;
+      dPdYi = dF * dY / Dcur;
+    }
+
+    if (isnan(Dtgt) || isnan(dX) || isnan(dY) || isnan(Dcur) || isnan(dF) || isnan(dPdXi) || isnan(dPdYi)) abort();
+
+    // if we are too close, then dX/Dcur is too ill-defined, just jump away by 5% of Dtgt
+    if (Dcur < 0.01*Dtgt) {
+      dPdXi = (drand48() - 0.5)*0.1*Dtgt;
+      dPdYi = (drand48() - 0.5)*0.1*Dtgt;
+    }
+
+    if (i >= nCloseMax) {
+      fprintf (stderr, "Dcur,Dtgt : %f %f : dX,dY,dP : %f %f : %f : %f %f\n", Dcur, Dtgt, dX, dY, dF, dPdXi, dPdYi);
+    }
+    dPdX += dPdXi;
+    dPdY += dPdYi;
+  }
+  object[iObj].dPdX = dPdX;
+  object[iObj].dPdY = dPdY;
+
+  return;
+}
+
+static void move_object (int iObj) {
+
+  object[iObj].Xo += 0.25*object[iObj].dPdX;
+  object[iObj].Yo += 0.25*object[iObj].dPdY;
+  return;
+}
+
+int spex2dgas (int argc, char **argv) {
+  
+  int i, j, iter, IoMax;
+  Vector *index1, *index2, *distance;
+  float XoMax, YoMax;
+
+  { 
+    // init random numbers
+    long A, B;
+    A = time(NULL);
+    for (B = 0; A == time(NULL); B++);
+    srand48(B);
+  }
+ 
+  if (argc != 9) goto usage;
+
+  if ((index1   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((index2   = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((distance = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
+  int Niter = atoi (argv[4]);
+  int nCloseMax = atoi (argv[5]);
+  int nCloseIter = atoi (argv[6]);
+  float farFrac = atof (argv[7]);
+  float maxPressure = atof (argv[8]);
+  // XXX enforce matching lengths on the three vectors
+
+  CastVector (index1, OPIHI_INT);
+  CastVector (index2, OPIHI_INT);
+
+  // how many objects do we have?
+  Nobject = 0;
+  for (i = 0; i < index1->Nelements; i++) {
+    Nobject = MAX (Nobject, index1->elements.Int[i]);
+    Nobject = MAX (Nobject, index2->elements.Int[i]);
+  }
+  Nobject ++;  // after the loop, Nobject has the value of the highest index, not the count
+
+  // allocate the list of Object -- these list the possible friends of the given object[i] 
+  ALLOCATE (object, Object, Nobject);
+  for (i = 0; i < Nobject; i++) {
+    ALLOCATE (object[i].index, int, Nobject);
+    memset (object[i].index, 0, Nobject*sizeof(int));
+    ALLOCATE (object[i].Dtgt, float, Nobject);
+    memset (object[i].Dtgt, 0, Nobject*sizeof(float));
+    ALLOCATE (object[i].Dcur, float, Nobject);
+    memset (object[i].Dcur, 0, Nobject*sizeof(float));
+    object[i].NINDEX = Nobject;
+    object[i].Nindex = 0;
+  }
+
+  // generate the set of vectors of all distances for each entry
+  // the object seq number (object[i]) will be used to match the index1 and index2 values
+  for (i = 0; i < index1->Nelements; i++) {
+    int N, m;
+
+    N = index1->elements.Int[i];
+    if (N >= Nobject) abort();
+    m = object[N].Nindex;
+    object[N].index[m] = index2->elements.Int[i];
+    object[N].Dtgt[m] = distance->elements.Flt[i];
+    object[N].Dcur[m] = 0.0;
+    object[N].Nindex ++;
+    if (object[N].Nindex == object[N].NINDEX) {
+      object[N].NINDEX += 100;
+      REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
+      REALLOCATE (object[N].Dcur, float, object[N].NINDEX);
+      REALLOCATE (object[N].index, int, object[N].NINDEX);
+    }
+
+    N = index2->elements.Int[i];
+    if (N >= Nobject) abort();
+    m = object[N].Nindex;
+    object[N].index[m] = index1->elements.Int[i];
+    object[N].Dtgt[m] = distance->elements.Flt[i];
+    object[N].Dcur[m] = 0.0;
+    object[N].Nindex ++;
+    if (object[N].Nindex == object[N].NINDEX) {
+      object[N].NINDEX += 100;
+      REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
+      REALLOCATE (object[N].Dcur, float, object[N].NINDEX);
+      REALLOCATE (object[N].index, int, object[N].NINDEX);
+    }
+  }
+
+  for (i = 0; i < Nobject; i++) {
+    // sort so closest friends are first
+    sortfriends (object[i].Dtgt, object[i].index, object[i].Nindex);
+  }
+
+  // find and save the max distance
+  float Dmax = 0;
+  for (i = 0; i < distance->Nelements; i++) {
+    Dmax = MAX (Dmax, distance->elements.Flt[i]);
+  }
+
+  // place the objects at the initial guess locations
+  // XXX let's try with just a simple grid dividing up the max range
+  // int Ngrid = sqrt(Nobject);
+  // float dgrid = 1.5 * Dmax / Ngrid; // XXX remove the fudge factor
+  // XXX for (i = 0; i < Nobject; i++) {
+  // XXX   // object[i].Xo = dgrid * (int) (i % Ngrid);
+  // XXX   // object[i].Yo = dgrid * (int) (i / Ngrid);
+  // XXX   object[i].Xo = Dmax*drand48();
+  // XXX   object[i].Yo = Dmax*drand48();
+  // XXX   object[i].dPdX = 0.0;
+  // XXX   object[i].dPdY = 0.0;
+  // XXX }
+
+  // place the objects at the initial guess locations
+  // use object 0 and its most distant friend to constrain:
+  
+  int idx1 = 0;
+  int idx2 = object[idx1].index[object[idx1].Nindex-1];
+  float A = object[idx1].Dtgt[object[idx1].Nindex-1];
+  object[idx1].Xo = 0.0;
+  object[idx1].Yo = 0.0;
+  object[idx1].dPdX = 0.0;
+  object[idx1].dPdY = 0.0;
+
+  object[idx2].Xo = A;
+  object[idx2].Yo = 0.0;
+  object[idx2].dPdX = 0.0;
+  object[idx2].dPdY = 0.0;
+		 
+  // choose the correct Y side by comparing to the distance from the most deviant
+  YoMax = 0;
+  XoMax = 0;
+  IoMax = 0;
+
+  for (i = 0; i < Nobject; i++) {
+    if (i == idx1) continue;
+    if (i == idx2) continue;
+    float B = NAN;
+    float C = NAN;
+    for (j = 0; (isnan(B) || isnan(C)) && (j < object[i].Nindex); j++) {
+      if (object[i].index[j] == idx1) { B = object[i].Dtgt[j]; }
+      if (object[i].index[j] == idx2) { C = object[i].Dtgt[j]; }
+    }
+    if (isnan(B) || isnan(C)) abort();
+
+    float Xo = (SQ(A) + SQ(B) - SQ(C)) / (2*A);
+    float Y2 = SQ(B) - SQ(Xo);
+
+    float Yo = (Y2 < 0) ? 0.0 : sqrt(Y2);
+    if (isnan(Yo)) abort();
+
+    object[i].Xo = Xo;
+    object[i].Yo = Yo;
+    object[i].dPdX = 0.0;
+    object[i].dPdY = 0.0;
+    if (Yo > YoMax) {
+      YoMax = Yo;
+      XoMax = Xo;
+      IoMax = i;
+    }
+  }
+
+  for (i = 0; i < Nobject; i++) {
+    if (i == idx1) continue;
+    if (i == idx2) continue;
+    if (i == IoMax) continue;
+    for (j = 0; j < object[i].Nindex; j++) {
+      if (object[i].index[j] == IoMax) { 
+	float d02 = SQ(object[i].Dtgt[j]);
+	float dX  = XoMax - object[i].Xo;
+	float dY1 = YoMax - object[i].Yo;
+	float dY2 = YoMax + object[i].Yo;
+	float d12 = SQ(dX) + SQ(dY1);
+	float d22 = SQ(dX) + SQ(dY2);
+	if (fabs(d12 - d02) > fabs(d22 - d02)) {
+	  object[i].Yo = -object[i].Yo;
+	  break;
+	}
+      }
+    }
+  }
+
+  for (iter = 0; iter < Niter; iter ++) {
+    fprintf (stderr, "iter %d\n", iter);
+
+    int nearNeighbors = (iter < nCloseIter);
+
+    // save the result
+    char name[64];
+    snprintf (name, 64, "output.%02d.dat", iter);
+    FILE *output = fopen (name, "w");
+    for (i = 0; i < Nobject; i++) {
+      fprintf (output, "%f %f : %f %f\n", object[i].Xo, object[i].Yo, object[i].dPdX, object[i].dPdY);
+    }
+    fclose (output);
+
+    // measure (dP/dX),(dP/dY) for all objects
+    for (i = 0; i < Nobject; i++) {
+      get_pressure_gradient (i, nCloseMax, farFrac, nearNeighbors, maxPressure);
+      move_object (i);
+    }
+
+    // given (dP/dX),(dP/dY), move each object 
+    // for (i = 0; i < Nobject; i++) {
+    // }
+  }
+
+  // save the result
+  FILE *output = fopen ("output.dat", "w");
+  for (i = 0; i < Nobject; i++) {
+    fprintf (output, "%f %f : %f %f\n", object[i].Xo, object[i].Yo, object[i].dPdX, object[i].dPdY);
+  }
+  fclose (output);
+
+  return TRUE;
+  
+escape: 
+  gprint (GP_ERR, "invalid vector\n");
+  return FALSE;
+  
+usage:
+  gprint (GP_ERR, "USAGE: spex2dgas (index1) (index2) (distance) (Niter) (nCloseMax) (nCloseIter) (farFrac) (maxPressure)\n");
+  return FALSE;
+}
+
+/* this function takes 3 vectors: a set of distances and a pair of indicies identifying the end points.  
+   
+   the goal is to generate groups of the end points that are relatively closer.
+
+   note that the distances may not be Euclidean: nothing can be assumed about relationships between d(1,2), d(2,3), and d(1,3)
+
+   the indicies should be sequential.  the algorithm does not require it, but storage is much cleaner if it is
+
+*/
Index: trunk/Ohana/src/opihi/cmd.astro/spexseq.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/spexseq.c	(revision 33662)
+++ trunk/Ohana/src/opihi/cmd.astro/spexseq.c	(revision 33662)
@@ -0,0 +1,381 @@
+# include "astro.h"
+
+/* We have N objects with N(N-1)/2 paired chisq.  We want to find the 1D sequence / index that
+   best describes the observed distances.  
+
+   A sequence is best described by the distances from objects which are only modestly close.
+   If we get too close to the reference, we (a) hit a non-euclidean (or N-D) minimum (A-B ~ A-C
+   ~ B-C) and (b) the positive and negative sides of the sequence are ambiguous.  If we get too
+   far away, the distances are all either equivalent or less consistent
+
+   1) find the dynamic range of distances
+
+   2) choose object (1) and select a friend (2) > a fixed fraction of the full dynamic range (say,
+   10-20%)
+
+   3) all nearby friends of (1) (D < f MaxDist) have sequences assigned based on the distance
+   to (2).  center based on the distance (1)-(2) ie, (1) gets a sequence value of 0.0
+
+   4) choose another object (1') and repeat [need to choose (2') on the same side as (2)?]
+
+   ** reconciliation is a bit tricky
+   ** edge effects imply that we prefer to start in the center
+   ** is this well-defined? 
+
+   */
+
+typedef struct {
+  int Nindex;   // number of tested relationships
+  int NINDEX;   // number of allocated relationships
+  int *index;  // index for distances to other objects: object[i].Dtgt[j] is distance from object[i] to object[idx] where idx = object[i].index[j] (sorted by Dtgt)
+  int *rindex;  // reverse index for Dtgt: distance from object[i] to object[idx] is object[i].Dtgt[j] where object[i].rindex[idx] = j
+  int shifted; // has this object already had the sequence adjusted?
+  float *Dtgt; // target distance for this relationship
+  float *seq;  // sequence value for this friend
+  float So, dSo;
+  float Xo, dXo;
+  float Yo, dYo;
+  float Srange;
+  int nSeq;
+} Object;
+
+void sortfriends (float *X, int *IDX1, int N);
+
+void get_sequence (Object *object, int Nobject, int idx1, float Dmax, float f1, float f2, int pin1, int pin2) {
+
+  int i, j, N;
+
+  // object (1) is provided above; choose object (2)
+  int idx2 = -1;
+  float A = NAN;
+  for (i = 0; (idx2 == -1) && (i < object[idx1].Nindex); i++) {
+    if (object[idx1].Dtgt[i] > f1*Dmax) {
+      idx2 = object[idx1].index[i];
+    }
+    A = object[idx1].Dtgt[i];
+  }
+  // fprintf (stderr, "idx1: %d, idx2: %d, A: %f\n", idx1, idx2, A);
+  // set_variable ("SP1D_A", A);
+
+  // is idx2 closer to pin1 or pin2?
+  // XX N = object[idx2].rindex[pin1];
+  // XX float D1 = object[idx2].Dtgt[N];
+  // XX N = object[idx2].rindex[pin2];
+  // XX float D2 = object[idx2].Dtgt[N];
+  // XX int parity = (D1 < D2) ? +1 : -1;
+
+  int Nfriends = 0;
+  for (i = 0; (object[idx1].Dtgt[i] < f2*Dmax) && (i < object[idx1].Nindex); i++) {
+    // what is the distance from this friend to object (2)?
+    N = object[idx1].index[i];
+    j = object[idx2].rindex[N];
+    assert (object[idx2].index[j] == N);
+    // XX object[idx1].seq[i] = parity*(object[idx2].Dtgt[j] - A);
+    object[idx1].seq[i] = (object[idx2].Dtgt[j] - A);
+    Nfriends ++;
+  }
+  assert (Nfriends <= object[idx1].Nindex);
+  assert (Nfriends <= Nobject);
+  fprintf (stderr, "object %d sequenced by object %d, %d friends sequenced\n", idx1, idx2, Nfriends);
+  return;
+}
+
+int spexseq (int argc, char **argv) {
+  
+  int i, j, k;
+  Vector *index1, *index2, *distance;
+  Object *object = NULL;
+  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;
+
+  // XXX enforce matching lengths on the three vectors
+  if ((index1   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((index2   = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((distance = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
+
+  float f1 = atof (argv[4]);
+  float f2 = atof (argv[5]);
+
+  // int idx1 = atoi(argv[9]);
+  // int idx2 = atoi(argv[10]);
+
+  CastVector (index1, OPIHI_INT);
+  CastVector (index2, OPIHI_INT);
+
+  // how many objects do we have?
+  Nobject = 0;
+  for (i = 0; i < index1->Nelements; i++) {
+    Nobject = MAX (Nobject, index1->elements.Int[i]);
+    Nobject = MAX (Nobject, index2->elements.Int[i]);
+  }
+  Nobject ++;  // after the loop, Nobject has the value of the highest index, not the count
+
+  // allocate the list of Object -- these list the possible friends of the given object[i] 
+  ALLOCATE (object, Object, Nobject);
+  for (i = 0; i < Nobject; i++) {
+    ALLOCATE (object[i].index, int, Nobject);
+    memset (object[i].index, 0, Nobject*sizeof(int));
+    ALLOCATE (object[i].Dtgt, float, Nobject);
+    memset (object[i].Dtgt, 0, Nobject*sizeof(float));
+    object[i].NINDEX = Nobject;
+    object[i].Nindex = 0;
+  }
+
+  // generate the set of vectors of all distances for each entry
+  // the object seq number (object[i]) will be used to match the index1 and index2 values
+  for (i = 0; i < index1->Nelements; i++) {
+    int N, m;
+
+    N = index1->elements.Int[i];
+    if (N >= Nobject) abort();
+    m = object[N].Nindex;
+    object[N].index[m] = index2->elements.Int[i];
+    object[N].Dtgt[m] = distance->elements.Flt[i];
+    object[N].Nindex ++;
+    if (object[N].Nindex == object[N].NINDEX) {
+      object[N].NINDEX += 100;
+      REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
+      REALLOCATE (object[N].index, int, object[N].NINDEX);
+    }
+
+    // Above we used the set (N,M,D_N,M) to assigned the distance of obj M to obj N.
+    // Below we use that set to assign the distance of obj N to obj M.  IF we supply all
+    // pairwise distances, then this duplicates all entries, so we need to skip that step
+    // in such a case.
+
+    // XX N = index2->elements.Int[i];
+    // XX if (N >= Nobject) abort();
+    // XX m = object[N].Nindex;
+    // XX object[N].index[m] = index1->elements.Int[i];
+    // XX object[N].Dtgt[m] = distance->elements.Flt[i];
+    // XX object[N].Nindex ++;
+    // XX if (object[N].Nindex == object[N].NINDEX) {
+    // XX   object[N].NINDEX += 100;
+    // XX   REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
+    // XX   REALLOCATE (object[N].index, int, object[N].NINDEX);
+    // XX }
+  }
+
+  for (i = 0; i < Nobject; i++) {
+    // sort so closest friends are first
+    sortfriends (object[i].Dtgt, object[i].index, object[i].Nindex);
+    object[i].shifted = FALSE;
+    object[i].So = NAN;
+    object[i].dSo = NAN;
+    object[i].Xo = NAN;
+    object[i].dXo = NAN;
+    // generate the reverse index
+    REALLOCATE (object[i].rindex, int, object[i].Nindex);
+    REALLOCATE (object[i].seq, float, object[i].Nindex);
+    for (j = 0; j < object[i].Nindex; j++) {
+      int N = object[i].index[j];
+      object[i].rindex[N] = j;
+      object[i].seq[j] = NAN;
+    }
+  }
+
+  // find and save the max distance
+  float Dmax = 0;
+  for (i = 0; i < distance->Nelements; i++) {
+    Dmax = MAX (Dmax, distance->elements.Flt[i]);
+  }
+
+  // find two objects at the opposite extremes
+  float Dm1 = object[0].Dtgt[object[0].Nindex - 1];
+  for (i = 0; (object[0].Dtgt[i] < 0.75*Dm1) && (i < object[0].Nindex); i++);
+  int pin1 = object[0].index[i];
+  assert (pin1 >= 0);
+  assert (pin1 < Nobject);
+
+  float Dm2 = object[pin1].Dtgt[object[pin1].Nindex - 1];
+  for (i = 0; (object[pin1].Dtgt[i] < 0.75*Dm2) && (i < object[pin1].Nindex); i++);
+  int pin2 = object[pin1].index[i];
+  assert (pin2 >= 0);
+  assert (pin2 < Nobject);
+
+  // generate sequences for the near friends of all objects based on a somewhat distant friend
+  for (i = 0; i < Nobject; i++) {
+    get_sequence (object, Nobject, i, Dmax, f1, f2, pin1, pin2);
+  }
+
+  // now we need to reconcile these sequences.
+  // start with object 0, find shift of all friends relative to 0
+  object[0].shifted = TRUE;
+  for (i = 0; i < Nobject; i++) {
+
+    // don't use unreconciled objects as a reference (we'll catch them eventually)
+    if (!object[i].shifted) continue;
+
+    int Nfixed = 0;
+
+    // shift all of the close friends of object (i) wrt object (i)
+    for (j = 0; (j < object[i].Nindex) && (object[i].Dtgt[j] < f2*Dmax); j++) {
+
+      // loop over the sequence for object Nj
+      int Nj = object[i].index[j];
+
+      // only reconcile once
+      if (object[Nj].shifted == TRUE) continue;
+
+      float dS1p = 0.0;
+      float dS2p = 0.0;
+      float dS1m = 0.0;
+      float dS2m = 0.0;
+      int nS = 0;
+
+      float Smin = +1000.0;
+      float Smax = -1000.0;
+
+      for (k = 0; k < object[Nj].Nindex; k++) {
+	int Nsj = object[Nj].index[k];
+	int Nsi = object[i].rindex[Nsj];
+	
+	if (isnan(object[Nj].seq[k])) continue;
+	if (isnan(object[i].seq[Nsi])) continue;
+	
+	Smin = MIN(object[Nj].seq[k], Smin);
+	Smax = MAX(object[Nj].seq[k], Smax);
+
+	float dSp = object[Nj].seq[k] - object[i].seq[Nsi];
+	float dSm = -object[Nj].seq[k] - object[i].seq[Nsi];
+	dS1p += dSp;
+	dS2p += dSp*dSp;
+	dS1m += dSm;
+	dS2m += dSm*dSm;
+	nS ++;
+      }
+      float Sp = dS1p / nS;
+      float Sm = dS1m / nS;
+      float dSp = sqrt(dS2p / nS - Sp*Sp);
+      float dSm = sqrt(dS2m / nS - Sm*Sm);
+      fprintf (stderr, "%d: %f +/- %f or %f +/- %f : %f range, %d matches\n", Nj, Sp, dSp, Sm, dSm, Smax - Smin, nS);
+
+      if (dSm > dSp) {
+	// parity is always 1 if we use the pins (but it does not seem to work)
+	for (k = 0; k < object[Nj].Nindex; k++) {
+	  if (isnan(object[Nj].seq[k])) continue;
+	  object[Nj].seq[k] -= Sp;
+	}
+	object[Nj].So = Sp;
+	object[Nj].dSo = dSp;
+      } else {
+	// flip the sequence and offset
+	for (k = 0; k < object[Nj].Nindex; k++) {
+	  if (isnan(object[Nj].seq[k])) continue;
+	  object[Nj].seq[k] = -object[Nj].seq[k] - Sm; 
+	}
+	object[Nj].So = Sm;
+	object[Nj].dSo = dSm;
+      }
+      object[Nj].Srange = Smax - Smin;
+      object[Nj].nSeq = nS;
+      object[Nj].shifted = TRUE;
+      Nfixed ++;
+    }
+    // fprintf (stderr, "%d used to reconcile %d friends\n", i, Nfixed);
+  }
+
+  int Nfixed = 0;
+  int Nunfix = 0;
+  for (i = 0; i < Nobject; i++) {
+    if (object[i].shifted) { 
+      Nfixed ++;
+    } else {
+      Nunfix ++;
+    }
+  }
+  fprintf (stderr, "%d fixed, %d unfixed\n", Nfixed, Nunfix);
+      
+
+  // calculate the mean Xo for each object using the values of the sequences calculated
+  for (i = 0; i < Nobject; i++) {
+
+    if (!object[i].shifted) {
+      fprintf (stderr, "this object was not reconciled\n");
+    }
+    float XoSum = 0.0;
+    float XoS2 = 0.0;
+    int XoNum = 0;
+
+    // find all the measurements of this object's sequence (from its friends values)
+    for (j = 0; j < Nobject; j++) {
+      int Ns = object[j].rindex[i];
+      if (isnan(object[j].seq[Ns])) continue;
+      XoSum += object[j].seq[Ns];
+      XoS2  += SQ(object[j].seq[Ns]);
+      XoNum ++;
+    }
+    float Xo = XoSum / XoNum;
+    float dXo = sqrt(XoS2 / XoNum - Xo*Xo);
+    object[i].Xo = Xo;
+    object[i].dXo = dXo;
+  }
+
+  // XXX seq is effectively Xo, calculate Yo for the objects based on the local 
+  // objects
+
+  // save the result
+  {
+    Vector *outindex = SelectVector ("sp1d_idx", ANYVECTOR, TRUE); if (!outindex) goto escape;
+    Vector *outXo    = SelectVector ("sp1d_Xo",  ANYVECTOR, TRUE); if (!outXo) goto escape;
+    Vector *outYo    = SelectVector ("sp1d_Yo",  ANYVECTOR, TRUE); if (!outYo) goto escape;
+
+    Vector *outSo    = SelectVector ("sp1d_So",  ANYVECTOR, TRUE); if (!outXo) goto escape;
+    Vector *outdSo   = SelectVector ("sp1d_dSo", ANYVECTOR, TRUE); if (!outYo) goto escape;
+    Vector *outdXo   = SelectVector ("sp1d_dXo", ANYVECTOR, TRUE); if (!outXo) goto escape;
+
+    Vector *outSr   = SelectVector ("sp1d_Sr", ANYVECTOR, TRUE); if (!outSr) goto escape;
+    Vector *outSn   = SelectVector ("sp1d_nS", ANYVECTOR, TRUE); if (!outSn) goto escape;
+
+    ResetVector (outindex, OPIHI_INT, Nobject);
+    ResetVector (outXo, OPIHI_FLT, Nobject);
+    ResetVector (outYo, OPIHI_FLT, Nobject);
+
+    ResetVector (outSo, OPIHI_FLT, Nobject);
+    ResetVector (outdSo, OPIHI_FLT, Nobject);
+    ResetVector (outdXo, OPIHI_FLT, Nobject);
+
+    ResetVector (outSr, OPIHI_FLT, Nobject);
+    ResetVector (outSn, OPIHI_INT, Nobject);
+
+    for (i = 0; i < Nobject; i++) {
+      outindex->elements.Int[i] = i;
+      outXo->elements.Flt[i] = object[i].Xo;
+      outYo->elements.Flt[i] = object[i].Yo;
+      outSo->elements.Flt[i] = object[i].So;
+      outdSo->elements.Flt[i] = object[i].dSo;
+      outdXo->elements.Flt[i] = object[i].dXo;
+      outSr->elements.Flt[i] = object[i].Srange;
+      outSn->elements.Int[i] = object[i].nSeq;
+    }
+  }
+
+  return TRUE;
+  
+escape: 
+  gprint (GP_ERR, "invalid vector\n");
+  return FALSE;
+  
+usage:
+  gprint (GP_ERR, "USAGE: spexseq (index1) (index2) (distance) (Niter) (nCloseMax) (nCloseIter) (farFrac) (maxPressure) (idx1) (idx2)\n");
+  return FALSE;
+}
+
+/* this function takes 3 vectors: a set of distances and a pair of indicies identifying the end points.  
+   
+   the goal is to generate groups of the end points that are relatively closer.
+
+   note that the distances may not be Euclidean: nothing can be assumed about relationships between d(1,2), d(2,3), and d(1,3)
+
+   the indicies should be sequential.  the algorithm does not require it, but storage is much cleaner if it is
+
+*/