Index: trunk/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 24027)
+++ trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 25757)
@@ -21,4 +21,5 @@
 $(SRC)/cgrid.$(ARCH).o		   \
 $(SRC)/coords.$(ARCH).o	   \
+$(SRC)/cdot.$(ARCH).o		   \
 $(SRC)/cplot.$(ARCH).o		   \
 $(SRC)/csystem.$(ARCH).o	   \
@@ -36,4 +37,11 @@
 $(SRC)/medianmap.$(ARCH).o	   \
 $(SRC)/mkgauss.$(ARCH).o	   \
+$(SRC)/mksersic.$(ARCH).o	   \
+$(SRC)/galradius.$(ARCH).o	   \
+$(SRC)/galradbins.$(ARCH).o	   \
+$(SRC)/galsectors.$(ARCH).o	   \
+$(SRC)/galprofiles.$(ARCH).o	   \
+$(SRC)/elliprofile.$(ARCH).o	   \
+$(SRC)/petrosian.$(ARCH).o	   \
 $(SRC)/multifit.$(ARCH).o	   \
 $(SRC)/objload.$(ARCH).o	   \
Index: trunk/Ohana/src/opihi/cmd.astro/cdot.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/cdot.c	(revision 25757)
+++ trunk/Ohana/src/opihi/cmd.astro/cdot.c	(revision 25757)
@@ -0,0 +1,38 @@
+# include "data.h"
+
+int cdot (int argc, char **argv) {
+  
+  int kapa, status;
+  Graphdata graphmode;
+  float x, y;
+  double r, d, Rmin, Rmax;
+
+  if (!style_args (&graphmode, &argc, argv, &kapa)) return FALSE;
+
+  if (argc != 3) {
+    gprint (GP_ERR, "USAGE: dot <ra> <dec>\n");
+    return (FALSE);
+  }
+  r = atof(argv[1]);
+  d = atof(argv[2]);
+
+  Rmin = graphmode.coords.crval1 - 182.0;
+  Rmax = graphmode.coords.crval1 + 182.0;
+  // Rmid = 0.5*(Rmin + Rmax);
+
+  /* set point style and errorbar mode (these are NOT sticky) */
+  graphmode.style = 2;
+  graphmode.etype = 0;
+
+  r = ohana_normalize_angle (r);
+  while (r < Rmin) r += 360.0;
+  while (r > Rmax) r -= 360.0;
+
+  status = fRD_to_XY (&x, &y, r, d, &graphmode.coords);
+
+  if (!KapaPrepPlot (kapa, 1, &graphmode)) return (FALSE);
+  KapaPlotVector (kapa, 1, &x, "x");
+  KapaPlotVector (kapa, 1, &y, "y");
+  
+  return (TRUE);
+}
Index: trunk/Ohana/src/opihi/cmd.astro/cgrid.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/cgrid.c	(revision 24027)
+++ trunk/Ohana/src/opihi/cmd.astro/cgrid.c	(revision 25757)
@@ -37,5 +37,5 @@
   /* set spacings for RA */
   minorRA = minorDEC = 0.1;
-  range = MIN (fabs(graphmode.xmax-graphmode.xmin), fabs(graphmode.ymax-graphmode.ymin));
+  range = MIN (fabs(graphmode.coords.cdelt1*(graphmode.xmax-graphmode.xmin)), fabs(graphmode.coords.cdelt2*(graphmode.ymax-graphmode.ymin)));
   if (NorthPole || SouthPole) range = 360;
   lrange = log10(MAX(fabs(range), 1e-30));
@@ -58,6 +58,7 @@
   }
   dR = range / 100.0;
+
   /* set spacings for DEC */
-  range = MIN (fabs(graphmode.xmax-graphmode.xmin), fabs(graphmode.ymax-graphmode.ymin));
+  range = MIN (fabs(graphmode.coords.cdelt1*(graphmode.xmax-graphmode.xmin)), fabs(graphmode.coords.cdelt2*(graphmode.ymax-graphmode.ymin)));
   lrange = log10(MAX(fabs(range), 1e-30));
   factor = (int) (lrange);
Index: trunk/Ohana/src/opihi/cmd.astro/elliprofile.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/elliprofile.c	(revision 25757)
+++ trunk/Ohana/src/opihi/cmd.astro/elliprofile.c	(revision 25757)
@@ -0,0 +1,64 @@
+# include "astro.h"
+
+int elliprofile (int argc, char **argv) {
+  
+  int i, ix, iy, Nx, Ny, N, center;
+  float *in;
+  double Rmaj, Rmin, phi, alpha, Io, ARatio;
+  double root1, root2, R, A1, A2, A3;
+  double Sx, Sy, Sxy;
+  double x, y, r, f, Xo, Yo;
+  Buffer *buf;
+  Vector *rvec, *fvec;
+
+  if (argc != 9) {
+    gprint (GP_ERR, "USAGE: elliprofile (buffer) (rvec) (fvec) (Xo) (Yo) (Rmaj) (Rmin) (phi)\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 ((rvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((fvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  ResetVector (rvec, OPIHI_FLT, Nx*Ny);
+  ResetVector (fvec, OPIHI_FLT, Nx*Ny);
+
+  Xo = atof(argv[4]);
+  Yo = atof(argv[5]);
+
+  /* shape parameters */
+  Rmaj = atof (argv[6]);
+  Rmin = atof (argv[7]);
+  phi  = atof (argv[8]);
+
+  /* given Rmaj, Rmin, phi, 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*phi);
+  A2 = 0.25*(root1 + root2) + 0.5*R*cos(2*RAD_DEG*phi);
+  A3 = -R*sin(2*RAD_DEG*phi);
+
+  Sx = 0.5/A1;
+  Sy = 0.5/A2;
+  Sxy = A3;
+
+  /* f = exp (-r^alpha), r^2 = (x^2 / 2Sx) + (y^2 / 2Sy) + Sxy*x*y */
+  in = (float *) buf[0].matrix.buffer;
+  for (i = iy = 0; iy < Ny; iy++) {
+      for (ix = 0; ix < Nx; ix++, in++, i++) {
+      x = ix + 0.5 - Xo;
+      y = iy + 0.5 - Yo;
+      rvec[0].elements.Flt[i] = sqrt(0.5*x*x/Sx + 0.5*y*y/Sy + x*y*Sxy);
+      fvec[0].elements.Flt[i] = *in;
+    }
+  }
+  rvec[0].Nelements = fvec[0].Nelements = i;
+  return (TRUE);
+}
Index: trunk/Ohana/src/opihi/cmd.astro/galcontour.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/galcontour.c	(revision 25757)
+++ trunk/Ohana/src/opihi/cmd.astro/galcontour.c	(revision 25757)
@@ -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: trunk/Ohana/src/opihi/cmd.astro/galprofiles.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/galprofiles.c	(revision 25757)
+++ trunk/Ohana/src/opihi/cmd.astro/galprofiles.c	(revision 25757)
@@ -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: trunk/Ohana/src/opihi/cmd.astro/galradbins.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/galradbins.c	(revision 25757)
+++ trunk/Ohana/src/opihi/cmd.astro/galradbins.c	(revision 25757)
@@ -0,0 +1,84 @@
+# include "astro.h"
+
+int MeanSurfaceBrightness (Vector *rvec, Vector *fvec, Vector *Rvec, Vector *Fvec, Vector *Avec, float Rmin, float Rmax, int bin);
+
+// given a collection of r, f points sampled at pixels, generate a pair of vectors r, f
+// where f is defined as the mean surface brightness for \alpha r_i < r < \beta r_i.
+// sample r at r_i = i
+
+// this function can be much more efficient if the input vectors are sorted by R and that
+// fact is used in generating the radial bins...
+int galradbins (int argc, char **argv) {
+  
+  int i, Nbin;
+  double Rmin, Rmax;
+  Vector *fvec, *rvec, *Fvec, *Rvec, *Avec;
+
+  if (argc != 6) {
+    gprint (GP_ERR, "USAGE: galradius (r_in) (f_in) (r_out) (f_out) (area_out)\n");
+    return (FALSE);
+  }
+
+  /* select input / output buffers */
+  if ((rvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((fvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Rvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Fvec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Avec = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  Rmax = rvec[0].elements.Flt[0];
+  for (i = 0; i < rvec[0].Nelements; i++) {
+    Rmax = MAX(Rmax, rvec[0].elements.Flt[i]);
+  }
+  Nbin = 0.8*Rmax + 2;
+  ResetVector (Rvec, OPIHI_FLT, Nbin);
+  ResetVector (Fvec, OPIHI_FLT, Nbin);
+  ResetVector (Avec, OPIHI_FLT, Nbin);
+
+  // the first three bins are specially defined:
+
+  // 0 : r < 1.0 -- Area is pi
+  MeanSurfaceBrightness (rvec, fvec, Rvec, Fvec, Avec, 0.00, 1.00, 0);
+
+  // area is \int 2 \pi r dr = \pi (r2^2 - r1^2)
+
+  // 1 : 1.0 < r < 1.25 -- Area is \pi (1.25^2 - 1^2)
+  MeanSurfaceBrightness (rvec, fvec, Rvec, Fvec, Avec, 1.00, 1.25, 1);
+
+  // 2 : 1.25 < r < 1.60 -- Area is \pi (1.6^2 - 1.25^2)
+  MeanSurfaceBrightness (rvec, fvec, Rvec, Fvec, Avec, 1.25, 1.60, 2);
+
+  // from bin 3 on out, r = i - 1 
+  for (i = 3; i < Nbin; i++) {
+    MeanSurfaceBrightness (rvec, fvec, Rvec, Fvec, Avec, 0.8*(i-1), 1.25*(i-1), i);
+  }
+
+  return (TRUE);
+}
+
+int MeanSurfaceBrightness (Vector *rvec, Vector *fvec, Vector *Rvec, Vector *Fvec, Vector *Avec, float Rmin, float Rmax, int bin) {
+
+  int i, Npts;
+  double Fsum;
+
+  Fsum = Npts = 0;
+
+  for (i = 0; i < rvec[0].Nelements; i++) {
+    if (rvec[0].elements.Flt[i] < Rmin) continue;
+    if (rvec[0].elements.Flt[i] > Rmax) continue;
+    Fsum += fvec[0].elements.Flt[i];
+    Npts ++;
+  }
+  if (Rmin > 0.0) {
+    Rvec[0].elements.Flt[bin] = sqrt(Rmin * Rmax); // XXX what is the correct mean radius?
+  } else {
+    Rvec[0].elements.Flt[bin] = 0.5 * Rmax; // XXX what is the correct mean radius?
+  }
+  if (Npts > 0) {
+    Fvec[0].elements.Flt[bin] = Fsum / Npts; // XXX what is the correct mean radius?
+  } else {
+    Fvec[0].elements.Flt[bin] = NAN; // XXX what is the correct mean radius?
+  }
+  Avec[0].elements.Flt[bin] = M_PI * (SQ(Rmax) - SQ(Rmin));
+  return (TRUE);
+}
Index: trunk/Ohana/src/opihi/cmd.astro/galradius.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/galradius.c	(revision 25757)
+++ trunk/Ohana/src/opihi/cmd.astro/galradius.c	(revision 25757)
@@ -0,0 +1,123 @@
+# include "astro.h"
+
+// given an image and associated galaxy profile, find:
+// a) a characteristic flux level that defines a useful isophot
+// b) radius at which the profile crosses that isophot
+int galradius (int argc, char **argv) {
+  
+  int i, j, Nx, Ny, N, Nsec, NELEMENTS, Rbin, Nout, above;
+  double theta, dtheta, x, y, Xo, Yo, r, Rmax, Rmin, value;
+  double Fmin, Fmax, dF, Fm, Fp, Fo, Rsum, Rnpt, Ro;
+  opihi_flt *flux, *radius, *values;
+  Buffer *buf;
+  Vector *fvec, *rvec, *tvec, *Fvec, *Rvec;
+  char name[128];
+
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: galradius (radius) (flux) (min_flux) (max_flux)\n");
+    return (FALSE);
+  }
+
+  /* select input / output buffers */
+  if ((rvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((fvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+
+  if ((Rvec = SelectVector ("rg", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Fvec = SelectVector ("fg", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  Fmin = atof(argv[3]);
+  Fmax = atof(argv[4]);
+
+  // fvec is a noise sample of the galaxy radial profile at points rvec
+  // rebin fvec into samples defined by the isophot 
+
+  // base selections on fluxes defined by the flux range dF
+  dF = Fmax - Fmin;
+
+  // select all points with flux in the range Fmin + 0.25*dF to Fmin + 0.75*dF
+  Fm = Fmin + 0.25*dF;
+  Fp = Fmin + 0.75*dF;
+  Fo = Fmin + 0.50*dF;
+      
+  Rsum = 0;
+  Rnpt = 0;
+  for (i = 0; i < fvec[0].Nelements; i++) {
+    if (fvec[0].elements.Flt[i] < Fm) continue;
+    if (fvec[0].elements.Flt[i] > Fp) continue;
+      
+    Rsum += rvec[0].elements.Flt[i];
+    Rnpt ++;
+  }
+
+  // determine the binned radius values
+  if (Rnpt == 0) {
+    Rbin = 1;
+  } else {
+    Rbin = MAX(1, 0.5*(Rsum / Rnpt));
+  }
+
+  // do not bother rebinning if the bin size is only 2 or less
+  if (Rbin <= 2) {
+    Nout   = fvec[0].Nelements;
+    ResetVector (Rvec, OPIHI_FLT, Nout);
+    ResetVector (Fvec, OPIHI_FLT, Nout);
+    memcpy (Fvec[0].elements.Flt, fvec[0].elements.Flt, Fvec[0].Nelements*sizeof(opihi_flt));
+    memcpy (Rvec[0].elements.Flt, rvec[0].elements.Flt, Rvec[0].Nelements*sizeof(opihi_flt));
+    flux   = fvec[0].elements.Flt;
+    radius = rvec[0].elements.Flt;
+  } else {
+    // rebin the vectors by Rbin values:
+    Nout = fvec[0].Nelements / Rbin + 1;
+    ALLOCATE (values, opihi_flt, fvec[0].Nelements);
+  
+    ResetVector (Rvec, OPIHI_FLT, Nout);
+    ResetVector (Fvec, OPIHI_FLT, Nout);
+    flux   = Fvec[0].elements.Flt;
+    radius = Rvec[0].elements.Flt;
+
+    for (i = 0; i < Nout; i++) {
+      radius[i] = (i + 0.5)*Rbin;
+      Rmin = radius[i] - 0.5*Rbin;
+      Rmax = radius[i] + 0.5*Rbin;
+
+      N = 0;
+      for (j = 0; j < fvec[0].Nelements; j++) {
+	if (rvec[0].elements.Flt[j] < Rmin) continue;
+	if (rvec[0].elements.Flt[j] > Rmax) continue;
+	values[N] = fvec[0].elements.Flt[j];
+	N++;
+      }
+
+      // take the median of the values
+      dsort (values, N);
+      if (N > 1) {
+	flux[i] = (N % 2) ? values[(int)(0.5*N)] : 0.5*(values[N/2] + values[N/2 + 1]);
+      } else {
+	flux[i] = values[0];
+      }
+    }
+    free (values);
+  }
+
+  above = TRUE;
+  for (i = 0; i < Nout; i++) {
+
+    // find the largest radius that matches the flux transition
+    if (above && (flux[i] < Fo)) {
+      if (i == 0) { 
+	// assume Fmax @ R = 0.0
+	Ro = radius[i] * (Fo - Fmax) / (flux[i] - Fmax);
+      } else {
+	Ro = radius[i-1] + (radius[i] - radius[i-1]) * (Fo - flux[i-1]) / (flux[i] - flux[i-1]);
+      }
+      above = FALSE;
+    }
+  
+    if (!above && (flux[i] >= Fo)) {
+      above = TRUE;
+    }
+  }
+
+  set_variable ("G_R50", Ro);
+  return (TRUE);
+}
Index: trunk/Ohana/src/opihi/cmd.astro/galsectors.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/galsectors.c	(revision 25757)
+++ trunk/Ohana/src/opihi/cmd.astro/galsectors.c	(revision 25757)
@@ -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: trunk/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 24027)
+++ trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 25757)
@@ -5,4 +5,5 @@
 int cgrid                   PROTO((int, char **));
 int coords                  PROTO((int, char **));
+int cdot                    PROTO((int, char **));
 int cplot                   PROTO((int, char **));
 int csystem                 PROTO((int, char **));
@@ -22,5 +23,12 @@
 int imsub                   PROTO((int, char **));
 int medianmap               PROTO((int, char **));
+int galsectors              PROTO((int, char **));
+int galprofiles             PROTO((int, char **));
+int galradius               PROTO((int, char **));
+int galradbins              PROTO((int, char **));
+int elliprofile             PROTO((int, char **));
+int petrosian               PROTO((int, char **));
 int mkgauss                 PROTO((int, char **));
+int mksersic                PROTO((int, char **));
 int multifit                PROTO((int, char **));
 int objload                 PROTO((int, char **));
@@ -44,4 +52,5 @@
   {1, "cgrid",       cgrid,        "plot sky coordinate grid"},
   {1, "coords",      coords,       "load coordinates for buffer from file"},
+  {1, "cdot",        cdot,         "plot point in sky coordinates"},
   {1, "cplot",       cplot,        "plot vectors in sky coordinates"},
   {1, "csystem",     csystem,      "convert between coordinate systems"},
@@ -61,4 +70,11 @@
   {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, "galradius",   galradius,    "generate radial vectors with interpolation along paths"},
+  {1, "galradbins",  galradbins,   "generate radial vectors with interpolation along paths"},
+  {1, "elliprofile", elliprofile,  "generate radial vectors with interpolation along paths"},
+  {1, "petrosian",   petrosian,    "petrosian parameters given radial bins"},
   {1, "multifit",    multifit,     "fit multi-order spectrum"},
   {1, "objload",     objload,      "plot obj data on Ximage "},
Index: trunk/Ohana/src/opihi/cmd.astro/mksersic.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/mksersic.c	(revision 25757)
+++ trunk/Ohana/src/opihi/cmd.astro/mksersic.c	(revision 25757)
@@ -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: trunk/Ohana/src/opihi/cmd.astro/petrosian.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/petrosian.c	(revision 25757)
+++ trunk/Ohana/src/opihi/cmd.astro/petrosian.c	(revision 25757)
@@ -0,0 +1,90 @@
+# include "astro.h"
+
+int petrosian (int argc, char **argv) {
+  
+  int i, above;
+  float *in;
+  double Fsum, Asum, Area, R_90, rad_90, flux_90;
+  Vector *rvec, *avec, *fvec, *Rvec, *Fvec, *Svec, *Avec;
+
+  if (argc != 8) {
+    gprint (GP_ERR, "USAGE: petrosian (radius) (area) (flux) (P_ratio) (P_flux) (P_sb) (P_area)\n");
+    return (FALSE);
+  }
+
+  /* select input vectors */
+  if ((rvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((avec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((fvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Rvec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Fvec = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Svec = SelectVector (argv[6], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((Avec = SelectVector (argv[7], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  // difficult work goes into cleaning the input galaxy profile:
+  // a) extract radial profiles applying 180 degree symmetry requirement
+  // b) determine the axial ratio and position angle (at r_50)
+  // c) extract normalized radial profiles applying outlier clipping
+  // 
+  // the resulting vectors r, f represent the mean surface brightness at a radius r_i, 
+  // measured between \alpha*r_i and \beta*r_i: f(r_i) = \sum_\alpha*r_i^\beta*r_i flux(r)
+
+  // given a clean set of vectors representing the f_i, r_i, determine F_P, R_P such that:
+  // f(r_i) / \sum_0^r_i f(r) = f_P 
+
+  // generate a vector representing f_P(r_i):
+
+  ResetVector (Fvec, OPIHI_FLT, fvec[0].Nelements);
+  ResetVector (Rvec, OPIHI_FLT, fvec[0].Nelements);
+  ResetVector (Svec, OPIHI_FLT, fvec[0].Nelements);
+  ResetVector (Avec, OPIHI_FLT, fvec[0].Nelements);
+
+  above = TRUE;
+  Fsum = 0.0;
+  Asum = 0.0;
+  Area = avec[0].elements.Flt[0];
+  for (i = 0; i < fvec[0].Nelements; i++) {
+    // for nan bins, we keep the area for use with the next valid bin
+    if (isnan(fvec[0].elements.Flt[i])) {
+      Area += avec[0].elements.Flt[i];
+      continue;
+    } 
+    Fsum += fvec[0].elements.Flt[i] * Area;
+    Asum += Area;
+    if (i+1 < fvec[0].Nelements) {
+      Area = avec[0].elements.Flt[i+1];
+    }
+
+    Rvec[0].elements.Flt[i] = Asum * fvec[0].elements.Flt[i] / Fsum;
+    Fvec[0].elements.Flt[i] = Fsum;
+    Svec[0].elements.Flt[i] = Fsum / Asum;
+    Avec[0].elements.Flt[i] = Asum;
+
+    // anytime we transition below the petrosian ratio R_90, calculate the radius and flux
+    // we will keep and report the last (largest radius) value
+    if (above && (Rvec[0].elements.Flt[i] < R_90)) {
+      // interpolate Rvec between i-1 and i to R_90 to get flux (Fvec) and radius (rvec)
+
+      if (i == 0) { 
+	// assume Fmax @ R = 0.0
+	rad_90  = rvec[0].elements.Flt[i] * (R_90 - 1.0) / (Rvec[0].elements.Flt[i] - 1.0);
+	flux_90 = Fvec[0].elements.Flt[i] * (R_90 - 1.0) / (Rvec[0].elements.Flt[i] - 1.0);
+      } else {
+	rad_90  = rvec[0].elements.Flt[i-1] + (rvec[0].elements.Flt[i] - rvec[0].elements.Flt[i-1]) * (R_90 - Rvec[0].elements.Flt[i-1]) / (Rvec[0].elements.Flt[i] - Rvec[0].elements.Flt[i-1]);
+	flux_90 = Fvec[0].elements.Flt[i-1] + (Fvec[0].elements.Flt[i] - Fvec[0].elements.Flt[i-1]) * (R_90 - Rvec[0].elements.Flt[i-1]) / (Rvec[0].elements.Flt[i] - Rvec[0].elements.Flt[i-1]);
+      }
+
+      above = FALSE;
+    }
+    
+    // reset on transitions up, but do not re-calculate rad_90, flux_90
+    if (!above && (Rvec[0].elements.Flt[i] >= R_90)) {
+      above = TRUE;
+    }
+  }
+
+  set_variable ("P_R90", rad_90);
+  set_variable ("P_F90", flux_90);
+
+  return (TRUE);
+}
Index: trunk/Ohana/src/opihi/cmd.astro/region.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/region.c	(revision 24027)
+++ trunk/Ohana/src/opihi/cmd.astro/region.c	(revision 25757)
@@ -6,5 +6,5 @@
   double Ra, Dec, Radius;
   float dx, dy;
-  int N, kapa;
+  int N, kapa, NoClear;
   char *name;
   Graphdata graphmode;
@@ -46,4 +46,10 @@
     remove_argument (N, &argc, argv);
     graphmode.flipnorth = FALSE;
+  }
+
+  NoClear = FALSE;
+  if ((N = get_argument (argc, argv, "-no-clear"))) {
+    remove_argument (N, &argc, argv);
+    NoClear = TRUE;
   }
 
@@ -122,10 +128,10 @@
   graphmode.coords.cdelt1 = graphmode.coords.cdelt2 = 1.0;
 
-  KapaClearSections (kapa);
+  if (!NoClear) KapaClearSections (kapa);
   KapaSetLimits (kapa, &graphmode);
 
   /* drop this? */
-  sprintf (string, "%8.4f %8.4f (%f)", Ra, Dec, Radius);
-  KapaSendLabel (kapa, string, 2);
+  // sprintf (string, "%8.4f %8.4f (%f)", Ra, Dec, Radius);
+  // KapaSendLabel (kapa, string, 2);
 
   // XXX is this the right thing to be doing?