Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/Makefile	(revision 25017)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/Makefile	(revision 25018)
@@ -37,6 +37,9 @@
 $(SRC)/mkgauss.$(ARCH).o	   \
 $(SRC)/mksersic.$(ARCH).o	   \
+$(SRC)/galradius.$(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: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/elliprofile.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/elliprofile.c	(revision 25018)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/elliprofile.c	(revision 25018)
@@ -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: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galradius.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galradius.c	(revision 25018)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galradius.c	(revision 25018)
@@ -0,0 +1,112 @@
+# 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;
+  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);
+  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, Rsum / Rnpt);
+  }
+
+  // do not bother rebinning if the bin size is only 2 or less
+  if (Rbin <= 2) {
+    Nout   = fvec[0].Nelements;
+    flux   = fvec[0].elements.Flt;
+    radius = rvec[0].elements.Flt;
+    // XXX handle the necessary free below
+  } else {
+    // rebin the vectors by Rbin values:
+    Nout = fvec[0].Nelements / Rbin + 1;
+    ALLOCATE (flux,   opihi_flt, Nout);
+    ALLOCATE (radius, opihi_flt, Nout);
+    ALLOCATE (values, opihi_flt, fvec[0].Nelements);
+  
+    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[i] < Rmin) continue;
+	if (rvec[0].elements.Flt[i] > Rmax) continue;
+	values[N] = fvec[0].elements.Flt[i];
+	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];
+      }
+    }
+  }
+
+  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: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/init.c	(revision 25017)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/init.c	(revision 25018)
@@ -24,4 +24,7 @@
 int galsectors              PROTO((int, char **));
 int galprofiles             PROTO((int, char **));
+int galradius               PROTO((int, char **));
+int elliprofile             PROTO((int, char **));
+int petrosian               PROTO((int, char **));
 int mkgauss                 PROTO((int, char **));
 int mksersic                PROTO((int, char **));
@@ -67,4 +70,7 @@
   {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, "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: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/petrosian.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/petrosian.c	(revision 25018)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/petrosian.c	(revision 25018)
@@ -0,0 +1,71 @@
+# include "astro.h"
+
+int petrosian (int argc, char **argv) {
+  
+  int i, above;
+  float *in;
+  double Fsum, R_90, rad_90, flux_90;
+  Vector *rvec, *fvec, *Rvec, *Fvec;
+
+  if (argc != 2) {
+    gprint (GP_ERR, "USAGE: petrosian (radius) (flux) (P_ratio) (P_flux)\n");
+    return (FALSE);
+  }
+
+  /* select input vectors */
+  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);
+
+  // 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);
+
+  above = TRUE;
+  Fsum = 0.0;
+  for (i = 0; i < fvec[0].Nelements; i++) {
+    Fsum += fvec[0].elements.Flt[i] * M_PI * SQ(rvec[0].elements.Flt[i]);
+    Rvec[0].elements.Flt[i] = Fsum / fvec[0].elements.Flt[i];
+    Fvec[0].elements.Flt[i] = Fsum;
+
+    // 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);
+}