Index: trunk/psphot/src/psphotPetrosianStats.c
===================================================================
--- trunk/psphot/src/psphotPetrosianStats.c	(revision 25755)
+++ trunk/psphot/src/psphotPetrosianStats.c	(revision 27819)
@@ -10,8 +10,9 @@
 bool psphotPetrosianStats (pmSource *source) {
 
-    pmSourceRadialProfile *profile = source->extpars->profile;
-
-    float petRadius = NAN;
-    float petFlux = NAN;
+    psAssert (source, "missing source");
+    psAssert (source->extpars, "missing extpars");
+    psAssert (source->extpars->petProfile, "missing petProfile");
+
+    pmSourceRadialProfile *profile = source->extpars->petProfile;
 
     psVector *binSB      = profile->binSB;
@@ -28,4 +29,7 @@
     psVector *areaSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
 
+    float petRadius = NAN;
+    float petFlux = NAN;
+
     bool anyPetro = false;
     bool manyPetro = false;
@@ -38,4 +42,6 @@
     int lowestSignificantRadius = 0;
     float lowestSignificantRatio = 1.0;
+
+    // find the Petrosian Radius and Petrosian Flux
 
     int nOut = 0;
@@ -142,14 +148,48 @@
     }
 
-    if (!source->extpars->petrosian_80) {
-        source->extpars->petrosian_80 = pmSourceExtendedFluxAlloc ();
-    }
-    pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80;
+    // now measure the radii R90 and R50 where flux = 0.9 (or 0.5) * petFlux;
+    float flux90 = 0.9 * petFlux;
+    float flux50 = 0.5 * petFlux;
+    float R50 = NAN;
+    float R90 = NAN;
+    bool found50 = false;
+    bool found90 = false;
+    // XXX use bisection to do this faster:
+    for (int i = 0; !(found50 && found90) && i < refRadius->n; i++) {
+	if (!found50 && (fluxSum->data.F32[i] > flux50)) {
+	    if (i == 0) {
+		psWarning ("does this case make any sense? (fluxSum[0] > flux50)");
+		continue;
+	    } else {
+		R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
+		found50 = true;
+	    }
+	}
+	if (!found90 && (fluxSum->data.F32[i] > flux90)) {
+	    if (i == 0) {
+		psWarning ("does this case make any sense? (fluxSum[0] > flux90)");
+		continue;
+	    } else {
+		R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);
+		found90 = true;
+	    }
+	}
+    }
+
 
     // XXX save flags (anyPetro, manyPetro)
-    petrosian->radius = petRadius;
-    petrosian->flux   = petFlux;
-
-    // psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);
+    source->extpars->petrosianRadius = petRadius;
+    source->extpars->petrosianFlux   = petFlux;
+    source->extpars->petrosianR50    = R50;
+    source->extpars->petrosianR90    = R90;
+    
+    // XXX add the errors
+    source->extpars->petrosianRadiusErr = NAN;
+    source->extpars->petrosianFluxErr   = NAN;
+    source->extpars->petrosianR50Err    = NAN;
+    source->extpars->petrosianR90Err    = NAN;
+
+    fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf);
+    psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);
 
     psFree(fluxSum);