Index: branches/eam_branches/20090715/psphot/src/Makefile.am
===================================================================
--- branches/eam_branches/20090715/psphot/src/Makefile.am	(revision 25032)
+++ branches/eam_branches/20090715/psphot/src/Makefile.am	(revision 25105)
@@ -25,5 +25,5 @@
 libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 
-bin_PROGRAMS = psphot psphotTest psphotMomentsStudy
+bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy 
 
 psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
@@ -38,4 +38,8 @@
 psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 psphotMomentsStudy_LDADD = libpsphot.la
+
+psphotPetrosianStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psphotPetrosianStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+psphotPetrosianStudy_LDADD = libpsphot.la
 
 psphot_SOURCES = \
@@ -61,4 +65,15 @@
 psphotMomentsStudy_SOURCES = \
         psphotMomentsStudy.c
+
+psphotPetrosianStudy_SOURCES = \
+        psphotPetrosianStudy.c \
+        psphotPetrosianProfile.c \
+        psphotRadialProfileByAngle.c \
+	psphotRadiiFromProfiles.c \
+        psphotEllipticalContour.c \
+        psphotEllipticalProfile.c \
+        psphotPetrosianRadialBins.c \
+        psphotPetrosianVisual.c \
+        psphotPetrosianStats.c
 
 libpsphot_la_SOURCES = \
@@ -128,5 +143,6 @@
 	psphotCheckStarDistribution.c  \
 	psphotThreadTools.c  	       \
-	psphotAddNoise.c
+	psphotAddNoise.c               \
+	pmPetrosian.c
 
 # dropped? psphotGrowthCurve.c
@@ -134,4 +150,5 @@
 include_HEADERS = \
 	psphot.h \
+	pmPetrosian.h \
 	psphotErrorCodes.h
 
Index: branches/eam_branches/20090715/psphot/src/pmPetrosian.c
===================================================================
--- branches/eam_branches/20090715/psphot/src/pmPetrosian.c	(revision 25032)
+++ branches/eam_branches/20090715/psphot/src/pmPetrosian.c	(revision 25105)
@@ -11,13 +11,20 @@
 # include "psphotInternal.h"
 
-static void pmPetrosianFree(pmPetrosian *petro)
+static void pmPetrosianFree(pmPetrosian *petrosian)
 {
-   if (!petro) {
+    if (!petrosian) {
         return;
-   }
-   psFree(petrosian->radii);
-   psFree(petrosian->fluxes);
-   psFree(petrosian->theta);
-   psFree(petrosian->isophotalRadii);
+    }
+    psFree(petrosian->radii);
+    psFree(petrosian->fluxes);
+    psFree(petrosian->theta);
+    psFree(petrosian->isophotalRadii);
+
+    psFree(petrosian->radiusElliptical);
+    psFree(petrosian->fluxElliptical);
+
+    psFree(petrosian->binnedFlux);
+    psFree(petrosian->radialBins);
+    psFree(petrosian->area);
 }
 
@@ -32,7 +39,65 @@
     petrosian->isophotalRadii = NULL;
 
-    petrosian->axes = {0.0, 0.0, 0.0};
+    petrosian->radiusElliptical = NULL;
+    petrosian->fluxElliptical = NULL;
 
-    return(petrosian);
+    petrosian->radialBins = NULL;
+    petrosian->area = NULL;
+    petrosian->binnedFlux = NULL;
+
+    petrosian->petrosianRadius = NAN;
+    petrosian->petrosianFlux = NAN;
+
+    // petrosian->axes = {0.0, 0.0, 0.0};
+
+    return petrosian;
 }
 
+bool psphotPetrosianFreeVectors(pmPetrosian *petrosian) {
+
+    psFree(petrosian->radii);
+    psFree(petrosian->fluxes);
+    psFree(petrosian->theta);
+    psFree(petrosian->isophotalRadii);
+
+    psFree(petrosian->radiusElliptical);
+    psFree(petrosian->fluxElliptical);
+
+    psFree(petrosian->binnedFlux);
+    psFree(petrosian->radialBins);
+    psFree(petrosian->area);
+
+    petrosian->radii = NULL;
+    petrosian->fluxes = NULL;
+    petrosian->theta = NULL;
+    petrosian->isophotalRadii = NULL;
+
+    petrosian->radiusElliptical = NULL;
+    petrosian->fluxElliptical = NULL;
+
+    petrosian->radialBins = NULL;
+    petrosian->area = NULL;
+    petrosian->binnedFlux = NULL;
+    
+    return true;
+}
+
+# define COMPARE_INDEX(A,B) (index->data.F32[A] < index->data.F32[B])
+# define SWAP_INDEX(TYPE,A,B) { \
+  float tmp; \
+  if (A != B) { \
+    tmp = index->data.F32[A]; \
+    index->data.F32[A] = index->data.F32[B]; \
+    index->data.F32[B] = tmp; \
+    tmp = extra->data.F32[A]; \
+    extra->data.F32[A] = extra->data.F32[B]; \
+    extra->data.F32[B] = tmp; \
+  } \
+}
+
+bool psphotPetrosianSortPair (psVector *index, psVector *extra) {
+
+    // sort the vector set by the radius
+    PSSORT (index->n, COMPARE_INDEX, SWAP_INDEX, NONE);
+    return true;
+}
Index: branches/eam_branches/20090715/psphot/src/pmPetrosian.h
===================================================================
--- branches/eam_branches/20090715/psphot/src/pmPetrosian.h	(revision 25032)
+++ branches/eam_branches/20090715/psphot/src/pmPetrosian.h	(revision 25105)
@@ -12,13 +12,51 @@
 
 typedef struct {
-    psArray *radii;			// radii for raw radial profiles at evenly-spaced angles
-    psArray *fluxes;			// fluxes measured at above radii
+    psArray  *radii;			// radii for raw radial profiles at evenly-spaced angles
+    psArray  *fluxes;			// fluxes measured at above radii
     psVector *theta;			// angles corresponding to above radial profiles
-    psVector *isophotalRadii;
-    psEllipseAxes axes;
+    psVector *isophotalRadii;		// isophotal radius for the above angles
+
+    psVector *radiusElliptical;		// normalized radial coordinates for all relevant pixels
+    psVector *fluxElliptical;		// flux for the above radial coordinates
+
+    psVector *binnedFlux;		// mean surface brightness within radial bins
+    psVector *radialBins;		// radii corresponding to above binnedBlux
+    psVector *area;			// differential area of the non-overlapping radial bins
+
+    psEllipseAxes axes;			// shape of elliptical contour
+
+    float petrosianRadius;
+    float petrosianFlux;
+
 } pmPetrosian;
 
 pmPetrosian *pmPetrosianAlloc();
+bool psphotPetrosianFreeVectors(pmPetrosian *petrosian);
+
+bool psphotPetrosianProfile (pmSource *source);
+bool psphotRadialProfilesByAngles (pmPetrosian *petro, pmSource *source, int Nsec, float Rmax);
+float psphotRadiusFromProfile (psVector *radius, psVector *flux, float fluxMin, float fluxMax);
+bool psphotRadiiFromProfiles (pmPetrosian *petrosian, float fluxMin, float fluxMax);
+bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian);
+bool psphotEllipticalContour (pmPetrosian *petrosian);
+bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax);
+bool psphotPetrosianStats (pmPetrosian *petrosian);
+
+bool psphotPetrosianSortPair (psVector *index, psVector *extra);
+
+
+bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux);
+bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float RadiusRef);
+bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian);
+
+bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, 
+				 psVector *refRadius, psVector *meanSB, 
+				 psVector *petRatio, psVector *fluxSum, 
+				 float petRadius, float petFlux);
+
+bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec);
 
 /// @}
+
+
 # endif /* PM_PETROSIAN_H */
Index: branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c
===================================================================
--- branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c	(revision 25032)
+++ branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c	(revision 25105)
@@ -5,4 +5,5 @@
 // model parameters
 enum {PAR_PHI, PAR_EPSILON, PAR_RMIN};
+psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord);
 
 bool psphotEllipticalContour (pmPetrosian *petrosian) {
@@ -10,5 +11,5 @@
     // use LMM to fit theta vs radius to an ellipse
     psVector *theta = petrosian->theta;
-    psVector *radius = petrosian->isophotalRadius;
+    psVector *radius = petrosian->isophotalRadii;
 
     // find Rmin and Rmax for the initial guess
@@ -20,5 +21,7 @@
     psArray *x = psArrayAlloc(2*radius->n);
     psVector *y = psVectorAlloc(2*radius->n, PS_TYPE_F32);
+    psVector *yErr = psVectorAlloc(2*radius->n, PS_TYPE_F32);
 
+    int n = 0;
     for (int i = 0; i < radius->n; i++) {
 
@@ -28,7 +31,8 @@
 	coord = psVectorAlloc (2, PS_TYPE_F32);
 	coord->data.F32[1] = 0.0;
-	coord->data.F32[0] = radius->data.F32[i]*cos(theta->data.F32[i]);
+	coord->data.F32[0] = theta->data.F32[i];
 	x->data[n] = coord;
-	y->data.F32[n] = theta->data.F32[i];
+	y->data.F32[n] = radius->data.F32[i]*cos(theta->data.F32[i]);
+	yErr->data.F32[n] = 1000.0;
 	n++;
 
@@ -36,7 +40,8 @@
 	coord = psVectorAlloc (2, PS_TYPE_F32);
 	coord->data.F32[1] = 1.0;
-	coord->data.F32[0] = radius->data.F32[i]*sin(theta->data.F32[i]);
+	coord->data.F32[0] = theta->data.F32[i];
 	x->data[n] = coord;
-	y->data.F32[n] = theta->data.F32[i];
+	y->data.F32[n] = radius->data.F32[i]*sin(theta->data.F32[i]);
+	yErr->data.F32[n] = 1000.0;
 	n++;
 
@@ -49,4 +54,6 @@
 
     psVector *params = psVectorAlloc (3, PS_TYPE_F32);
+    
+    // psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 7);
     
     // create the minimization constraints
@@ -64,8 +71,8 @@
     
     // XXX skip the weights for now
-    psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psphotEllipticalContourFunc);
+    psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, psphotEllipticalContourFunc);
 
     fprintf (stderr, "# fitted values:\n");
-    fprintf (stderr, "Po:  %f\n", params->data.F32[PAR_PHI]);
+    fprintf (stderr, "Po:  %f\n", params->data.F32[PAR_PHI]*PS_DEG_RAD);
     fprintf (stderr, "Ep:  %f\n", params->data.F32[PAR_EPSILON]);
     fprintf (stderr, "Rm:  %f\n", params->data.F32[PAR_RMIN]);
@@ -75,4 +82,11 @@
     petrosian->axes.minor = params->data.F32[PAR_RMIN];
     petrosian->axes.theta = params->data.F32[PAR_PHI];
+
+    // show the results
+    psphotPetrosianVisualEllipticalContour (petrosian);
+
+    psFree (x);
+    psFree (y);
+    psFree (yErr);
 
     return true;
@@ -88,10 +102,13 @@
 psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord) {
 
+    static int pass = 0;
+
     psF32 *par = params->data.F32;
+
+    float alpha = coord->data.F32[0];
 
     float cs_alpha = cos(alpha);
     float sn_alpha = sin(alpha);
 
-    float alpha = coord->data.F32[0];
     float cs_phi = cos(alpha - par[PAR_PHI]);
     float sn_phi = sin(alpha - par[PAR_PHI]);
@@ -103,5 +120,7 @@
 
     // value is X
-    if (coord->data.F32[1] == 0) {
+    // if (coord->data.F32[1] == 0) {
+    if (pass == 0) {
+	pass = 1;
 
 	float value = par[PAR_RMIN]*cs_alpha*r;
@@ -109,7 +128,7 @@
 	if (deriv) {
 	    psF32 *dPAR = deriv->data.F32;
-	    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;
+	    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);
@@ -117,5 +136,7 @@
 
     // value is Y
-    if (coord->data.F32[1] == 1) {
+    // if (coord->data.F32[1] == 1) {
+    if (pass == 1) {
+	pass = 0;
 
 	float value = par[PAR_RMIN]*sn_alpha*r;
@@ -123,7 +144,7 @@
 	if (deriv) {
 	    psF32 *dPAR = deriv->data.F32;
-	    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;
+	    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);
Index: branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c
===================================================================
--- branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c	(revision 25032)
+++ branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c	(revision 25105)
@@ -3,9 +3,9 @@
 bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian) {
 
-    petrosian->radiusElli = psArrayAllocEmpty(100);
-    petrosian->fluxElli = psArrayAllocEmpty(100);
+    petrosian->radiusElliptical = psVectorAllocEmpty(100, PS_TYPE_F32);
+    petrosian->fluxElliptical = psVectorAllocEmpty(100, PS_TYPE_F32);
 
-    psVector *radius = petrosian->radiusElli;
-    psVector *flux = petrosian->fluxElli;
+    psVector *radius = petrosian->radiusElliptical;
+    psVector *flux = petrosian->fluxElliptical;
 
     // the psEllipse functions use z = 0.5(x/Sxx)^2 + 0.5(y/Syy)^2 + x y Sxy
@@ -14,9 +14,17 @@
     // a = A / sqrt(2)
 
+    // we have the shape parameters of the elliptical contour at the reference isophote.
+    // use the axis ratio (major/minor) to rescale the radial profile so that 1 pixel
+    // along the major axis is 1 pixel, and a smaller amount on the minor axis
+
     psEllipseAxes axes;
-    axes.major = petrosian->axes.major * M_SQRT1_2;
-    axes.minor = petrosian->axes.minor * M_SQRT1_2;
+    axes.major = M_SQRT1_2;
+    axes.minor = M_SQRT1_2 * (petrosian->axes.minor / petrosian->axes.major);
+
+    // axes.major = 1.0;
+    // axes.minor = petrosian->axes.minor / petrosian->axes.major;
+
     axes.theta = petrosian->axes.theta;
-    psEllipseShape shape = psEllipseShapeFromAxes (petrosian->axes);
+    psEllipseShape shape = psEllipseAxesToShape (axes);
 
     float Sxx = shape.sx;
@@ -24,15 +32,38 @@
     float Syy = shape.sy;
 
-    for (int iy = 0; iy < Ny; iy++) {
-	for (int ix = 0; ix < Nx; ix++) {
+    psVector *radiusRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *fluxRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+    for (int iy = 0; iy < source->pixels->numRows; iy++) {
+	for (int ix = 0; ix < source->pixels->numCols; ix++) {
 
 	    float x = ix - source->peak->xf + source->pixels->col0;
 	    float y = iy - source->peak->yf + source->pixels->row0;
 
-	    psVectorAppend(radius, sqrt(0.5*PS_SQR(x/Sxx) + 0.5*PS_SQR(y/Syy) + x*y*Sxy));
-	    psVectorAppend(flux, pixels[iy][ix]);
+	    float r2 = 0.5*PS_SQR(x/Sxx) + 0.5*PS_SQR(y/Syy) + x*y*Sxy;
+	    // float r2 = PS_SQR(x/Sxx) + PS_SQR(y/Syy) + x*y*Sxy;
+
+	    psVectorAppend(radius, sqrt(r2));
+	    psVectorAppend(flux, source->pixels->data.F32[iy][ix]);
+
+	    float Rraw = hypot(x, y);
+	    psVectorAppend(radiusRaw, Rraw);
+	    psVectorAppend(fluxRaw, source->pixels->data.F32[iy][ix]);
 	}
     }
 
+    // psVector *radiusRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
+    // psVector *fluxRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
+    // for (int i = 0; i < petrosian->radii->n; i++) {
+    //   psVector *r = petrosian->radii->data[i];
+    //   psVector *f = petrosian->fluxes->data[i];
+    //   for (int j = 0; j < r->n; j++) {
+    // 	psVectorAppend(radiusRaw, r->data.F32[j]);
+    // 	psVectorAppend(fluxRaw, f->data.F32[j]);
+    //   }
+    // }
+
+    // psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, 0.0);
+    psphotPetrosianVisualProfileByAngle (radius, flux);
     return true;
 }
Index: branches/eam_branches/20090715/psphot/src/psphotInternal.h
===================================================================
--- branches/eam_branches/20090715/psphot/src/psphotInternal.h	(revision 25032)
+++ branches/eam_branches/20090715/psphot/src/psphotInternal.h	(revision 25105)
@@ -14,4 +14,5 @@
 #include <psmodules.h>
 #include "psphot.h"
+#include "pmPetrosian.h"
 
 #endif
Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c
===================================================================
--- branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c	(revision 25032)
+++ branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c	(revision 25105)
@@ -1,13 +1,31 @@
 # include "psphotInternal.h"
+
+// generate the Petrosian radius and flux using elliptical contours
+
+// XXX much of this function is focused on generating the clean contours, which can be used by 
+// any number of aperture-like measurements.  probably will want to rename the pmPetrosian
+// structure to something the pmRadialProfile
 
 bool psphotPetrosianProfile (pmSource *source) {
 
+  // container to hold results from the radial profile analysis
   pmPetrosian *petrosian = pmPetrosianAlloc();
 
-  if (!psphotRadialProfilesByAngle (petrosian, source, Nsec, Rmax)) {
+  int Nsec = 24;
+  float Rmax = 64;
+  float fluxMin = 0.0;
+  float fluxMax = 1000.0;
+
+  // generate a series of radial profiles at Nsec evenly spaced angles.  the profile flux
+  // is measured by interpolation for small radii; for large radii, the pixels in a box
+  // are averaged to increase the S/N (XXX not yet done)
+  if (!psphotRadialProfilesByAngles (petrosian, source, Nsec, Rmax)) {
     psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian");
     return false;
   }
 
+  // use the radial profiles to determine the radius of a given isophote.  this isophote
+  // is used to determine the elliptical shape of the object, so it has a relatively high
+  // value (nominally 50% of the peak)
   if (!psphotRadiiFromProfiles (petrosian, fluxMin, fluxMax)) {
     psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles");
@@ -15,4 +33,5 @@
   }
 
+  // convert the isophotal radius vs angle measurements to an elliptical contour
   if (!psphotEllipticalContour (petrosian)) {
     psError (PS_ERR_UNKNOWN, false, "failed to measure elliptical contour");
@@ -20,4 +39,29 @@
   }
   
+  // generate a single, normalized radial profile following the elliptical contours.
+  // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel
+  if (!psphotEllipticalProfile (source, petrosian)) {
+    psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
+    return false;
+  }
+  
+  // integrate the radial profile for radial bins defined for the petrosian measurement:
+  // SB_i (r_i) where \alpha r_i < r < \beta r_i
+  if (!psphotPetrosianRadialBins (source, petrosian, Rmax)) {
+    psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
+    return false;
+  }
+  
+  // use the SB_i from above to calculate the petrosian radius and the flux within that radius
+  if (!psphotPetrosianStats (petrosian)) {
+    psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
+    return false;
+  }
+  
+  psphotPetrosianFreeVectors(petrosian);
 
+  fprintf (stderr, "petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f\n",
+	   petrosian->petrosianRadius, petrosian->petrosianFlux, petrosian->axes.minor/petrosian->axes.major, PS_DEG_RAD*petrosian->axes.theta);
+
+  return true;
 }
Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c
===================================================================
--- branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c	(revision 25032)
+++ branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c	(revision 25105)
@@ -2,34 +2,67 @@
 
 // convert the flux vs elliptical radius to annular bins
-bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian) {
 
-    psVector *radius = petrosian->radiusElli;
-    psVector *flux = petrosian->fluxElli;
+// we are guaranteed to be limited by either the seeing (1 - few pixels) or by the pixels
+// themselves.  this function does not attempt to measure the radial profiles accurately
+// for radii that are smaller than a minimum (currently 1.0 pixels).  
 
-    sort (radius, flux);
+// for small radii, we are measuring the mean surface brightness in non-overlapping radial
+// bins.  for large radii (r > 2 pixels), we are measuring the surface brightness for a
+// radius range \alpha r_i < i < \beta r_i, but performing this measurement for radii more
+// finely spaced than r_{i+1} = r_i * \beta / \alpha.  for the integration, we need to
+// track the non-overlapping radius values.
+
+# define PETROSIAN_ALPHA 0.8
+# define PETROSIAN_BETA 1.25
+
+bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax) {
+
+    psVector *radius = petrosian->radiusElliptical;
+    psVector *flux = petrosian->fluxElliptical;
+
+    // sort incoming vectors by radius
+    psphotPetrosianSortPair (radius, flux);
+
+    int nMax = radiusMax;
 
     // radBin stores the centers of the radial bins, 
     // radMin, radMax store the bounds
-    psVector *radBin = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *radMin = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *radMax = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *radMin  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+    psVector *radMax  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+    psVector *radAlp  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+    psVector *radBet  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+
+    psVector *fluxBin = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+    psVector *radBin  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+    psVector *area    = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+
+    psVectorInit (fluxBin, 0.0);
+    psVectorInit (radBin, 0.0);
 
     // generate radial bin bounds
     radMin->data.F32[0] = 0.0;
     radMax->data.F32[0] = 1.0;
+    radAlp->data.F32[0] = 0.0;
+    radBet->data.F32[0] = 1.0;
     
     radMin->data.F32[1] = 1.0;
     radMax->data.F32[1] = 1.5;
+    radAlp->data.F32[1] = 1.0;
+    radBet->data.F32[1] = 1.5;
     
     radMin->data.F32[2] = 1.5;
     radMax->data.F32[2] = 2.0;
+    radAlp->data.F32[2] = 1.5;
+    radBet->data.F32[2] = 2.0;
     
-    for (int i = 3; i < rmax; i++) {
-	radMin->data.F32[i] = (i - 1);
-	radMax->data.F32[i] = i;
+    int nPts = 3;
+    for (int i = 3; i < radiusMax; i++) {
+	radMin->data.F32[nPts] = (i - 1);
+	radMax->data.F32[nPts] = i;
+	nPts++;
     }
+    radMin->n = radMax->n = radAlp->n = radBet->n = nPts;
 
     // generate radial area-weighted mean radius & non-overlapping areas
-    float Rprev = 0.0;
     for (int i = 0; i < radMin->n; i++) {
 	float rMin = radMin->data.F32[i];
@@ -42,10 +75,69 @@
 	float rMax3 = rMax2*rMax;
 
-	float rBin = (rMax3 - rMin3) / (rMax2 - rMin2) / 3.0;
+	float rBin = 2.0 * (rMax3 - rMin3) / (rMax2 - rMin2) / 3.0;
 	
+	// XXX calculate area-weighted radius rather than asserting?
 	radBin->data.F32[i] = rBin;
+	area->data.F32[i] = M_PI * (rMax2 - rMin2);
 
-	area->data.F32[i] = 
+	if (i > 2) {
+	    radAlp->data.F32[i] = rBin*PETROSIAN_ALPHA;
+	    radBet->data.F32[i] = rBin*PETROSIAN_BETA;
+	}
     }
+
+    // storage vector for stats
+    psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+
+    // integrate flux, radius for each of these bins.  since flux is sorted by radius, 
+    // we can do this fairly quickly
+
+    bool done = false;
+    int nOut = 0;
+    float Rmin = radAlp->data.F32[nOut];
+    float Rmax = radBet->data.F32[nOut];
+    float Rnxt = radAlp->data.F32[nOut+1];  // minimum radius for next range
+    int iNext = 0;
+    for (int i = 0; !done && (i < radius->n); i++) {
+	if (radius->data.F32[i] < Rnxt) {
+	  iNext = i;
+	}
+	if (radius->data.F32[i] > Rmax) {
+	    // calculate the value for the nOut bin
+	    psVectorStats (stats, values, NULL, NULL, 0);
+	    fluxBin->data.F32[nOut] = stats->sampleMedian;
+	    nOut ++;
+	    if (nOut >= nMax) break;
+	    Rmin = radAlp->data.F32[nOut];
+	    Rmax = radBet->data.F32[nOut];
+	    Rnxt = (nOut < nMax - 1) ? radAlp->data.F32[nOut+1] : Rmax;  // minimum radius for next range
+	    values->n = 0;
+	    psStatsInit(stats);
+	    i = iNext;
+	}
+	if (radius->data.F32[i] < Rmin) {
+	    continue;
+	}
+	psVectorAppend (values, flux->data.F32[i]);
+    }
+    fluxBin->n = radBin->n = area->n = nOut;
+    // XXX I think this misses the last radial bin -- do we care?
+
+    // save the vectors
+    petrosian->radialBins = radBin;
+    petrosian->area = area;
+    petrosian->binnedFlux = fluxBin;
+
+    psphotPetrosianVisualProfileRadii (radius, flux, radBin, fluxBin, 0.0);
+
+    psFree(radMin);
+    psFree(radMax);
+    psFree(radAlp);
+    psFree(radBet);
+    psFree(values);
+    psFree(stats);
+
+    return true;
 }
 
Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c
===================================================================
--- branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c	(revision 25105)
+++ branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c	(revision 25105)
@@ -0,0 +1,79 @@
+# include "psphotInternal.h"
+
+# define PETROSIAN_RATIO 0.05
+
+// generate the Petrosian radius and flux from the mean surface brightness (r_i)
+
+bool psphotPetrosianStats (pmPetrosian *petrosian) {
+
+  float petRadius, petFlux;
+
+  psVector *fluxBin = petrosian->binnedFlux;
+  psVector *radBin  = petrosian->radialBins;
+  psVector *area    = petrosian->area;
+
+  psVector *fluxSum  = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32);
+  psVector *refRadius= psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32);
+  psVector *petRatio = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32);
+  psVector *meanSB   = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32);
+  psVector *areaSum  = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32);
+
+  bool above = true;
+  float Fsum = 0.0;
+  float Asum = 0.0;
+  float Area = area->data.F32[0];
+
+  int nOut = 0;
+  for (int i = 0; i < fluxBin->n; i++) {
+    // for nan bins, we keep the area for use with the next valid bin
+    if (!isfinite(fluxBin->data.F32[i])) {
+      Area += area->data.F32[i];
+      continue;
+    } 
+    Fsum += fluxBin->data.F32[i] * Area;
+    Asum += Area;
+    if (i+1 < fluxBin->n) {
+      Area = area->data.F32[i+1];
+    }
+
+    psVectorAppend(petRatio, Asum * fluxBin->data.F32[i] / Fsum);
+    psVectorAppend(fluxSum, Fsum);
+    psVectorAppend(meanSB, Fsum / Asum);
+    psVectorAppend(areaSum, Asum);
+    psVectorAppend(refRadius, radBin->data.F32[i]);
+
+    // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
+    // we will keep and report the last (largest radius) value
+    if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO)) {
+      // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
+      if (i == 0) { 
+	// assume Fmax @ R = 0.0
+	petRadius = radBin->data.F32[nOut]  * (PETROSIAN_RATIO - 1.0) / (petRatio->data.F32[nOut] - 1.0);
+	petFlux   = fluxSum->data.F32[nOut] * (PETROSIAN_RATIO - 1.0) / (petRatio->data.F32[nOut] - 1.0);
+      } else {
+	petRadius = radBin->data.F32[nOut-1]  + (radBin->data.F32[nOut]  - radBin->data.F32[nOut-1])  * (PETROSIAN_RATIO - petRatio->data.F32[nOut-1]) / (petRatio->data.F32[nOut] - petRatio->data.F32[nOut-1]);
+	petFlux   = fluxSum->data.F32[nOut-1] + (fluxSum->data.F32[nOut] - fluxSum->data.F32[nOut-1]) * (PETROSIAN_RATIO - petRatio->data.F32[nOut-1]) / (petRatio->data.F32[nOut] - petRatio->data.F32[nOut-1]);
+      }
+      above = false;
+    }
+    
+    // reset on transitions up, but do not re-calculate rad_90, flux_90
+    if (!above && (petRatio->data.F32[nOut] >= PETROSIAN_RATIO)) {
+      above = true;
+    }
+    nOut ++;
+  }
+
+  // save petRadius, petFlux
+  petrosian->petrosianRadius = petRadius;
+  petrosian->petrosianFlux   = petFlux;
+
+  psphotPetrosianVisualStats (radBin, fluxBin, refRadius, meanSB, petRatio, fluxSum, petRadius, petFlux);
+
+  psFree(fluxSum);
+  psFree(petRatio);
+  psFree(meanSB);
+  psFree(areaSum);
+
+  return true;
+}
Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c
===================================================================
--- branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c	(revision 25105)
+++ branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c	(revision 25105)
@@ -0,0 +1,111 @@
+# include "psphotInternal.h"
+
+int main (int argc, char **argv) {
+
+  pmErrorRegister();                  // register psModule's error codes/messages
+  pmModelClassInit();
+
+  int N;
+
+  // XXX add noise and seeing.
+  // XXX double check on sersic functional form
+  // XXX modify ratio if ratio > 1.0 (swap major and minor)
+
+  float peak = 1000.0;
+  float sigma = 2.0;	      // major axis size 
+  float ARatio = 1.0;
+  float angle = 0.0;
+  float sersic = 0.5;
+  float skynoise = 0.0;
+
+  if ((N = psArgumentGet (argc, argv, "-peak"))) {
+    psArgumentRemove (N, &argc, argv);
+    peak = atof(argv[N]);
+    psArgumentRemove (N, &argc, argv);
+  }
+  if ((N = psArgumentGet (argc, argv, "-sigma"))) {
+    psArgumentRemove (N, &argc, argv);
+    sigma = atof(argv[N]);
+    psArgumentRemove (N, &argc, argv);
+  }
+  if ((N = psArgumentGet (argc, argv, "-aratio"))) {
+    psArgumentRemove (N, &argc, argv);
+    ARatio = atof(argv[N]);
+    psArgumentRemove (N, &argc, argv);
+  }
+  if ((N = psArgumentGet (argc, argv, "-angle"))) {
+    psArgumentRemove (N, &argc, argv);
+    angle = PS_RAD_DEG*atof(argv[N]);
+    psArgumentRemove (N, &argc, argv);
+  }
+  if ((N = psArgumentGet (argc, argv, "-sersic"))) {
+    psArgumentRemove (N, &argc, argv);
+    sersic = atof(argv[N]);
+    psArgumentRemove (N, &argc, argv);
+  }
+  if ((N = psArgumentGet (argc, argv, "-skynoise"))) {
+    psArgumentRemove (N, &argc, argv);
+    skynoise = atof(argv[N]);
+    psArgumentRemove (N, &argc, argv);
+  }
+  if ((N = psArgumentGet (argc, argv, "-visual"))) {
+    psArgumentRemove (N, &argc, argv);
+    pmVisualSetVisual(true);
+  }
+
+  if (argc != 1) {
+    fprintf (stderr, "USAGE: psphotPetrosianStudy\n");
+    exit (2);
+  }
+
+  // create a containing image & associated readout
+  pmReadout *readout = pmReadoutAlloc(NULL);
+  readout->image = psImageAlloc(128, 128, PS_TYPE_F32);
+
+  // create a dummy variance, but don't populate -- it is not used by pmSourceMoments if the sigma parameters is set to 0.0
+  readout->variance = psImageAlloc(128, 128, PS_TYPE_F32);
+
+  // create a model & associated source
+  pmModelType type = pmModelClassGetType("PS_MODEL_SERSIC");
+  pmModel *model = pmModelAlloc(type);
+
+  // set the model parameters
+  model->params->data.F32[PM_PAR_SKY]  = 0.0;
+  model->params->data.F32[PM_PAR_I0]   = peak;
+  model->params->data.F32[PM_PAR_XPOS] = 64.0;
+  model->params->data.F32[PM_PAR_YPOS] = 64.0;
+
+  psEllipseAxes axes;
+  axes.major = sigma;
+  axes.minor = sigma*ARatio;
+  axes.theta = angle;
+
+  psEllipseShape shape = psEllipseAxesToShape (axes);
+
+  // XXX set the sigma with user input
+  model->params->data.F32[PM_PAR_SXX]  = shape.sx * M_SQRT2;
+  model->params->data.F32[PM_PAR_SYY]  = shape.sy * M_SQRT2;
+  model->params->data.F32[PM_PAR_SXY]  = shape.sxy;
+
+  if (model->params->n > 7) {
+    model->params->data.F32[PM_PAR_7]  = sersic;
+  }
+
+  pmSource *source = pmSourceFromModel(model, readout, 32.0, PM_SOURCE_TYPE_STAR);
+
+  // generate the modelFlux 
+  pmSourceCacheModel(source, 0);
+
+  // instantiate the source
+  pmSourceAdd(source, PM_MODEL_OP_FUNC, 0); 
+
+  // XXX add noise here...
+
+  psphotSaveImage(NULL, readout->image, "sersic.fits");
+
+  psphotPetrosianProfile (source);
+
+  psFree (source);
+
+  exit (0);
+}
Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c
===================================================================
--- branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c	(revision 25105)
+++ branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c	(revision 25105)
@@ -0,0 +1,393 @@
+# include "psphotInternal.h"
+
+// this function displays representative images as the psphot analysis progresses:
+// 0 : image, 1 : variance
+// 0 : backsub, 1 : variance, 2 : backgnd
+// 0 : backsub, 1 : variance, 2 : signif
+// (overlay peaks on images)
+// (overlay footprints on images)
+// (overlay moments on images)
+// (overlay rough class on images)
+// 0 : backsub, 1 : psfpos, 2: psfsub
+// 0 : backsub, 1 : lin_resid, 2: psfsub
+
+# if (HAVE_KAPA)
+# include <kapa.h>
+
+// functions used to visualize the analysis as it goes
+// these are invoked by the -visual options
+
+static int kapa = -1;
+
+// if no valid data is supplied (NULL or n <- 0), leave limits as they were
+bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec) {
+
+    if (xVec && xVec->n > 0) {
+	graphdata->xmin = graphdata->xmax = xVec->data.F32[0];
+	for (int i = 1; i < xVec->n; i++) {
+	    graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]);
+	    graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]);
+	}
+	float range = graphdata->xmax - graphdata->xmin;
+	graphdata->xmax += 0.05*range;
+	graphdata->xmin -= 0.05*range;
+    }
+    if (yVec && yVec->n > 0) {
+	graphdata->ymin = graphdata->ymax = yVec->data.F32[0];
+	for (int i = 1; i < yVec->n; i++) {
+	    graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]);
+	    graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]);
+	}
+	float range = graphdata->ymax - graphdata->ymin;
+	graphdata->ymax += 0.05*range;
+	graphdata->ymin -= 0.05*range;
+    }
+    return true;
+}
+
+bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux) {
+
+    Graphdata graphdata;
+
+    // return true;
+    if (!pmVisualIsVisual()) return true;
+
+    if (kapa == -1) {
+        kapa = KapaOpenNamedSocket ("kapa", "psphot:plots");
+        if (kapa == -1) {
+            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+            pmVisualSetVisual(false);
+            return false;
+        }
+    }
+
+    KapaClearPlots (kapa);
+    KapaInitGraph (&graphdata);
+    KapaSetFont (kapa, "courier", 14);
+
+    pmVisualLimitsFromVectors (&graphdata, radius, flux);
+    KapaSetLimits (kapa, &graphdata);
+
+    KapaBox (kapa, &graphdata);
+    KapaSendLabel (kapa, "radius", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "flux", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.style = 2;
+    graphdata.ptype = 0;
+    graphdata.size = 1.0;
+    KapaPrepPlot (kapa, radius->n, &graphdata);
+    KapaPlotVector (kapa, radius->n, radius->data.F32, "x");
+    KapaPlotVector (kapa, radius->n, flux->data.F32, "y");
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+        psWarning("Unable to read option");
+    }
+    return true;
+}
+
+bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float RadiusRef) {
+
+    float FluxRef = 500.0;
+
+    Graphdata graphdata;
+
+    // return true;
+    if (!pmVisualIsVisual()) return true;
+
+    if (kapa == -1) {
+        kapa = KapaOpenNamedSocket ("kapa", "psphot:plots");
+        if (kapa == -1) {
+            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+            pmVisualSetVisual(false);
+            return false;
+        }
+    }
+
+    KapaClearPlots (kapa);
+    KapaInitGraph (&graphdata);
+    KapaSetFont (kapa, "courier", 14);
+
+    pmVisualLimitsFromVectors (&graphdata, radius, flux);
+    KapaSetLimits (kapa, &graphdata);
+
+    KapaBox (kapa, &graphdata);
+    KapaSendLabel (kapa, "radius", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "flux", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.style = 2;
+    graphdata.ptype = 0;
+    graphdata.size = 1.0;
+    KapaPrepPlot (kapa, radius->n, &graphdata);
+    KapaPlotVector (kapa, radius->n, radius->data.F32, "x");
+    KapaPlotVector (kapa, radius->n, flux->data.F32, "y");
+
+    // do this with log-r, log-flux?
+    graphdata.color = KapaColorByName ("blue");
+    graphdata.style = 2;
+    graphdata.ptype = 2;
+    graphdata.size = 2.0;
+    KapaPrepPlot   (kapa, radiusBin->n, &graphdata);
+    KapaPlotVector (kapa, radiusBin->n, radiusBin->data.F32, "x");
+    KapaPlotVector (kapa, radiusBin->n, fluxBin->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.style = 2;
+    graphdata.ptype = 7;
+    graphdata.size = 3.0;
+    KapaPrepPlot (kapa, 1, &graphdata);
+    KapaPlotVector (kapa, 1, &RadiusRef, "x");
+    KapaPlotVector (kapa, 1, &FluxRef, "y");
+
+    fprintf (stderr, "radius: %f\n", RadiusRef);
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+        psWarning("Unable to read option");
+    }
+    return true;
+}
+
+bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, 
+				 psVector *refRadius, psVector *meanSB, 
+				 psVector *petRatio, psVector *fluxSum, 
+				 float petRadius, float petFlux)
+{
+    Graphdata graphdata;
+    KapaSection section;
+
+    if (!pmVisualIsVisual()) return true;
+
+    if (kapa == -1) {
+        kapa = KapaOpenNamedSocket ("kapa", "psphot:plots");
+        if (kapa == -1) {
+            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+            pmVisualSetVisual(false);
+            return false;
+        }
+    }
+
+    KapaClearPlots (kapa);
+    KapaInitGraph (&graphdata);
+    KapaSetFont (kapa, "courier", 14);
+
+    // radius vs flux
+    // radius vs mean SB
+    // radius vs petRatio
+
+    // *** section 1: radius vs mean SB
+    section.dx = 1.00;
+    section.dy = 0.33;
+    section.x  = 0.00;
+    section.y  = 0.00;
+    section.name = psStringCopy ("meanSB");
+    KapaSetSection (kapa, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    pmVisualLimitsFromVectors (&graphdata, radBin, fluxBin);
+    KapaSetLimits (kapa, &graphdata);
+
+    KapaBox (kapa, &graphdata);
+    KapaSendLabel (kapa, "radius", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "mean SB", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.style = 2;
+    graphdata.ptype = 0;
+    graphdata.size = 1.0;
+    KapaPrepPlot (kapa, radBin->n, &graphdata);
+    KapaPlotVector (kapa, radBin->n, radBin->data.F32, "x");
+    KapaPlotVector (kapa, radBin->n, fluxBin->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.style = 2;
+    graphdata.ptype = 1;
+    graphdata.size = 2.0;
+    KapaPrepPlot (kapa, refRadius->n, &graphdata);
+    KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x");
+    KapaPlotVector (kapa, refRadius->n, meanSB->data.F32, "y");
+
+    // *** section 2: radius vs petrosian ratio
+    section.dx = 1.00;
+    section.dy = 0.33;
+    section.x  = 0.00;
+    section.y  = 0.33;
+    section.name = psStringCopy ("ratio");
+    KapaSetSection (kapa, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    pmVisualLimitsFromVectors (&graphdata, radBin, petRatio);
+    KapaSetLimits (kapa, &graphdata);
+
+    KapaBox (kapa, &graphdata);
+    KapaSendLabel (kapa, "radius", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "ratio", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.style = 2;
+    graphdata.ptype = 0;
+    graphdata.size = 1.0;
+    KapaPrepPlot (kapa, refRadius->n, &graphdata);
+    KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x");
+    KapaPlotVector (kapa, refRadius->n, petRatio->data.F32, "y");
+
+    float ratio_90 = 0.1;
+    graphdata.color = KapaColorByName ("red");
+    graphdata.style = 2;
+    graphdata.ptype = 2;
+    graphdata.size = 2.0;
+    KapaPrepPlot   (kapa, 1, &graphdata);
+    KapaPlotVector (kapa, 1, &petRadius, "x");
+    KapaPlotVector (kapa, 1, &ratio_90, "y");
+
+    // *** section 3: radius vs integrated flux
+    section.dx = 1.00;
+    section.dy = 0.33;
+    section.x  = 0.00;
+    section.y  = 0.66;
+    section.name = psStringCopy ("flux");
+    KapaSetSection (kapa, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    pmVisualLimitsFromVectors (&graphdata, radBin, fluxSum);
+    KapaSetLimits (kapa, &graphdata);
+
+    KapaBox (kapa, &graphdata);
+    KapaSendLabel (kapa, "radius", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "integrated flux", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.style = 2;
+    graphdata.ptype = 0;
+    graphdata.size = 1.0;
+    KapaPrepPlot   (kapa, refRadius->n, &graphdata);
+    KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x");
+    KapaPlotVector (kapa, refRadius->n, fluxSum->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 2;
+    graphdata.style = 2;
+    graphdata.size = 2.0;
+    KapaPrepPlot   (kapa, 1, &graphdata);
+    KapaPlotVector (kapa, 1, &petRadius, "x");
+    KapaPlotVector (kapa, 1, &petFlux, "y");
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+        psWarning("Unable to read option");
+    }
+    return true;
+}
+
+bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian) {
+
+    Graphdata graphdata;
+
+    if (!pmVisualIsVisual()) return true;
+
+    if (kapa == -1) {
+        kapa = KapaOpenNamedSocket ("kapa", "psphot:plots");
+        if (kapa == -1) {
+            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+            pmVisualSetVisual(false);
+            return false;
+        }
+    }
+
+    KapaClearPlots (kapa);
+    KapaInitGraph (&graphdata);
+    KapaSetFont (kapa, "courier", 14);
+
+    psVector *theta = petrosian->theta;
+    psVector *radius = petrosian->isophotalRadii;
+
+    // find Rmin and Rmax for the initial guess
+    float Rmin = radius->data.F32[0];
+    float Rmax = radius->data.F32[0];
+
+    psVector *Rx = psVectorAlloc(radius->n, PS_TYPE_F32);
+    psVector *Ry = psVectorAlloc(radius->n, PS_TYPE_F32);
+
+    for (int i = 0; i < theta->n; i++) {
+	Rx->data.F32[i] = radius->data.F32[i]*cos(theta->data.F32[i]);
+	Ry->data.F32[i] = radius->data.F32[i]*sin(theta->data.F32[i]);
+
+	// check the radius range
+	Rmin = MIN (Rmin, radius->data.F32[i]);
+	Rmax = MAX (Rmax, radius->data.F32[i]);
+    }	
+
+    psVector *rx = psVectorAlloc(361, PS_TYPE_F32);
+    psVector *ry = psVectorAlloc(361, PS_TYPE_F32);
+
+    float epsilon = petrosian->axes.minor / petrosian->axes.major;
+
+    for (int i = 0; i < 361; i++) {
+
+	float alpha = PS_RAD_DEG * i;
+
+	float cs_alpha = cos(alpha);
+	float sn_alpha = sin(alpha);
+
+	float cs_phi = cos(alpha - petrosian->axes.theta);
+	float sn_phi = sin(alpha - petrosian->axes.theta);
+
+	float r = 1.0 / sqrt(SQ(sn_phi) + SQ(epsilon*cs_phi));
+
+	// generate the model fit here
+	rx->data.F32[i] = petrosian->axes.minor * cs_alpha * r;
+	ry->data.F32[i] = petrosian->axes.minor * sn_alpha * r;
+    }	
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.xmin = -1.1*Rmax;
+    graphdata.ymin = -1.1*Rmax;
+    graphdata.xmax = +1.1*Rmax;
+    graphdata.ymax = +1.1*Rmax;
+    KapaSetLimits (kapa, &graphdata);
+
+    KapaBox (kapa, &graphdata);
+    KapaSendLabel (kapa, "R_x", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "R_y", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.style = 2;
+    graphdata.ptype = 2;
+    graphdata.size = 1.0;
+    KapaPrepPlot (kapa, Rx->n, &graphdata);
+    KapaPlotVector (kapa, Rx->n, Rx->data.F32, "x");
+    KapaPlotVector (kapa, Rx->n, Ry->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.style = 0;
+    graphdata.ptype = 0;
+    graphdata.size = 1.0;
+    KapaPrepPlot (kapa, rx->n, &graphdata);
+    KapaPlotVector (kapa, rx->n, rx->data.F32, "x");
+    KapaPlotVector (kapa, rx->n, ry->data.F32, "y");
+
+    // pause and wait for user input:
+    // continue, save (provide name), ??
+    char key[10];
+    fprintf (stdout, "[c]ontinue? ");
+    if (!fgets(key, 8, stdin)) {
+        psWarning("Unable to read option");
+    }
+    return true;
+}
+
+# endif
Index: branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c
===================================================================
--- branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c	(revision 25032)
+++ branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c	(revision 25105)
@@ -8,15 +8,15 @@
 // XXX how does this fit in context?
 
-bool psphotRadialProfilesByAngles (pmPetrosianData *petro, pmSource *source, int Nsec, float Rmax) {
+bool psphotRadialProfilesByAngles (pmPetrosian *petrosian, pmSource *source, int Nsec, float Rmax) {
 
     float dtheta = 2.0*M_PI / Nsec;
 
-    psFree(petro->radii);
-    psFree(petro->fluxes);
-    psFree(petro->theta);
+    psFree(petrosian->radii);
+    psFree(petrosian->fluxes);
+    psFree(petrosian->theta);
 
-    petro->radii = psArrayAllocEmpty(Nsec);
-    petro->fluxes = psArrayAllocEmpty(Nsec);
-    petro->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32);
+    petrosian->radii = psArrayAllocEmpty(Nsec);
+    petrosian->fluxes = psArrayAllocEmpty(Nsec);
+    petrosian->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32);
 
 
@@ -25,6 +25,6 @@
 	float theta = i*dtheta;
 
-	psVector *radius = psVectorAllocEmpty(PS_TYPE_F32, Rmax);
-	psVector *flux   = psVectorAllocEmpty(PS_TYPE_F32, Rmax);
+	psVector *radius = psVectorAllocEmpty(Rmax, PS_TYPE_F32);
+	psVector *flux   = psVectorAllocEmpty(Rmax, PS_TYPE_F32);
 
 	// start at Xo,Yo and find the x,y locations for r_i, theta where r_i increments by 1 pixel
@@ -32,13 +32,18 @@
 	for (float r = 0; r < Rmax; r += 1.0) {
 
-	    float Xo = source->peak->xf - source->pixels->col0;
-	    float Yo = source->peak->yf - source->pixels->row0;
+	    float Xo = source->peak->xf + 0.5;
+	    float Yo = source->peak->yf + 0.5;
 
 	    // Xo,Yo are referenced to pixels with bounds i+0.0, i+1.0
-	    x = r * cos (theta) + Xo
-	    y = r * sin (theta) + Yo;
+	    float x = r * cos (theta) + Xo;
+	    float y = r * sin (theta) + Yo;
+
+	    if (x < 0) continue;
+	    if (y < 0) continue;
+	    if (x >= source->pixels->parent->numCols) continue;
+	    if (y >= source->pixels->parent->numRows) continue;
 
 	    // value is NAN if we run off the image
-	    value = psImageInterpolatePixelBilinear(x, y, image);
+	    float value = psImageInterpolatePixelBilinear(x, y, source->pixels);
 	    if (isnan(value)) continue;
 
@@ -47,7 +52,12 @@
 	}
 
-	psArrayAdd (petro->radii, 100, radius);
-	psArrayAdd (petro->fluxes, 100, flux);
-	psVectorAppend (petro->theta, theta);
+	psArrayAdd (petrosian->radii, 100, radius);
+	psArrayAdd (petrosian->fluxes, 100, flux);
+	psVectorAppend (petrosian->theta, theta);
+
+	// psphotPetrosianVisualProfileByAngle (radius, flux);
+
+	psFree(radius);
+	psFree(flux);
     }
 
Index: branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c
===================================================================
--- branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c	(revision 25032)
+++ branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c	(revision 25105)
@@ -10,5 +10,7 @@
       psVector *radii = petrosian->radii->data[i];
       psVector *fluxes = petrosian->fluxes->data[i];
-      float radius =  = psphotRadiusFromProfile (radii, fluxes, fluxMin, fluxMax);
+      float radius = psphotRadiusFromProfile (radii, fluxes, fluxMin, fluxMax);
+
+      // psphotPetrosianVisualProfileByAngle (radii, fluxes, radius);
 
       // warn on NAN?
@@ -44,4 +46,5 @@
 
     // rebin with bins 1/2 the size of the mean radius
+    int Rbin = 0;
     if (Rnpt == 0) {
 	Rbin = 1;
@@ -58,18 +61,19 @@
 	radiusBinned = psMemIncrRefCounter (radius);
     } else {
-	// storage vector for sorting
-	psVector *values = psVectorAllocEmpty (PS_TYPE_F32, flux->n);
+	// storage vector for stats
+	psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
+	psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
   
 	// rebinned vectors
-	fluxBinned = psVectorAllocEmpty (PS_TYPE_F32, flux->n);
-	radiusBinned = psVectorAllocEmpty (PS_TYPE_F32, flux->n);
+	fluxBinned = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
+	radiusBinned = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
 
 	// sort the flux by the radius
-	sort2vec (radius, flux);
+	psphotPetrosianSortPair (radius, flux);
 
 	int nOut = 0;
 	radiusBinned->data.F32[nOut] = (nOut + 0.5)*Rbin;
-	Rmin = radiusBinned->data.F32[nOut] - 0.5*Rbin;
-	Rmax = radiusBinned->data.F32[nOut] + 0.5*Rbin;
+	float Rmin = radiusBinned->data.F32[nOut] - 0.5*Rbin;
+	float Rmax = radiusBinned->data.F32[nOut] + 0.5*Rbin;
 
 	for (int i = 0; i < flux->n; i++) {
@@ -80,13 +84,13 @@
 	    if (radius->data.F32[i] > Rmax) {
 		// calculate the value for the nOut bin
-		psStats (stats, values);
-		fluxBinned->data.F32[nOut] = stats->value;
+		// XXX need to fix this as well psStats (stats, values);
+		psVectorStats (stats, values, NULL, NULL, 0);
+		fluxBinned->data.F32[nOut] = stats->sampleMedian;
 		nOut ++;
 		radiusBinned->data.F32[nOut] = (nOut + 0.5)*Rbin;
 		Rmin = radiusBinned->data.F32[nOut] - 0.5*Rbin;
 		Rmax = radiusBinned->data.F32[nOut] + 0.5*Rbin;
-		i--; 
-		continue;
-		// XXX the sequencing here is a bit ugly -- can we unfold this?
+		values->n = 0;
+		psStatsInit(stats);
 	    }
 	    psVectorAppend (values, flux->data.F32[i]);
@@ -95,8 +99,9 @@
 	radiusBinned->n = nOut;
 	psFree (values);
+	psFree(stats);
     }
 
-    Ro = NAN;
-    above = TRUE;
+    float Ro = NAN;
+    bool above = true;
     for (int i = 0; i < fluxBinned->n; i++) {
 
@@ -105,6 +110,5 @@
 	    // XXX is there a macro in psLib that does this interpolation?
 	    if (i == 0) { 
-		// assume Fmax @ R = 0.0
-		Ro = radiusBinned->data.F32[i] * (Fo - Fmax) / (fluxBinned->data.F32[i] - Fmax);
+	      psAbort ("impossible condition f[0] < Fo");
 	    } else {
 		Ro = radiusBinned->data.F32[i-1] + (radiusBinned->data.F32[i] - radiusBinned->data.F32[i-1]) * (Fo - fluxBinned->data.F32[i-1]) / (fluxBinned->data.F32[i] - fluxBinned->data.F32[i-1]);
@@ -118,5 +122,6 @@
     }
 
-    // XXX call a visual function to show the unbinned & binned vectors + result radius
+    // show the results
+    // psphotPetrosianVisualProfileRadii (radius, flux, radiusBinned, fluxBinned, Ro);
 
     psFree(fluxBinned);