Index: trunk/psphot/src/psphotRadialApertures.c
===================================================================
--- trunk/psphot/src/psphotRadialApertures.c	(revision 30624)
+++ trunk/psphot/src/psphotRadialApertures.c	(revision 31154)
@@ -37,4 +37,10 @@
     int Nradial = 0;
 
+    // perform full non-linear fits / extended source analysis?
+    if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) {
+	psLogMsg ("psphot", PS_LOG_INFO, "skipping radial apertures\n");
+	return true;
+    }
+
     psTimerStart ("psphot.radial");
 
@@ -62,4 +68,5 @@
     psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
     psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
+    float outerRadius = radMax->data.F32[radMax->n - 1];
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -67,7 +74,4 @@
     assert (maskVal);
 
-    // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
-    float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
-
     // S/N limit to perform full non-linear fits
     float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
@@ -75,5 +79,5 @@
     // source analysis is done in S/N order (brightest first)
     // XXX are we getting the objects out of order? does it matter?
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // option to limit analysis to a specific region
@@ -95,5 +99,5 @@
 	// limit selection to some SN limit
 	assert (source->peak); // how can a source not have a peak?
-	if (source->peak->SN < SN_LIM) continue;
+	if (sqrt(source->peak->detValue) < SN_LIM) continue;
 
 	// limit selection by analysis region
@@ -112,14 +116,18 @@
 	    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
 	}
+
+	// we need to change the view for the radial aperture analysis, but we want to recover exactly
+	// the original view; the following elements get destroyed by pmSourceRedefinePixels so save them:
+	psImage *oldMaskObj   = psMemIncrRefCounter(source->maskObj);
+	psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux);
+	psImage *oldPSFimage  = psMemIncrRefCounter(source->psfImage);
+	psRegion oldRegion    = source->region;
+
 	Nradial ++;
 
 	// force source image to be a bit larger...
-	float radius = source->peak->xf - source->pixels->col0;
-	radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
-	radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
-	radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
-	pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
-
-	if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, 0)) {
+	pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
+
+	if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, 0)) {
 	    psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
 	} else {
@@ -127,4 +135,9 @@
 	}
 
+	pmSourceRedefinePixelsByRegion (source, readout, oldRegion);
+	psFree(source->maskObj);   source->maskObj   = oldMaskObj;
+	psFree(source->modelFlux); source->modelFlux = oldModelFlux;
+	psFree(source->psfImage);  source->psfImage  = oldPSFimage;
+	
 	// re-subtract the object, leave local sky
 	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
@@ -135,5 +148,5 @@
 }
 
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
+bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
 
     // if we are a child source, save the results to the parent source radial aperture array
@@ -163,7 +176,12 @@
     }
 
-    // center of the apertures
-    float xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
-    float yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+    float xCM = NAN, yCM = NAN;
+    if (pmSourcePositionUseMoments(source)) {
+	xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
+	yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+    } else {
+	xCM = source->peak->xf - 0.5 - source->pixels->col0; // coord of peak in subimage
+	yCM = source->peak->yf - 0.5 - source->pixels->row0; // coord of peak in subimage
+    }
 
     // one pass through the pixels to select the valid pixels and calculate R^2
@@ -203,7 +221,9 @@
     psVector *flux    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
     psVector *fluxErr = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
+    psVector *fluxStd = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
     psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
 
     psVectorInit (flux,    0.0);
+    psVectorInit (fluxStd, 0.0);
     psVectorInit (fluxErr, 0.0);
     psVectorInit (fill,    0.0);
@@ -212,22 +232,40 @@
     for (int i = 0; i < pixRadius2->n; i++, rPix2++) {
 
+	int j = 0;
 	float *aRad2 = aperRadii2->data.F32;
-	for (int j = 0; (*rPix2 < *aRad2) && (j < aperRadii2->n); j++, aRad2++) {
+	for (; (*aRad2 < *rPix2) && (j < aperRadii2->n); j++, aRad2++);
+	for (; j < aperRadii2->n; j++, aRad2++) {
 	    flux->data.F32[j]    += pixFlux->data.F32[i];
+	    fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]);
 	    fluxErr->data.F32[j] += pixVar->data.F32[i];
 	    fill->data.F32[j]    += 1.0;
 	}
     }
+
+    /* for each radial bin, R(i), we measure:
+       1) the flux within that aperture: F(i) = \sum_{r_j<R_i}(F_j)
+       2) the fractional fill factor (count of valid pixels / effective area of the aperture
+       3) the error on the flux within that aperture
+     */
 
     for (int i = 0; i < flux->n; i++) {
 	// calculate the total flux for bin 'nOut'
 	float Area = M_PI*aperRadii2->data.F32[i];
-	fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
-	fill->data.F32[i] /= Area;
-	psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %4.2f %6.1f\n", 
-		 i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], fill->data.F32[i], Area);
+
+	int nPix = fill->data.F32[i];
+	float SBmean = flux->data.F32[i] / nPix;
+	float SBstdv = sqrt((fluxStd->data.F32[i] / nPix) - PS_SQR(SBmean));
+
+	// flux->data.F32[i]    = SBmean * Area;
+	fluxStd->data.F32[i] = SBstdv * Area;
+	fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
+
+	// fill->data.F32[i] /= Area;
+	psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %8.1f +/- %8.1f : %4.2f %6.1f\n", 
+		 i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], SBmean, SBstdv, fill->data.F32[i], Area);
     }
     
     radialAper->flux = flux;
+    radialAper->fluxStdev = fluxStd;
     radialAper->fluxErr = fluxErr;
     radialAper->fill = fill;
@@ -245,5 +283,5 @@
 static int nPix = 0;
 
-bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry) {
+bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry) {
 
     psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?");