Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/Makefile	(revision 25032)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/Makefile	(revision 25033)
@@ -38,4 +38,5 @@
 $(SRC)/mksersic.$(ARCH).o	   \
 $(SRC)/galradius.$(ARCH).o	   \
+$(SRC)/galradbins.$(ARCH).o	   \
 $(SRC)/galsectors.$(ARCH).o	   \
 $(SRC)/galprofiles.$(ARCH).o	   \
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galradbins.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galradbins.c	(revision 25033)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galradbins.c	(revision 25033)
@@ -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: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galradius.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galradius.c	(revision 25032)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/galradius.c	(revision 25033)
@@ -11,5 +11,5 @@
   opihi_flt *flux, *radius, *values;
   Buffer *buf;
-  Vector *fvec, *rvec, *tvec;
+  Vector *fvec, *rvec, *tvec, *Fvec, *Rvec;
   char name[128];
 
@@ -22,4 +22,8 @@
   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]);
@@ -50,5 +54,5 @@
     Rbin = 1;
   } else {
-    Rbin = MAX(1, Rsum / Rnpt);
+    Rbin = MAX(1, 0.5*(Rsum / Rnpt));
   }
 
@@ -56,14 +60,20 @@
   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;
-    // 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);
   
+    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;
@@ -73,7 +83,7 @@
       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];
+	if (rvec[0].elements.Flt[j] < Rmin) continue;
+	if (rvec[0].elements.Flt[j] > Rmax) continue;
+	values[N] = fvec[0].elements.Flt[j];
 	N++;
       }
@@ -87,4 +97,5 @@
       }
     }
+    free (values);
   }
 
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/init.c	(revision 25032)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/init.c	(revision 25033)
@@ -25,4 +25,5 @@
 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 **));
@@ -71,4 +72,5 @@
   {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"},
Index: /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/petrosian.c
===================================================================
--- /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/petrosian.c	(revision 25032)
+++ /branches/eam_branches/20090715/Ohana/src/opihi/cmd.astro/petrosian.c	(revision 25033)
@@ -5,9 +5,9 @@
   int i, above;
   float *in;
-  double Fsum, R_90, rad_90, flux_90;
-  Vector *rvec, *fvec, *Rvec, *Fvec;
+  double Fsum, Asum, Area, R_90, rad_90, flux_90;
+  Vector *rvec, *avec, *fvec, *Rvec, *Fvec, *Svec, *Avec;
 
-  if (argc != 2) {
-    gprint (GP_ERR, "USAGE: petrosian (radius) (flux) (P_ratio) (P_flux)\n");
+  if (argc != 8) {
+    gprint (GP_ERR, "USAGE: petrosian (radius) (area) (flux) (P_ratio) (P_flux) (P_sb) (P_area)\n");
     return (FALSE);
   }
@@ -15,7 +15,10 @@
   /* 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);
+  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:
@@ -34,11 +37,27 @@
   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++) {
-    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];
+    // 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