Index: /branches/eam_branches/psphot.20100506/src/Makefile.am
===================================================================
--- /branches/eam_branches/psphot.20100506/src/Makefile.am	(revision 27931)
+++ /branches/eam_branches/psphot.20100506/src/Makefile.am	(revision 27932)
@@ -191,4 +191,5 @@
         psphotEllipticalProfile.c      \
 	psphotRadialBins.c	       \
+	psphotRadialApertures.c	       \
 	psphotPetrosian.c	       \
         psphotPetrosianRadialBins.c    \
Index: /branches/eam_branches/psphot.20100506/src/psphot.h
===================================================================
--- /branches/eam_branches/psphot.20100506/src/psphot.h	(revision 27931)
+++ /branches/eam_branches/psphot.20100506/src/psphot.h	(revision 27932)
@@ -458,3 +458,7 @@
 bool psphotCopySourcesReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
 
+bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax);
+
 #endif
Index: /branches/eam_branches/psphot.20100506/src/psphotRadialApertures.c
===================================================================
--- /branches/eam_branches/psphot.20100506/src/psphotRadialApertures.c	(revision 27932)
+++ /branches/eam_branches/psphot.20100506/src/psphotRadialApertures.c	(revision 27932)
@@ -0,0 +1,208 @@
+# include "psphotInternal.h"
+
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    // 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;
+    }
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+	if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i);
+	    return false;
+	}
+    }
+    return true;
+}
+
+// aperture-like measurements for extended sources
+// flux in simple, circular apertures
+bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+
+    bool status;
+    int Nradial = 0;
+
+    psTimerStart ("psphot.radial");
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
+
+    if (!sources->n) {
+	psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
+	return true;
+    }
+
+    // radMax stores the upper bounds of the annuli
+    // XXX keep the same name here as for the petrosian / elliptical apertures?
+    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    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");
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    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");
+
+    // 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);
+
+    // option to limit analysis to a specific region
+    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
+    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
+
+    // choose the sources of interest
+    for (int i = 0; i < sources->n; i++) {
+
+	pmSource *source = sources->data[i];
+
+	// skip PSF-like and non-astronomical objects
+	if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+	if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+	if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
+	if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+
+	// limit selection to some SN limit
+	assert (source->peak); // how can a source not have a peak?
+	if (source->peak->SN < SN_LIM) continue;
+
+	// limit selection by analysis region
+	if (source->peak->x < AnalysisRegion.x0) continue;
+	if (source->peak->y < AnalysisRegion.y0) continue;
+	if (source->peak->x > AnalysisRegion.x1) continue;
+	if (source->peak->y > AnalysisRegion.y1) continue;
+
+	// replace object in image
+	if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+	    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+	}
+	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)) {
+	    psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+	} else {
+	    source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
+	}
+
+	// re-subtract the object, leave local sky
+	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+    }
+
+    psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
+    return true;
+}
+
+bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax) {
+
+    // allocate pmSourceExtendedParameters, if not already defined
+    if (!source->radial) {
+        source->radial = pmSourceRadialAperturesAlloc ();
+    }
+    
+    psVector *radius  = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixFlux = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixVar  = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+    for (int iy = 0; iy < source->pixels->numRows; iy++) {
+	for (int ix = 0; ix < source->pixels->numCols; ix++) {
+
+	    // 0.5 PIX: get radius as a function of pixel coord
+	    float x = ix + 0.5 - source->peak->xf + source->pixels->col0;
+	    float y = iy + 0.5 - source->peak->yf + source->pixels->row0;
+
+	    float r = hypot(x, y);
+
+	    psVectorAppend(radius, r);
+	    psVectorAppend(pixFlux, source->pixels->data.F32[iy][ix]);
+	    psVectorAppend(pixVar, source->variance->data.F32[iy][ix]);
+	}
+    }
+
+    psVector *flux    = psVectorAllocEmpty(radMax->n, PS_TYPE_F32); // surface brightness of radial bin
+    psVector *fluxErr = psVectorAllocEmpty(radMax->n, PS_TYPE_F32); // surface brightness of radial bin
+    psVector *fill    = psVectorAllocEmpty(radMax->n, PS_TYPE_F32); // surface brightness of radial bin
+
+    psVectorInit (flux,    0.0);
+    psVectorInit (fluxErr, 0.0);
+    psVectorInit (fill,    0.0);
+
+    float fluxSum = 0.0;
+    float varSum = 0.0;
+    int nPixSum = 0;
+
+    bool done = false;
+    int nOut = 0;
+    float Rmax = radMax->data.F32[nOut];
+
+    // XXX assume (or enforce) that the bins are contiguous and non-overlapping (Rmax[i] = Rmin[i+1])
+    for (int i = 0; !done && (i < radius->n); i++) {
+	if (radius->data.F32[i] > Rmax) {
+	    // calculate the total flux for bin 'nOut'
+	    float Area = M_PI*PS_SQR(Rmax);
+	    flux->data.F32[nOut] = fluxSum;
+	    fluxErr->data.F32[nOut] = sqrt(varSum);
+	    fill->data.F32[nOut] = nPixSum / Area;
+
+	    psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %4.2f %6.1f\n", 
+		     nOut, radMax->data.F32[nOut], flux->data.F32[nOut], fluxErr->data.F32[nOut], fill->data.F32[nOut], Area);
+
+	    nOut ++;
+	    if (nOut >= radMax->n) break;
+	    Rmax = radMax->data.F32[nOut];
+	}
+	fluxSum += pixFlux->data.F32[i];
+	varSum += pixVar->data.F32[i];
+	nPixSum ++;
+    }
+    flux->n = fluxErr->n = fill->n = nOut;
+    
+    psFree(source->radial->flux);
+    psFree(source->radial->fluxErr);
+    psFree(source->radial->fill);
+
+    source->radial->flux = flux;
+    source->radial->fluxErr = fluxErr;
+    source->radial->fill = fill;
+
+    psFree (radius);
+    psFree (pixFlux);
+    psFree (pixVar);
+
+    return true;
+}
Index: /branches/eam_branches/psphot.20100506/src/psphotRadialBins.c
===================================================================
--- /branches/eam_branches/psphot.20100506/src/psphotRadialBins.c	(revision 27931)
+++ /branches/eam_branches/psphot.20100506/src/psphotRadialBins.c	(revision 27932)
@@ -145,5 +145,5 @@
 	}
     }
-    binSB->n = binSBstdev->n = binRad->n = binArea->n = nOut;
+    binSB->n = binSBstdev->n = binSum->n = binRad->n = binArea->n = nOut;
 
     // interpolate any bins that were empty (extrapolate to center if needed)
Index: /branches/eam_branches/psphot.20100506/src/psphotStackReadout.c
===================================================================
--- /branches/eam_branches/psphot.20100506/src/psphotStackReadout.c	(revision 27931)
+++ /branches/eam_branches/psphot.20100506/src/psphotStackReadout.c	(revision 27932)
@@ -161,4 +161,5 @@
     }
 
+    psphotRadialApertures (config, view, STACK_OUT); 
     psphotExtendedSourceAnalysis (config, view, STACK_OUT); // pass 1 (detections->allSources)
     psphotExtendedSourceFits (config, view, STACK_OUT); // pass 1 (detections->allSources)