Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c
===================================================================
--- /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c	(revision 25127)
+++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c	(revision 25128)
@@ -13,7 +13,7 @@
 
   int Nsec = 24;
-  float Rmax = 64;
+  float Rmax = 200;
   float fluxMin = 0.0;
-  float fluxMax = 1000.0;
+  float fluxMax = source->peak->flux;
 
   // generate a series of radial profiles at Nsec evenly spaced angles.  the profile flux
Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c
===================================================================
--- /branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c	(revision 25127)
+++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c	(revision 25128)
@@ -1,3 +1,12 @@
 # include "psphotInternal.h"
+
+# define DX 512
+# define DY 512
+
+// XXX add noise and seeing.
+// XXX double check on sersic functional form
+// XXX modify ratio if ratio > 1.0 (swap major and minor)
+
+pmPeak *psphotLocalPeak(pmReadout *readout, int Xo, int Yo);
 
 int main (int argc, char **argv) {
@@ -7,8 +16,8 @@
 
   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 Xo = 0.5*DX;
+  float Yo = 0.5*DY;
+  char *image = NULL;
+  pmSource *source = NULL;
 
   float peak = 1000.0;
@@ -18,5 +27,5 @@
   float sersic = 0.5;
   float skynoise = 0.0;
-
+  
   if ((N = psArgumentGet (argc, argv, "-peak"))) {
     psArgumentRemove (N, &argc, argv);
@@ -53,4 +62,16 @@
     pmVisualSetVisual(true);
   }
+  if ((N = psArgumentGet (argc, argv, "-coords"))) {
+    psArgumentRemove (N, &argc, argv);
+    Xo = atof(argv[N]);
+    psArgumentRemove (N, &argc, argv);
+    Yo = atof(argv[N]);
+    psArgumentRemove (N, &argc, argv);
+  }
+  if ((N = psArgumentGet (argc, argv, "-image"))) {
+    psArgumentRemove (N, &argc, argv);
+    image = psStringCopy (argv[N]);
+    psArgumentRemove (N, &argc, argv);
+  }
 
   if (argc != 1) {
@@ -61,46 +82,57 @@
   // create a containing image & associated readout
   pmReadout *readout = pmReadoutAlloc(NULL);
-  readout->image = psImageAlloc(128, 128, PS_TYPE_F32);
+  if (!image) {
+      readout->image = psImageAlloc(DX, DY, 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 dummy variance, but don't populate -- it is not used by pmSourceMoments if the sigma parameters is set to 0.0
+      readout->variance = psImageAlloc(DX, DY, PS_TYPE_F32);
 
-  // create a model & associated source
-  pmModelType type = pmModelClassGetType("PS_MODEL_SERSIC");
-  pmModel *model = pmModelAlloc(type);
+      // 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;
+      // 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] = Xo;
+      model->params->data.F32[PM_PAR_YPOS] = Yo;
 
-  psEllipseAxes axes;
-  axes.major = sigma;
-  axes.minor = sigma*ARatio;
-  axes.theta = angle;
+      psEllipseAxes axes;
+      axes.major = sigma;
+      axes.minor = sigma*ARatio;
+      axes.theta = angle;
 
-  psEllipseShape shape = psEllipseAxesToShape (axes);
+      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;
+      // 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;
+      if (model->params->n > 7) {
+	  model->params->data.F32[PM_PAR_7]  = sersic;
+      }
+
+      // generate source container & populate image
+      source = pmSourceFromModel(model, readout, Xo, 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");
+
+  } else {
+      psRegion full = psRegionSet(0,0,0,0);
+      psFits *fits = psFitsOpen(image, "r");
+      readout->image = psFitsReadImage(fits, full, 0);
+
+      source = pmSourceAlloc();
+      source->peak = psphotLocalPeak(readout, Xo, Yo);
+      pmSourceDefinePixels (source, readout, Xo, Yo, 64);
   }
-
-  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);
@@ -110,2 +142,43 @@
   exit (0);
 }
+
+// Xo, Yo are in parent coords
+pmPeak *psphotLocalPeak(pmReadout *readout, int Xo, int Yo) {
+
+    int Xp = Xo;
+    int Yp = Yo;
+    float peakFlux = readout->image->data.F32[Yp][Xp];
+
+    // find local peak within +/- 3 pix of the given coordinate
+    for (int iy = Yo - 3; iy <= Yo + 3; iy++) {
+	for (int ix = Xo - 3; ix <= Xo + 3; ix++) {
+	    if (peakFlux < readout->image->data.F32[iy][ix]) {
+		Xp = ix;
+		Yp = iy;
+		peakFlux = readout->image->data.F32[Yp][Xp];
+	    }
+	}
+    }
+
+    pmPeak *peak = pmPeakAlloc(Xp, Yp, peakFlux, PM_PEAK_LONE);
+
+    // calculate fractional peak position relative to Xp,Yp
+    psPolynomial2D *bicube = psImageBicubeFit (readout->image, Xp, Yp);
+    psPlane min = psImageBicubeMin (bicube);
+    psFree (bicube);
+
+    // if min point is too deviant, use the peak value
+    if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
+        peak->xf = min.x + Xp;
+        peak->yf = min.y + Yp;
+
+	// xf,yf must land on image with 0 pixel border
+	peak->xf = PS_MAX (PS_MIN (peak->xf, readout->image->numCols - 1), readout->image->col0);
+	peak->yf = PS_MAX (PS_MIN (peak->yf, readout->image->numRows - 1), readout->image->row0);
+    } else {
+        peak->xf = Xp;
+        peak->yf = Yp;
+    }
+
+    return peak;
+}
Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c
===================================================================
--- /branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c	(revision 25127)
+++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c	(revision 25128)
@@ -26,4 +26,5 @@
 	graphdata->xmin = graphdata->xmax = xVec->data.F32[0];
 	for (int i = 1; i < xVec->n; i++) {
+	    if (!isfinite(xVec->data.F32[i])) continue;
 	    graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]);
 	    graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]);
@@ -36,4 +37,5 @@
 	graphdata->ymin = graphdata->ymax = yVec->data.F32[0];
 	for (int i = 1; i < yVec->n; i++) {
+	    if (!isfinite(yVec->data.F32[i])) continue;
 	    graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]);
 	    graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]);
Index: /branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c
===================================================================
--- /branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c	(revision 25127)
+++ /branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c	(revision 25128)
@@ -5,9 +5,9 @@
 
 // XXX choices for rMax and Nsec?
-// XXX where does the output go?
-// XXX how does this fit in context?
 
 bool psphotRadialProfilesByAngles (pmPetrosian *petrosian, pmSource *source, int Nsec, float Rmax) {
 
+    // we want to have an even number of sectors so we can do 180 deg symmetrizing
+    Nsec = (Nsec % 2) ? Nsec + 1 : Nsec;
     float dtheta = 2.0*M_PI / Nsec;
 
@@ -46,5 +46,7 @@
 	    // value is NAN if we run off the image
 	    float value = psImageInterpolatePixelBilinear(x, y, source->pixels);
-	    if (isnan(value)) continue;
+
+	    // keep the nan values so all vectors are matched
+	    // if (isnan(value)) continue;
 
 	    psVectorAppend (radius, r);
@@ -62,4 +64,41 @@
     }
 
+    for (int i = 0; i < Nsec / 2; i++) {
+
+	psVector *r1 = petrosian->radii->data[i];
+	psVector *r2 = petrosian->radii->data[i+Nsec/2];
+
+	psVector *f1 = petrosian->fluxes->data[i];
+	psVector *f2 = petrosian->fluxes->data[i+Nsec/2];
+
+	psAssert (r1->n == r2->n, "mis-matched vectors");
+	psAssert (f1->n == f2->n, "mis-matched vectors");
+
+	// we have a pair of vectors i, i+Nsec/2; replace them with the finite minimum of the pair
+	for (int j = 0; j < r1->n; j++) {
+	    
+	    float flux;
+
+	    if (!isfinite(f1->data.F32[j]) && !isfinite(f2->data.F32[j])) {
+		flux = NAN;
+		goto setflux;
+	    }
+
+	    if (!isfinite(f1->data.F32[j])) {
+		flux = f2->data.F32[j];
+		goto setflux;
+	    }
+	    if (!isfinite(f2->data.F32[j])) {
+		flux = f1->data.F32[j];
+		goto setflux;
+	    }
+
+	    flux = PS_MIN(f1->data.F32[j], f2->data.F32[j]);
+
+	setflux:
+	    f1->data.F32[j] = flux;
+	    f2->data.F32[j] = flux;
+	}
+    }    
     return true;
 }
Index: /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c
===================================================================
--- /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c	(revision 25127)
+++ /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c	(revision 25128)
@@ -30,25 +30,29 @@
 
     // examine data in the two ranges Fm - Fo and Fo - Fp to define the bin size
+    // XXX reconsider the fractional isophote value
     float Fm = fluxMin + 0.25*fluxRange;
     float Fp = fluxMin + 0.75*fluxRange;
     float Fo = fluxMin + 0.50*fluxRange;
+    int Rbin = 1;
       
-    // find the mean radius of the points in the flux range Fm - Fp:
-    float Rsum = 0;
-    float Rnpt = 0;
-    for (int i = 0; i < flux->n; i++) {
-	if (flux->data.F32[i] < Fm) continue;
-	if (flux->data.F32[i] > Fp) continue;
-      
-	Rsum += radius->data.F32[i];
-	Rnpt ++;
-    }
+    // find the median radius of the points in the flux range Fm - Fp:
+    { 
+	// storage vector for stats
+	psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
+	psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
 
-    // rebin with bins 1/2 the size of the mean radius
-    int Rbin = 0;
-    if (Rnpt == 0) {
-	Rbin = 1;
-    } else {
-	Rbin = MAX(1, 0.5*(Rsum / Rnpt));
+	for (int i = 0; i < flux->n; i++) {
+	    if (!isfinite(flux->data.F32[i])) continue;
+	    if (flux->data.F32[i] < Fm) continue;
+	    if (flux->data.F32[i] > Fp) continue;
+	    
+	    psVectorAppend (values, radius->data.F32[i]);
+	}
+	psVectorStats (stats, values, NULL, NULL, 0);
+
+	// if we have a valid range, rebin with bin size 1/2 of median radius
+	if (isfinite(stats->sampleMedian)) {
+	    Rbin = MAX(1, 0.5*stats->sampleMedian);
+	}
     }
 
@@ -94,4 +98,5 @@
 		psStatsInit(stats);
 	    }
+	    if (!isfinite(flux->data.F32[i])) continue;
 	    psVectorAppend (values, flux->data.F32[i]);
 	}
@@ -105,4 +110,6 @@
     bool above = true;
     for (int i = 0; i < fluxBinned->n; i++) {
+
+	if (!isfinite(fluxBinned->data.F32[i])) continue;
 
 	// find the largest radius that matches the flux transition