Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/Makefile	(revision 24991)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/Makefile	(revision 24992)
@@ -36,4 +36,7 @@
 $(SRC)/medianmap.$(ARCH).o	   \
 $(SRC)/mkgauss.$(ARCH).o	   \
+$(SRC)/mksersic.$(ARCH).o	   \
+$(SRC)/galsectors.$(ARCH).o	   \
+$(SRC)/galprofiles.$(ARCH).o	   \
 $(SRC)/multifit.$(ARCH).o	   \
 $(SRC)/objload.$(ARCH).o	   \
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galcontour.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galcontour.c	(revision 24992)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galcontour.c	(revision 24992)
@@ -0,0 +1,89 @@
+# include "astro.h"
+
+XXX probably not a great algorithm...
+
+
+int galcontour (int argc, char **argv) {
+  
+  int i, j, Nx, Ny, N, Nsec, NELEMENTS, Npt;
+  float *in;
+  double theta, dtheta, x, y, Xo, Yo;
+  Buffer *buf;
+  Vector *xvec, *yvec, *tvec;
+  char name[128];
+
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: galcontour (buffer) (Xo) (Yo) (flux)\n");
+    gprint (GP_ERR, "  generate contour points (Cx) (Cy) (Ct)\n");
+    return (FALSE);
+  }
+
+  /* select input / output buffers */
+  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  Nx = buf[0].header.Naxis[0];
+  Ny = buf[0].header.Naxis[1];
+  
+  Xo = atof(argv[2]);
+  Yo = atof(argv[3]);
+  flux = atof(argv[4]);
+
+  if ((xvec = SelectVector ("Cx", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((yvec = SelectVector ("Cy", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((tvec = SelectVector ("Ct", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  N = 0;
+  NELEMENTS = 100;
+  ResetVector (xvec, OPIHI_FLT, NELEMENTS);
+  ResetVector (yvec, OPIHI_FLT, NELEMENTS);
+  ResetVector (tvec, OPIHI_FLT, NELEMENTS);
+
+  in = (float *) buf[0].matrix.buffer;
+
+  // find the flux transition for each row until not found
+
+  // first go up:
+  for (j = Yo; j < Ny; j++) {
+    
+    // first go left
+    for (i = Xo; i >= 0; i--) {
+      F = in[i + j*Nx];
+      if (F < flux) {
+	// interpolate to flux value between f[i,j] and f[i+1,j]
+      }
+    }
+
+  for (j = 0; j < Ny; j++) {
+    for (i = 0; i < Nx; i++, in++) {
+      x = i - Xo;
+      y = j - Yo;
+      theta = DEG_RAD*atan2(y, x);
+      if (theta < -0.5*dtheta) theta += 360.0; // need to allow -0.5dtheta for first bin
+      if (theta > 360 - 0.5*dtheta) theta -= 360.0; // need to allow -0.5dtheta for first bin
+      N = (int)(theta / dtheta + 0.5);
+      if (N < 0) {
+	fprintf (stderr, "?");
+	continue;
+      }
+      if (N >= Nsec) {
+	fprintf (stderr, "!");
+	continue;
+      }
+      
+      fvec[N][0].elements.Flt[fvec[N][0].Nelements] = *in;
+      rvec[N][0].elements.Flt[rvec[N][0].Nelements] = hypot(x, y);
+      fvec[N][0].Nelements ++;
+      rvec[N][0].Nelements ++;
+
+      *in = N;
+
+      if (fvec[N][0].Nelements >= NELEMENTS) {
+	NELEMENTS += 1000;
+	for (N = 0; N < Nsec; N++) {
+	  REALLOCATE (fvec[N][0].elements.Flt, opihi_flt, NELEMENTS);
+	  REALLOCATE (rvec[N][0].elements.Flt, opihi_flt, NELEMENTS);
+	}
+      }
+    }
+  }
+  return (TRUE);
+}
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galprofiles.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galprofiles.c	(revision 24992)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galprofiles.c	(revision 24992)
@@ -0,0 +1,187 @@
+# include "astro.h"
+
+double Interpolate (const double x, const double y, float *in, int Nx, int Ny);
+
+int galprofiles (int argc, char **argv) {
+  
+  int i, j, Nx, Ny, N, Nsec, NELEMENTS;
+  float *in;
+  double theta, dtheta, x, y, Xo, Yo, r, Rmax, value;
+  Buffer *buf;
+  Vector **fvec, **rvec, *tvec;
+  char name[128];
+
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: galprofiles (buffer) (Xo) (Yo) Nsec\n");
+    gprint (GP_ERR, "  generates Nsec profiles around the circle; the first is centered on 0.0 degrees\n");
+    return (FALSE);
+  }
+
+  /* select input / output buffers */
+  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  Nx = buf[0].header.Naxis[0];
+  Ny = buf[0].header.Naxis[1];
+  Rmax = 0.5*hypot(Nx, Ny);
+  
+  Xo = atof(argv[2]);
+  Yo = atof(argv[3]);
+  Nsec = atof(argv[4]);
+
+  dtheta = 360.0 / Nsec;
+  // sector i ranges from (i-0.5)*dtheta to (i+0.5)*dtheta
+  // i = theta / dtheta + 0.5;
+
+  ALLOCATE (fvec, Vector *, Nsec);
+  ALLOCATE (rvec, Vector *, Nsec);
+  if ((tvec = SelectVector ("theta", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  ResetVector (tvec, OPIHI_FLT, Nsec);
+
+  // dereference pointer
+  in = (float *) buf[0].matrix.buffer;
+
+  for (i = 0; i < Nsec; i++) {
+
+    // generate the output vector names & vectors
+    sprintf (name, "flux_%d", i);
+    if ((fvec[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) return (FALSE);
+    sprintf (name, "rad_%d", i);
+    if ((rvec[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+    // allocate initial data space
+    N = 0;
+    NELEMENTS = 1000;
+    ResetVector (fvec[i], OPIHI_FLT, NELEMENTS);
+    ResetVector (rvec[i], OPIHI_FLT, NELEMENTS);
+
+    // angle for this profile
+    tvec[0].elements.Flt[i] = i*dtheta;
+    theta = RAD_DEG*i*dtheta;
+
+    // start at Xo,Yo and find the x,y locations for r_i, theta where r_i increments by 1 pixel
+    for (r = 0; r < Rmax; r++) {
+
+      // XXX Xo,Yo are referenced to pixels with bounds i+0.0, i+1.0
+      x = r * cos (theta) + Xo;
+      y = r * sin (theta) + Yo;
+
+      // value is NAN if we run off the image
+      value = Interpolate(x, y, in, Nx, Ny);
+      if (isnan(value)) break;
+
+      rvec[i][0].elements.Flt[N] = r;
+      fvec[i][0].elements.Flt[N] = value;
+      N++;
+
+      if (N >= NELEMENTS) {
+	NELEMENTS += 100;
+	REALLOCATE (rvec[N][0].elements.Flt, opihi_flt, NELEMENTS);
+	REALLOCATE (fvec[N][0].elements.Flt, opihi_flt, NELEMENTS);
+      }
+    }
+    fvec[i][0].Nelements = N;
+    rvec[i][0].Nelements = N;
+  }
+  return (TRUE);
+}
+
+// fast & simple API to interpolate to a subpixel position using bilinear interpolation
+// x,y in parent image coordinates (pixel centers at 0.5, 0.5)
+// stolen from psLib
+double Interpolate (const double x, const double y, float *in, int Nx, int Ny) {
+
+  // allow extrapolation a small distance beyond the edge of valid pixels, but no
+  // further (this allows the nXskip,nYskip boundary areas to be used as well)
+  float nXedge = 0.125*Nx;
+  float nYedge = 0.125*Ny;
+
+  if ((x < -nXedge) || (x > Nx + nXedge) || (y < -nYedge) || (y > Ny + nYedge)) {
+    return NAN;
+  }
+
+  // limiting cases: Nx == 1 and/or Ny == 1
+
+  // if we have a single pixel, there is no spatial information
+  if ((Nx == 1) && (Ny == 1)) {
+    return in[0];
+  }
+
+  // handle edge cases with extrapolation
+  const int ix = x - 0.5; // index of reference pixel
+  const int iy = y - 0.5; // index of reference pixel
+
+  // do numCols,Rows first so we are never < 0
+  const int Xs = MAX (MIN (ix, Nx - 2), 0);
+  const int Ys = MAX (MIN (iy, Ny - 2), 0);
+
+  const int Xe = Xs + 1;
+  const int Ye = Ys + 1;
+
+  // dx,dy range from 0.0 to 1.0 for interpolated pixels, and -0.5 to 1.5 for extrapolation
+  const double dx = x - 0.5 - Xs;
+  const double dy = y - 0.5 - Ys;
+
+  const double rx = 1.0 - dx;
+  const double ry = 1.0 - dy;
+
+  // if Nx == 1, we have no x-dir spatial information
+  if (Nx == 1) {
+    double V0 = in[Ys*Nx + Xs];
+    double V1 = in[Ye*Nx + Xs];
+
+    const double value = V0*ry + V1*dy;
+    return value;
+  }
+
+  // if Ny == 1, we have no y-dir spatial information
+  if (Ny == 1) {
+    double V0 = in[Ys*Nx + Xs];
+    double V1 = in[Ys*Nx + Xe];
+
+    const double value = V0*rx + V1*dx;
+    return value;
+  }
+
+  // Vxy
+  double V00 = in[Ys*Nx + Xs];
+  double V01 = in[Ye*Nx + Xs];
+  double V10 = in[Ys*Nx + Xe];
+  double V11 = in[Ye*Nx + Xe];
+
+  double value;
+
+  // corners
+  if ((dx < 0.0) && (dy < 0.0)) {
+    return V00;
+  }
+  if ((dx > 1.0) && (dy < 0.0)) {
+    return V10;
+  }
+  if ((dx < 0.0) && (dy > 1.0)) {
+    return V01;
+  }
+  if ((dx > 1.0) && (dy > 1.0)) {
+    return V11;
+  }
+
+  // sides
+  if (dx < 0.0) {
+    value = V00*ry + V01*dy;
+    return value;
+  }
+  if (dy < 0.0) {
+    value = V00*rx + V10*dx;
+    return value;
+  }
+  if (dx > 1.0) {
+    value = V10*ry + V11*dy;
+    return value;
+  }
+  if (dy > 1.0) {
+    value = V01*rx + V11*dx;
+    return value;
+  }
+
+  // bilinear interpolation
+  value = V00*rx*ry + V10*dx*ry + V01*rx*dy + V11*dx*dy;
+  return value;
+}
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galsectors.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galsectors.c	(revision 24992)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galsectors.c	(revision 24992)
@@ -0,0 +1,107 @@
+# include "astro.h"
+
+int galsectors (int argc, char **argv) {
+  
+  int i, j, Nx, Ny, N, Nsec, NELEMENTS, Npt;
+  float *in;
+  double theta, dtheta, x, y, Xo, Yo;
+  Buffer *buf;
+  Vector **fvec, **rvec, *tvec, *Rv, *Tv, *Fv;
+  char name[128];
+
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: galsectors (buffer) (Xo) (Yo) Nsec\n");
+    gprint (GP_ERR, "  generates Nsec sectors around the circle; the first is centered on 0.0 degrees\n");
+    return (FALSE);
+  }
+
+  /* select input / output buffers */
+  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  Nx = buf[0].header.Naxis[0];
+  Ny = buf[0].header.Naxis[1];
+  
+  Xo = atof(argv[2]);
+  Yo = atof(argv[3]);
+  Nsec = atof(argv[4]);
+
+  // for a test, generate R, theta vectors for the whole image:
+  if ((Rv = SelectVector ("Rvec", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Tv = SelectVector ("Tvec", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Fv = SelectVector ("Fvec", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  ResetVector (Rv, OPIHI_FLT, Nx*Ny);
+  ResetVector (Tv, OPIHI_FLT, Nx*Ny);
+  ResetVector (Fv, OPIHI_FLT, Nx*Ny);
+
+  in = (float *) buf[0].matrix.buffer;
+  Npt = 0;
+  for (j = 0; j < Ny; j++) {
+    for (i = 0; i < Nx; i++, in++) {
+      x = i - Xo;
+      y = j - Yo;
+      Rv[0].elements.Flt[Npt] = hypot(x, y);
+      Tv[0].elements.Flt[Npt] = DEG_RAD*atan2(y, x);
+      Fv[0].elements.Flt[Npt] = *in;
+      Npt ++;
+    }
+  }
+  return (TRUE);
+
+  dtheta = 360.0 / Nsec;
+  // sector i ranges from (i-0.5)*dtheta to (i+0.5)*dtheta
+  // i = theta / dtheta + 0.5;
+
+  ALLOCATE (fvec, Vector *, Nsec);
+  ALLOCATE (rvec, Vector *, Nsec);
+  if ((tvec = SelectVector ("theta", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  ResetVector (tvec, OPIHI_FLT, Nsec);
+
+  // generate the output vector names
+  NELEMENTS = 1000;
+  for (i = 0; i < Nsec; i++) {
+    sprintf (name, "flux_%d", i);
+    if ((fvec[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) return (FALSE);
+    sprintf (name, "rad_%d", i);
+    if ((rvec[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) return (FALSE);
+    ResetVector (fvec[i], OPIHI_FLT, NELEMENTS);
+    ResetVector (rvec[i], OPIHI_FLT, NELEMENTS);
+    fvec[i][0].Nelements = 0;
+    rvec[i][0].Nelements = 0;
+    tvec[0].elements.Flt[i] = i*dtheta;
+  }
+
+  in = (float *) buf[0].matrix.buffer;
+  for (j = 0; j < Ny; j++) {
+    for (i = 0; i < Nx; i++, in++) {
+      x = i - Xo;
+      y = j - Yo;
+      theta = DEG_RAD*atan2(y, x);
+      if (theta < -0.5*dtheta) theta += 360.0; // need to allow -0.5dtheta for first bin
+      if (theta > 360 - 0.5*dtheta) theta -= 360.0; // need to allow -0.5dtheta for first bin
+      N = (int)(theta / dtheta + 0.5);
+      if (N < 0) {
+	fprintf (stderr, "?");
+	continue;
+      }
+      if (N >= Nsec) {
+	fprintf (stderr, "!");
+	continue;
+      }
+      
+      fvec[N][0].elements.Flt[fvec[N][0].Nelements] = *in;
+      rvec[N][0].elements.Flt[rvec[N][0].Nelements] = hypot(x, y);
+      fvec[N][0].Nelements ++;
+      rvec[N][0].Nelements ++;
+
+      *in = N;
+
+      if (fvec[N][0].Nelements >= NELEMENTS) {
+	NELEMENTS += 1000;
+	for (N = 0; N < Nsec; N++) {
+	  REALLOCATE (fvec[N][0].elements.Flt, opihi_flt, NELEMENTS);
+	  REALLOCATE (rvec[N][0].elements.Flt, opihi_flt, NELEMENTS);
+	}
+      }
+    }
+  }
+  return (TRUE);
+}
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/init.c	(revision 24991)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/init.c	(revision 24992)
@@ -22,5 +22,8 @@
 int imsub                   PROTO((int, char **));
 int medianmap               PROTO((int, char **));
+int galsectors              PROTO((int, char **));
+int galprofiles             PROTO((int, char **));
 int mkgauss                 PROTO((int, char **));
+int mksersic                PROTO((int, char **));
 int multifit                PROTO((int, char **));
 int objload                 PROTO((int, char **));
@@ -61,4 +64,7 @@
   {1, "medianmap",   medianmap,    "small median image"},
   {1, "mkgauss",     mkgauss,      "generate a 2-D gaussian centered in image"},
+  {1, "mksersic",    mksersic,     "generate a 2-D sersic profile"},
+  {1, "galsectors",  galsectors,   "generate radial vectors for sectors of width dtheta"},
+  {1, "galprofiles", galprofiles,  "generate radial vectors with interpolation along paths"},
   {1, "multifit",    multifit,     "fit multi-order spectrum"},
   {1, "objload",     objload,      "plot obj data on Ximage "},
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/mksersic.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/mksersic.c	(revision 24992)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/mksersic.c	(revision 24992)
@@ -0,0 +1,87 @@
+# include "astro.h"
+
+int mksersic (int argc, char **argv) {
+  
+  int i, j, Nx, Ny, N, center;
+  float *in;
+  double Rmaj, Rmin, Theta, alpha, Io, ARatio;
+  double root1, root2, R, A1, A2, A3;
+  double Sx, Sy, Sxy;
+  double x, y, r, f, Xo, Yo;
+  Buffer *buf;
+
+  Xo = Yo = 0;
+  center = TRUE;
+  if ((N = get_argument (argc, argv, "-coord"))) {
+    remove_argument (N, &argc, argv);
+    Xo = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    Yo = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    center = FALSE;
+  }
+
+  Theta = 0.0;
+  if ((N = get_argument (argc, argv, "-angle"))) {
+    remove_argument (N, &argc, argv);
+    Theta = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  ARatio = 1.0;
+  if ((N = get_argument (argc, argv, "-aratio"))) {
+    remove_argument (N, &argc, argv);
+    ARatio = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: mksersic (buffer) (Io) (Rmaj) (alpha)\n");
+    return (FALSE);
+  }
+
+  /* select input / output buffers */
+  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  Nx = buf[0].header.Naxis[0];
+  Ny = buf[0].header.Naxis[1];
+  if (center) {
+    Xo = 0.5*Nx;
+    Yo = 0.5*Ny;
+  }
+  
+  Io = atof (argv[2]);
+
+  /* shape parameters */
+  Rmaj = atof (argv[3]);
+  Rmin = Rmaj / ARatio;
+
+  alpha = atof (argv[4]);
+
+  /* given Rmaj, Rmin, Theta, find Sx, Sy, Sxy */
+  root1 = SQ(1.0 / Rmaj);
+  root2 = SQ(1.0 / Rmin);
+
+  // XXX check this
+  R = 0.5 * (root1 - root2);
+  A1 = 0.25*(root1 + root2) - 0.5*R*cos(2*RAD_DEG*Theta);
+  A2 = 0.25*(root1 + root2) + 0.5*R*cos(2*RAD_DEG*Theta);
+  A3 = -R*sin(2*RAD_DEG*Theta);
+
+  Sx = 0.5/A1;
+  Sy = 0.5/A2;
+  Sxy = A3;
+
+  /* f = exp (-r^alpha), r = (x^2 / 2Sx) + (y^2 / 2Sy) + Sxy*x*y */
+  in = (float *) buf[0].matrix.buffer;
+  for (j = 0; j < Ny; j++) {
+    for (i = 0; i < Nx; i++, in++) {
+      x = i + 0.5 - Xo;
+      y = j + 0.5 - Yo;
+      r = pow((0.5*x*x/Sx + 0.5*y*y/Sy + x*y*Sxy), alpha);
+      f = Io*exp (-r);
+      *in += f;
+    }
+  }
+
+  return (TRUE);
+}
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/Makefile
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/Makefile	(revision 24991)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/Makefile	(revision 24992)
@@ -110,5 +110,5 @@
 $(SRC)/style.$(ARCH).o		   \
 $(SRC)/subraster.$(ARCH).o	   \
-$(SRC)/subset.$(ARCH).o	   \
+$(SRC)/subset.$(ARCH).o	           \
 $(SRC)/svd.$(ARCH).o		   \
 $(SRC)/swapbytes.$(ARCH).o	   \
@@ -118,11 +118,13 @@
 $(SRC)/tvcolors.$(ARCH).o	   \
 $(SRC)/tvcontour.$(ARCH).o	   \
-$(SRC)/tvgrid.$(ARCH).o	   \
+$(SRC)/tvgrid.$(ARCH).o	           \
 $(SRC)/uniq.$(ARCH).o		   \
-$(SRC)/unsign.$(ARCH).o	   \
+$(SRC)/unsign.$(ARCH).o	           \
 $(SRC)/vbin.$(ARCH).o		   \
+$(SRC)/vgroup.$(ARCH).o		   \
 $(SRC)/vclip.$(ARCH).o		   \
-$(SRC)/vgauss.$(ARCH).o           \
-$(SRC)/vmaxwell.$(ARCH).o           \
+$(SRC)/vgauss.$(ARCH).o            \
+$(SRC)/vellipse.$(ARCH).o          \
+$(SRC)/vmaxwell.$(ARCH).o          \
 $(SRC)/vgrid.$(ARCH).o		   \
 $(SRC)/vload.$(ARCH).o		   \
@@ -130,5 +132,5 @@
 $(SRC)/vpop.$(ARCH).o		   \
 $(SRC)/vroll.$(ARCH).o		   \
-$(SRC)/vsmooth.$(ARCH).o	\
+$(SRC)/vsmooth.$(ARCH).o	   \
 $(SRC)/vstats.$(ARCH).o		   \
 $(SRC)/wd.$(ARCH).o		   \
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/init.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/init.c	(revision 24991)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/init.c	(revision 24992)
@@ -110,8 +110,10 @@
 int unsign           PROTO((int, char **));
 int vbin             PROTO((int, char **));
+int vgroup           PROTO((int, char **));
 int vclip            PROTO((int, char **));
 int vect_select      PROTO((int, char **));
 int vgrid            PROTO((int, char **));
 int vgauss           PROTO((int, char **));
+int vellipse         PROTO((int, char **));
 int vmaxwell         PROTO((int, char **));
 int vload            PROTO((int, char **));
@@ -243,8 +245,10 @@
   {1, "uniq",         uniq,             "create a uniq vector subset from a vector"},
   {1, "unsign",       unsign,           "toggle the UNSIGN status"},
-  {1, "vbin",         vbin,             "rebin vector data by a fector of N"},
+  {1, "vbin",         vbin,             "rebin vector data by a factor of N"},
+  {1, "vgroup",       vgroup,           "group y vector into bins defined by x vector values"},
   {1, "vclip",        vclip,            "clip values in a vector to be within a range"},
   {1, "vectors",      list_vectors,     "list vectors"},
   {1, "vgauss",       vgauss,           "fit a Gaussian to a vector"},
+  {1, "vellipse",     vellipse,         "fit a Ellipse to a vector pair"},
   {1, "vgrid",        vgrid,            "generate an image from a triplet of vectors"},
   {1, "vhistogram",   histogram,        "generate histogram from vector"},
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/vellipse.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/vellipse.c	(revision 24992)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/vellipse.c	(revision 24992)
@@ -0,0 +1,272 @@
+# include "data.h"
+
+/* local private functions */
+opihi_flt fellipseOD (opihi_flt theta, opihi_flt *par, int Npar, opihi_flt *dpar);
+
+// enum {PAR_X0, PAR_Y0, PAR_RX, PAR_RY, PAR_T0, PAR_P0};
+enum {PAR_RMIN, PAR_EPSILON, PAR_PHI};
+
+# define GET_VAR(V,A) { \
+  char *c; \
+  c = get_variable (A); \
+  if (c == NULL) { \
+    gprint (GP_ERR, "missing fit parameter A\n"); \
+    return (FALSE); \
+  } \
+  V = atof (c); \
+  free (c); }
+
+int vellipse (int argc, char **argv) {
+
+  int i, N, Npts, Npar, Quiet;
+  opihi_flt par[6], *pos, *dpos, *theta, **covar;
+  opihi_flt chisq, ochisq, dchisq;
+  // opihi_flt Rmaj;
+  opihi_flt Xmin, Xmax, Ymin, Ymax;
+  Vector *xvec, *yvec, *Xvec, *Yvec, *tvec;
+
+  Quiet = FALSE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    Quiet = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-quiet"))) {
+    Quiet = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 6) {
+    gprint (GP_ERR, "USAGE: vellipse <xobs> <yobs> <xfit> <yfit> (theta)\n");
+    // gprint (GP_ERR, " uses guesses: E_X0, E_Y0, E_RMAJ, E_RMIN, E_THETA\n");
+    gprint (GP_ERR, " uses guesses: E_RMAJ, E_RMIN, E_PHI\n");
+    return (FALSE);
+  }
+  
+  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Xvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Yvec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((tvec = SelectVector (argv[5], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+
+  REQUIRE_VECTOR_FLT (xvec, FALSE); 
+  REQUIRE_VECTOR_FLT (yvec, FALSE); 
+
+  // GET_VAR (par[PAR_X0], "E_X0");
+  // GET_VAR (par[PAR_Y0], "E_Y0");
+  // GET_VAR (par[PAR_RX], "E_RMAJ");
+  // GET_VAR (par[PAR_RY], "E_RMIN");
+  // GET_VAR (par[PAR_T0], "E_THETA");
+  // GET_VAR (par[PAR_P0], "E_PHI");
+  // Npar = 6;
+
+  // we need to generate a vector of alternating x,y values and the independent variable theta:
+  Xmax = Ymax = Xmin = Ymin = 0.0;
+  ALLOCATE (pos,   opihi_flt, 2*xvec[0].Nelements); 
+  ALLOCATE (dpos,  opihi_flt, 2*xvec[0].Nelements); 
+  ALLOCATE (theta, opihi_flt, 2*xvec[0].Nelements); 
+  for (i = 0; i < xvec[0].Nelements; i++) {
+    pos[2*i+0]   = xvec[0].elements.Flt[i];
+    pos[2*i+1]   = yvec[0].elements.Flt[i];
+    Xmin = MIN (Xmin, xvec[0].elements.Flt[i]);
+    Xmax = MAX (Xmax, xvec[0].elements.Flt[i]);
+    Ymin = MIN (Ymin, yvec[0].elements.Flt[i]);
+    Ymax = MAX (Ymax, yvec[0].elements.Flt[i]);
+    dpos[2*i+0]  = 1.0;
+    dpos[2*i+1]  = 1.0;
+    theta[2*i+0] = atan2(yvec[0].elements.Flt[i], xvec[0].elements.Flt[i]);
+    theta[2*i+1] = theta[2*i+0];
+    // fprintf (stderr, "x,y: %f, %f -> %f, %f\n", yvec[0].elements.Flt[i], xvec[0].elements.Flt[i], yvec[0].elements.Flt[i] - par[PAR_Y0], xvec[0].elements.Flt[i] - par[PAR_X0]);
+    // fprintf (stderr, "x,y: %f, %f -> %f, %f\n", yvec[0].elements.Flt[i], xvec[0].elements.Flt[i], yvec[0].elements.Flt[i] - par[PAR_Y0], xvec[0].elements.Flt[i] - par[PAR_X0]);
+    // fprintf (stderr, "theta: %f - %f : %f\n", DEG_RAD*theta[2*i+0], tvec[0].elements.Flt[i], DEG_RAD*theta[2*i+0] - tvec[0].elements.Flt[i]);
+  }
+  Npts = 2*xvec[0].Nelements;
+
+  // basic guess from range, PHI in degrees, converted to radians
+  // GET_VAR (Rmaj,          "E_RMAJ");
+  // GET_VAR (par[PAR_RMIN], "E_RMIN");
+  // GET_VAR (par[PAR_PHI],  "E_PHI");
+  par[PAR_PHI] = 0.0;
+  par[PAR_EPSILON] = 0.5;
+  par[PAR_RMIN] = 0.25*((Xmax - Xmin) + (Ymax - Ymin));
+  Npar = 3;
+
+  ochisq = mrqinit (theta, pos, dpos, Npts, par, Npar, fellipseOD, !Quiet);
+  dchisq = ochisq + 2*Npts;
+  chisq  = ochisq + dchisq;
+
+  for (i = 0; (i < 20) && (chisq > 1e-3) && ((dchisq > 0.001*(Npts - Npar)) || (dchisq <= 0.0)); i++) {
+    chisq = mrqmin (theta, pos, dpos, Npts, par, Npar, fellipseOD, !Quiet);
+    dchisq = ochisq - chisq;
+    ochisq = chisq;
+    if (!Quiet) gprint (GP_ERR, "dchisq: %f, Ndof: %d\n", dchisq, Npts - Npar);
+  }  
+  if (!Quiet) gprint (GP_ERR, "%d iterations\n", i); 
+
+  ResetVector (Xvec, OPIHI_FLT, xvec[0].Nelements);
+  ResetVector (Yvec, OPIHI_FLT, xvec[0].Nelements);
+  for (i = 0; i < xvec[0].Nelements; i++) {
+    Xvec[0].elements.Flt[i] = fellipseOD (theta[2*i+0], par, Npar, dpos);
+    Yvec[0].elements.Flt[i] = fellipseOD (theta[2*i+1], par, Npar, dpos);
+  }
+  Xvec[0].Nelements = xvec[0].Nelements;
+  Yvec[0].Nelements = xvec[0].Nelements;
+
+  covar = mrqcovar (Npar);
+  //set_variable ("E_X0",   par[PAR_X0]);
+  //set_variable ("E_Y0",   par[PAR_Y0]);
+  //set_variable ("E_RMAJ", par[PAR_RX]);
+  //set_variable ("E_RMIN", par[PAR_RY]);
+  //set_variable ("E_T0",   par[PAR_T0]);
+  //set_variable ("E_P0",   par[PAR_P0]);
+
+  set_variable ("E_RMAJ", par[PAR_RMIN] / par[PAR_EPSILON]);
+  set_variable ("E_RMIN", par[PAR_RMIN]);
+  set_variable ("E_PHI",  DEG_RAD*par[PAR_PHI]);
+
+  // set_variable ("dE_X0",   sqrt(covar[PAR_X0][PAR_X0]));
+  // set_variable ("dE_Y0",   sqrt(covar[PAR_Y0][PAR_Y0]));
+  // set_variable ("dE_RMAJ", sqrt(covar[PAR_RX][PAR_RX]));
+  // set_variable ("dE_RMIN", sqrt(covar[PAR_RY][PAR_RY]));
+  // set_variable ("dE_T0",   sqrt(covar[PAR_T0][PAR_T0]));
+  // set_variable ("dE_P0",   sqrt(covar[PAR_P0][PAR_P0]));
+
+  set_variable ("dE_RMAJ", sqrt(covar[PAR_EPSILON][PAR_EPSILON]*covar[PAR_RMIN][PAR_RMIN]));
+  set_variable ("dE_RMIN", sqrt(covar[PAR_RMIN][PAR_RMIN]));
+  set_variable ("dE_PHI",  DEG_RAD*sqrt(covar[PAR_PHI][PAR_PHI]));
+
+  // if (!Quiet) gprint (GP_ERR, "Xo : %f +/- %f\n", par[PAR_X0], sqrt(covar[PAR_X0][PAR_X0]));
+  // if (!Quiet) gprint (GP_ERR, "Yo : %f +/- %f\n", par[PAR_Y0], sqrt(covar[PAR_Y0][PAR_Y0]));
+  // if (!Quiet) gprint (GP_ERR, "RX : %f +/- %f\n", par[PAR_RX], sqrt(covar[PAR_RX][PAR_RX]));
+  // if (!Quiet) gprint (GP_ERR, "RY : %f +/- %f\n", par[PAR_RY], sqrt(covar[PAR_RY][PAR_RY]));
+  // if (!Quiet) gprint (GP_ERR, "To : %f +/- %f\n", par[PAR_T0], sqrt(covar[PAR_T0][PAR_T0]));
+  // if (!Quiet) gprint (GP_ERR, "Po : %f +/- %f\n", par[PAR_P0], sqrt(covar[PAR_P0][PAR_P0]));
+
+  if (!Quiet) gprint (GP_ERR, "RMAJ : %f +/- %f\n", par[PAR_RMIN] / par[PAR_EPSILON], 0.0);
+  if (!Quiet) gprint (GP_ERR, "RMIN : %f +/- %f\n", par[PAR_RMIN], sqrt(covar[PAR_RMIN][PAR_RMIN]));
+  if (!Quiet) gprint (GP_ERR, "PHI  : %f +/- %f\n", DEG_RAD*par[PAR_PHI],  sqrt(covar[PAR_PHI][PAR_PHI]));
+
+  free (pos);
+  free (dpos);
+  free (theta);
+
+  mrqfree (Npar);
+  return (TRUE);
+}
+
+# if (0)
+/**
+ * the full chisq is built of two associated sums over coordinates:
+ * chisq = sum ((X_obs - X_fit(t))^2 + (Y_obs - Y_fit(t))^2)
+ * 
+ * the independent variable is Theta:
+ * X_fit = X0 + RX*cos(theta)*cos(phi) - RY*sin(theta)*sin(phi)
+ * Y_fit = Y0 + RX*cos(theta)*sin(phi) + RY*sin(theta)*cos(phi)
+ *
+ * alternating calls to fellipseOD refer alternatively to X or Y
+ *
+ * parameters: Xo, Yo (center), RX, RY (axes), To (angle), Po (phase)
+ */
+
+opihi_flt fellipseOD (opihi_flt theta, opihi_flt *par, int Npar, opihi_flt *dpar) {
+  
+  static int pass = 0;
+
+  opihi_flt csphi = cos(par[PAR_P0]);
+  opihi_flt snphi = sin(par[PAR_P0]);
+
+  opihi_flt dtheta = theta - par[PAR_T0];
+  opihi_flt cstht = cos(dtheta);
+  opihi_flt sntht = sin(dtheta);
+
+  // value is X
+  if (pass == 0) {
+    pass = 1;
+
+    opihi_flt value = par[PAR_X0] + par[PAR_RX]*cstht*csphi - par[PAR_RY]*sntht*snphi;
+
+    dpar[PAR_X0] = 1.0;
+    dpar[PAR_Y0] = 0.0;
+    dpar[PAR_RX] = +cstht*csphi;
+    dpar[PAR_RY] = -sntht*snphi;
+    dpar[PAR_P0] = -par[PAR_RX]*cstht*snphi - par[PAR_RY]*sntht*csphi;
+    dpar[PAR_T0] = -par[PAR_RX]*sntht*csphi - par[PAR_RY]*cstht*snphi;
+
+    return (value);
+  }  
+
+  // value is Y
+  if (pass == 1) {
+    pass = 0;
+
+    opihi_flt value = par[PAR_Y0] + par[PAR_RX]*cstht*snphi + par[PAR_RY]*sntht*csphi;
+
+    dpar[PAR_X0]  = 0.0;
+    dpar[PAR_Y0]  = 1.0;
+    dpar[PAR_RX]  = +cstht*snphi;
+    dpar[PAR_RY]  = +sntht*csphi;
+    dpar[PAR_P0]  = +par[PAR_RX]*cstht*csphi - par[PAR_RY]*sntht*snphi;
+    dpar[PAR_T0]  = -par[PAR_RX]*sntht*snphi + par[PAR_RY]*cstht*csphi;
+    return (value);
+  }  
+  abort();
+}
+
+# endif
+
+/**
+ * the full chisq is built of two associated sums over coordinates:
+ * chisq = sum ((X_obs - X_fit(t))^2 + (Y_obs - Y_fit(t))^2)
+ * 
+ * the independent variable is theta (angle from 0,0 to point):
+ * X_fit = RMIN*R*cos(theta)
+ * Y_fit = RMIN*R*sin(theta)
+ * R = 1 / sqrt(sin^2(theta-phi) + epsilon^2 * cos^2(theta-phi))
+ *
+ * alternating calls to fellipseOD refer alternatively to X or Y
+ *
+ * RMIN, EPSILON, PHI
+ */
+
+opihi_flt fellipseOD (opihi_flt alpha, opihi_flt *par, int Npar, opihi_flt *dpar) {
+  
+  static int pass = 0;
+
+  opihi_flt cs_alpha = cos(alpha);
+  opihi_flt sn_alpha = sin(alpha);
+
+  opihi_flt cs_phi = cos(alpha - par[PAR_PHI]);
+  opihi_flt sn_phi = sin(alpha - par[PAR_PHI]);
+
+  opihi_flt r     = 1.0 / sqrt(SQ(sn_phi) + SQ(par[PAR_EPSILON]*cs_phi));
+  opihi_flt r3    = pow(r, 3.0);
+  opihi_flt drdE  = -0.5 * r3 * SQ(cs_phi) * 2.0 * par[PAR_EPSILON];
+  opihi_flt drdP  = -0.5 * r3 * (SQ(par[PAR_EPSILON]) - 1) * 2.0 * cs_phi * sn_phi;
+
+  // value is X
+  if (pass == 0) {
+    pass = 1;
+
+    opihi_flt value = par[PAR_RMIN]*cs_alpha*r;
+
+    dpar[PAR_RMIN]    = r*cs_alpha;
+    dpar[PAR_EPSILON] = par[PAR_RMIN]*cs_alpha*drdE;
+    dpar[PAR_PHI]     = 4.0*par[PAR_RMIN]*cs_alpha*drdP;
+
+    return (value);
+  }  
+
+  // value is Y
+  if (pass == 1) {
+    pass = 0;
+
+    opihi_flt value = par[PAR_RMIN]*sn_alpha*r;
+
+    dpar[PAR_RMIN]    = r*sn_alpha;
+    dpar[PAR_EPSILON] = par[PAR_RMIN]*sn_alpha*drdE;
+    dpar[PAR_PHI]     = 4.0*par[PAR_RMIN]*sn_alpha*drdP;
+
+    return (value);
+  }  
+  abort();
+}
+
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/vgroup.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/vgroup.c	(revision 24992)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.data/vgroup.c	(revision 24992)
@@ -0,0 +1,71 @@
+# include "data.h"
+
+int vgroup (int argc, char **argv) {
+  
+  int i, j, N;
+  Vector *xin, *yin, *xout, *yout;
+  opihi_flt xmin, xmax, sum;
+  double *values;
+
+  // Normalize = FALSE;
+  // if ((N = get_argument (argc, argv, "-norm"))) {
+  //   remove_argument (N, &argc, argv);
+  //   Normalize = TRUE;
+  // }
+  // 
+  // Ignore = FALSE;
+  // IgnoreValue = 0.0;
+  // if ((N = get_argument (argc, argv, "-ignore"))) {
+  //   Ignore = TRUE;
+  //   remove_argument (N, &argc, argv);
+  //   IgnoreValue = atof (argv[N]);
+  //   remove_argument (N, &argc, argv);
+  // }
+
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: vbin <xin> <yin> <xout> <yout>\n");
+    gprint (GP_ERR, " group x,y values in bins defined by <xout>\n");
+    return (FALSE);
+  }
+
+  if ((xin  = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((yin  = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((xout = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((yout = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  // re-binning creates a float vector
+  ResetVector (yout, OPIHI_FLT, xout[0].Nelements);
+
+  // storage vector
+  ALLOCATE (values, double, xin[0].Nelements);
+
+  for (i = 0; i < xout[0].Nelements - 1; i++) {
+    xmin = xout[0].elements.Flt[i];
+    xmax = xout[0].elements.Flt[i+1];
+
+    N = 0;
+    for (j = 0; j < xin[0].Nelements; j++) {
+      if (xin[0].elements.Flt[j] < xmin) continue;
+      if (xin[0].elements.Flt[j] > xmax) continue;
+      values[N] = yin[0].elements.Flt[j];
+      N++;
+    }
+    
+    dsort (values, N);
+    if (N > 1) {
+      sum = (N % 2) ? values[(int)(0.5*N)] : 0.5*(values[N/2] + values[N/2 + 1]);
+    } else {
+      sum = values[0];
+    }
+
+    // measure the stat for this bin
+    // sum = 0.0;
+    // for (j = 0; j < N; j++) {
+    //   sum += values[j];
+    // }
+    // sum /= N;
+    yout[0].elements.Flt[i] = sum;
+  }
+  free (values);
+  return (TRUE);
+}