Index: trunk/psphot/src/Makefile.am
===================================================================
--- trunk/psphot/src/Makefile.am	(revision 15405)
+++ trunk/psphot/src/Makefile.am	(revision 15562)
@@ -47,4 +47,5 @@
 	psphotModelWithPSF.c     \
 	psphotExtendedSources.c	 \
+	psphotRadialProfile.c	 \
 	psphotPetrosian.c	 \
 	psphotIsophotal.c	 \
Index: trunk/psphot/src/psphot.h
===================================================================
--- trunk/psphot/src/psphot.h	(revision 15405)
+++ trunk/psphot/src/psphot.h	(revision 15562)
@@ -116,4 +116,5 @@
 psKernel       *psphotKernelFromPSF (pmSource *source, int nPix);
 
+bool            psphotRadialProfile (pmSource *source, psMetadata *recipe, psMaskType maskVal);
 bool 		psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal);
 bool 		psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal);
Index: trunk/psphot/src/psphotAnnuli.c
===================================================================
--- trunk/psphot/src/psphotAnnuli.c	(revision 15405)
+++ trunk/psphot/src/psphotAnnuli.c	(revision 15562)
@@ -3,6 +3,57 @@
 bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
 
-  psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
+  assert (source->extpars);
+  assert (source->extpars->profile);
+  assert (source->extpars->profile->radius);
+  assert (source->extpars->profile->flux);
+
+  bool status;
+
+  psVector *radius = source->extpars->profile->radius;
+  psVector *weight = source->extpars->profile->weight;
+  psVector *flux = source->extpars->profile->flux;
+
+  // XXX how do I define the radii?  we can put a vector in the recipe...
+  // radialBins defines the bounds or start and stop (we can skip some that way...
+  psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
+  psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+  assert (radialBinsLower->n == radialBinsUpper->n);
+
+  psVector *fluxValues = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
+  psVector *fluxSquare = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
+  psVector *fluxWeight = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
+  psVector *pixelCount = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
+  psVectorInit (fluxValues, 0.0);
+  psVectorInit (fluxSquare, 0.0);
+  psVectorInit (fluxWeight, 0.0);
+  psVectorInit (pixelCount, 0.0);
+
+  // XXX this code assumes the radii are in pixels.  convert from arcsec with plate scale
+  // XXX assume the annulii above are not overlapping?  much faster...
+  // XXX this might be must faster in the reverse order: loop over annulii and use disection to 
+  // skip to the start of the annulus.
+  for (int i = 0; i < flux->n; i++) {
+    for (int j = 0; j < radialBinsLower->n; j++) {
+      if (radius->data.F32[i] < radialBinsLower->data.F32[j]) continue;
+      if (radius->data.F32[i] > radialBinsUpper->data.F32[j]) continue;
+      fluxValues->data.F32[j] += flux->data.F32[i];
+      fluxSquare->data.F32[j] += PS_SQR(flux->data.F32[i]);
+      fluxWeight->data.F32[j] += weight->data.F32[i];
+      pixelCount->data.F32[j] += 1.0;
+    }
+  }
+
+  for (int j = 0; j < radialBinsLower->n; j++) {
+    fluxValues->data.F32[j] /= pixelCount->data.F32[j];
+    fluxSquare->data.F32[j] /= pixelCount->data.F32[j];
+    fluxSquare->data.F32[j] -= PS_SQR(fluxValues->data.F32[j]);
+  }
+  
+  source->extpars->annuli = pmSourceAnnuliAlloc ();
+  source->extpars->annuli->flux    = fluxValues;
+  source->extpars->annuli->fluxErr = fluxWeight;
+  source->extpars->annuli->fluxVar = fluxSquare;
+
   return true;
+}
 
-}
Index: trunk/psphot/src/psphotExtendedSources.c
===================================================================
--- trunk/psphot/src/psphotExtendedSources.c	(revision 15405)
+++ trunk/psphot/src/psphotExtendedSources.c	(revision 15562)
@@ -66,7 +66,18 @@
 	    return false;
 	} else {
+	  // XXX why am I caching the model?
 	    pmSourceCacheModel (source, maskVal); // XXX put this in the source model function?
 	    psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
             Npsf ++;
+	}
+
+	// all of these below require the radial profile
+	// XXX push this as a test and call in each of the functions below?
+	// XXX this constructs a pmSourceExtendedParameters element
+	if (doPetrosian || doIsophotal || doAnnuli || doKron) {
+	  if (!psphotRadialProfile (source, recipe, maskVal)) {
+	    psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate radial profile");
+	    return false;
+	  }
 	}
 
Index: trunk/psphot/src/psphotIsophotal.c
===================================================================
--- trunk/psphot/src/psphotIsophotal.c	(revision 15405)
+++ trunk/psphot/src/psphotIsophotal.c	(revision 15562)
@@ -3,5 +3,57 @@
 bool psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
 
-  psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
+  assert (source->extpars);
+  assert (source->extpars->profile);
+  assert (source->extpars->profile->radius);
+  assert (source->extpars->profile->flux);
+
+  bool status;
+
+  psVector *radius = source->extpars->profile->radius;
+  psVector *flux = source->extpars->profile->flux;
+
+  // flux at which to measure isophotal parameters
+  // XXX cache this?
+  float ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX");
+
+  // find the first bin below the flux level and the last above the level
+  // XXX can this be done faster with disection?
+  // 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] > ISOPHOT_FLUX) lastAbove = i;
+    if ((firstBelow < 0) && (flux->data.F32[i] < ISOPHOT_FLUX)) firstBelow = i;
+  }
+
+  // need to examine pixels in this vicinity
+  float isophotalFluxFirst = 0;
+  float isophotalFluxLast = 0;
+  for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) {
+    if (i <= firstBelow) {
+      isophotalFluxFirst += flux->data.F32[i];
+    }
+    if (i <= lastAbove) {
+      isophotalFluxLast += flux->data.F32[i];
+    }
+  }
+  float isophotalFlux    = 0.5*(isophotalFluxLast + isophotalFluxFirst);
+  float isophotalFluxErr = 0.5*fabs(isophotalFluxLast - isophotalFluxFirst);
+
+  float isophotalRad     = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]);
+  float isophotalRadErr  = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]);
+
+  if (!source->extpars->isophot) {
+    source->extpars->isophot = pmSourceIsophotalValuesAlloc ();
+  }
+  
+  // these are uncalibrated: instrumental mags and pixel units
+  source->extpars->isophot->mag    = -2.5*log10(isophotalFlux);
+  source->extpars->isophot->magErr = isophotalFluxErr / isophotalFlux;
+
+  source->extpars->isophot->rad    = isophotalRad;
+  source->extpars->isophot->radErr = isophotalRadErr;
+
   return true;
 
Index: trunk/psphot/src/psphotMakeResiduals.c
===================================================================
--- trunk/psphot/src/psphotMakeResiduals.c	(revision 15405)
+++ trunk/psphot/src/psphotMakeResiduals.c	(revision 15562)
@@ -28,4 +28,7 @@
 
     float nSigma = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.NSIGMA");
+    PS_ASSERT (status, false);
+
+    float pixelSN = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.PIX.SN");
     PS_ASSERT (status, false);
 
@@ -194,4 +197,9 @@
                 resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
                 //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
+
+		if (resid->Ro->data.F32[oy][ox] < pixelSN*fluxStats->sampleStdev) {
+		  resid->mask->data.U8[oy][ox] = 1;
+		}
+
             } else {
                 assert (SPATIAL_ORDER == 1);
@@ -227,4 +235,9 @@
                 resid->Rx->data.F32[oy][ox] = B->data.F64[1];
                 resid->Ry->data.F32[oy][ox] = B->data.F64[2];
+
+		float dRo = sqrt(A->data.F32[0][0]);
+		if (resid->Ro->data.F32[oy][ox] < pixelSN*dRo) {
+		  resid->mask->data.U8[oy][ox] = 1;
+		}
                 //resid->weight->data.F32[oy][ox] = XXX;
             }
Index: trunk/psphot/src/psphotPetrosian.c
===================================================================
--- trunk/psphot/src/psphotPetrosian.c	(revision 15405)
+++ trunk/psphot/src/psphotPetrosian.c	(revision 15562)
@@ -3,5 +3,79 @@
 bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
 
-  psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
+  bool status;
+
+  assert (source->extpars);
+  assert (source->extpars->profile);
+  assert (source->extpars->profile->radius);
+  assert (source->extpars->profile->flux);
+
+  psVector *radius = source->extpars->profile->radius;
+  psVector *flux = source->extpars->profile->flux;
+
+  // flux at which to measure isophotal parameters
+  // XXX cache this?
+  float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
+  float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");
+
+  // first find flux at R0
+  int firstBelow = -1;
+  int lastAbove = -1;
+  for (int i = 0; i < radius->n; i++) {
+    if (radius->data.F32[i] > PETROSIAN_R0) lastAbove = i;
+    if ((firstBelow < 0) && (flux->data.F32[i] < PETROSIAN_R0)) firstBelow = i;
+  }
+
+  // average flux in this range
+  float fluxR0 = 0.0;
+  int fluxRn = 0;
+  for (int i = PS_MIN(firstBelow, lastAbove); i <= PS_MAX(firstBelow, lastAbove); 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 disection?
+  // XXX do I need to worry about crazy outliers?  
+  // XXX should i be smoothing or fitting the curve?
+  firstBelow = -1;
+  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;
+  }
+
+  // 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;
+
   return true;
 
Index: trunk/psphot/src/psphotRadialProfile.c
===================================================================
--- trunk/psphot/src/psphotRadialProfile.c	(revision 15405)
+++ trunk/psphot/src/psphotRadialProfile.c	(revision 15562)
@@ -5,9 +5,9 @@
     // allocate pmSourceExtendedParameters, if not already defined
     if (!source->extpars) {
-	source->extpars = pmSourceExtendedParametersAlloc ();
+	source->extpars = pmSourceExtendedParsAlloc ();
     }
 
     if (!source->extpars->profile) {
-	source->extpars->profile = pmSourceRadialProfile (); 
+	source->extpars->profile = pmSourceRadialProfileAlloc (); 
     }    
     
@@ -19,5 +19,5 @@
     psVector *radius = source->extpars->profile->radius;
     psVector *flux   = source->extpars->profile->flux;
-    psVector *weight = source->extpars->profile->radius;
+    psVector *weight = source->extpars->profile->weight;
 
     // XXX use the extended source model here for Xo, Yo?
@@ -40,5 +40,6 @@
     flux->n = n;
 
-    SortVectorsByRadius (radius, flux, weight);
+    // XXX need to sort here
+    // SortVectorsByRadius (radius, flux, weight);
 
     return true;