Index: trunk/psphot/src/psphotPetrosian.c
===================================================================
--- trunk/psphot/src/psphotPetrosian.c	(revision 21183)
+++ trunk/psphot/src/psphotPetrosian.c	(revision 25755)
@@ -1,109 +1,32 @@
 # include "psphotInternal.h"
 
-bool psphotPetrosian (pmSource *source, psMetadata *recipe, psImageMaskType maskVal) {
+bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) {
 
-  bool status;
+    // XXX these need to go into recipe values
+    float Rmax = 200;
 
-  assert (source->extpars);
-  assert (source->extpars->profile);
-  assert (source->extpars->profile->radius);
-  assert (source->extpars->profile->flux);
+    psAssert (source->extpars, "need to run psphotRadialProfile first");
+    psAssert (source->extpars->profile, "need to run psphotRadialProfile first");
 
-  psVector *radius = source->extpars->profile->radius;
-  psVector *flux = source->extpars->profile->flux;
+    // 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, Rmax, skynoise)) {
+	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 (source)) {
+	psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
+	return false;
+    }
+  
+    psTrace ("psphot", 3, "source at %f,%f: petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f",
+	     source->peak->xf, source->peak->yf, 
+	     source->extpars->petrosian_80->radius, 
+	     source->extpars->petrosian_80->flux, 
+	     source->extpars->profile->axes.minor/source->extpars->profile->axes.major, 
+	     source->extpars->profile->axes.theta*PS_DEG_RAD);
 
-  // flux at which to measure isophotal parameters
-  float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
-  float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");
-  assert (status);
-
-  // first find flux at R0
-  int firstAbove = -1;
-  int lastBelow = -1;
-  for (int i = 0; i < radius->n; i++) {
-    if (radius->data.F32[i] < PETROSIAN_R0) lastBelow = i;
-    if ((firstAbove < 0) && (radius->data.F32[i] > PETROSIAN_R0)) firstAbove = i;
-  }
-  // if we don't go out far enough, we have a problem...
-  if (lastBelow == radius->n - 1) {
-    psTrace ("psphot", 5, "did not go out far enough to reach petrosian reference radius...");
-    // XXX skip object? raise a flag ?
-    return false;
-  }
-  if (firstAbove < 0) {
-    psTrace ("psphot", 5, "did not go out far enough to bound petrosian reference radius");
-    // XXX raise a flag ?
-    return false;
-  }
-
-  // average flux in this range
-  float fluxR0 = 0.0;
-  int fluxRn = 0;
-  for (int i = PS_MIN(firstAbove, lastBelow); i <= PS_MAX(firstAbove, lastBelow); i++) {
-    fluxR0 += flux->data.F32[i];
-    fluxRn ++;
-  }
-  fluxR0 /= (float)(fluxRn);
-
-  // target flux for petrosian radius
-  float fluxRP = fluxR0 * PETROSIAN_RF;
-
-  // find the first bin below the flux level and the last above the level
-  // XXX can this be done faster with bisection?
-  // XXX do I need to worry about crazy outliers?
-  // XXX should i be smoothing or fitting the curve?
-  int firstBelow = -1;
-  int lastAbove = -1;
-  for (int i = 0; i < flux->n; i++) {
-    if (flux->data.F32[i] > fluxRP) lastAbove = i;
-    if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i;
-  }
-  // if we don't go out far enough, we have a problem...
-  if (lastAbove == radius->n - 1) {
-    psTrace ("psphot", 5, "did not go out far enough to reach petrosian radius...");
-    // XXX skip object? raise a flag ?
-    return false;
-  }
-  if (firstBelow < 0) {
-    psTrace ("psphot", 5, "did not go out far enough to bound petrosian radius");
-    // XXX raise a flag ?
-    return false;
-  }
-
-  // need to examine pixels in this vicinity
-  float fluxFirst = 0;
-  float fluxLast = 0;
-  for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) {
-    if (i <= firstBelow) {
-      fluxFirst += flux->data.F32[i];
-    }
-    if (i <= lastAbove) {
-      fluxLast += flux->data.F32[i];
-    }
-  }
-  float fluxRPSum    = 0.5*(fluxLast + fluxFirst);
-  float fluxRPSumErr = 0.5*fabs(fluxLast - fluxFirst);
-  // XXX need to use the weight appropriately here...
-
-  float rad     = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]);
-  float radErr  = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]);
-
-  if (!source->extpars->petrosian) {
-    source->extpars->petrosian = pmSourcePetrosianValuesAlloc ();
-  }
-
-  // these are uncalibrated: instrumental mags and pixel units
-  source->extpars->petrosian->mag    = -2.5*log10(fluxRPSum);
-  source->extpars->petrosian->magErr = fluxRPSumErr / fluxRPSum;
-
-  source->extpars->petrosian->rad    = rad;
-  source->extpars->petrosian->radErr = radErr;
-
-  psTrace ("psphot", 5, "Petrosian flux:%f +/- %f @ %f +/- %f for %f, %f\n",
-           source->extpars->petrosian->mag, source->extpars->petrosian->magErr,
-           source->extpars->petrosian->rad, source->extpars->petrosian->radErr,
-           source->peak->xf, source->peak->yf);
-
-  return true;
-
+    return true;
 }