Index: trunk/psphot/src/psphotRadialApertures.c
===================================================================
--- trunk/psphot/src/psphotRadialApertures.c	(revision 29936)
+++ trunk/psphot/src/psphotRadialApertures.c	(revision 30624)
@@ -102,4 +102,9 @@
 	if (source->peak->x > AnalysisRegion.x1) continue;
 	if (source->peak->y > AnalysisRegion.y1) continue;
+
+	// allocate pmSourceExtendedParameters, if not already defined
+	if (!source->radialAper) {
+	    source->radialAper = psArrayAlloc(1);
+	}
 
 	// replace object in image
@@ -116,5 +121,5 @@
 	pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
 
-	if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) {
+	if (!psphotRadialApertureSource (source, recipe, skynoise, 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 {
@@ -130,12 +135,122 @@
 }
 
-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 ();
+bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
+
+    // if we are a child source, save the results to the parent source radial aperture array
+    psArray *radialAperSet = source->radialAper;
+    if (source->parent) {
+	radialAperSet = source->parent->radialAper;
+    }
+    psAssert(radialAperSet, "this should be defined before calling");
+    psAssert(radialAperSet->data[entry] == NULL, "why is this already defined?");
+
+    pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc ();
+    radialAperSet->data[entry] = radialAper;
+
+    // storage for the derived pixel values
+    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+    // outer-most radius for initial truncation
+    float Rmax  = aperRadii->data.F32[aperRadii->n - 1];
+    float Rmax2 = PS_SQR(Rmax);
+
+    // store the R^2 values for the apertures
+    psVector *aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
+    for (int i = 0; i < aperRadii->n; i++) {
+	aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
+    }
+
+    // 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
+
+    // one pass through the pixels to select the valid pixels and calculate R^2
+    for (int iy = 0; iy < source->pixels->numRows; iy++) {
+
+	float yDiff = iy - yCM;
+	if (fabs(yDiff) > Rmax) continue;
+
+	float *vPix = source->pixels->data.F32[iy];
+	float *vWgt = source->variance->data.F32[iy];
+	psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy];
+
+	for (int ix = 0; ix < source->pixels->numCols; ix++, vPix++, vWgt++) {
+
+	    if (vMsk) {
+		if (*vMsk & maskVal) {
+		    vMsk++;
+		    continue;
+		}
+		vMsk++;
+	    }
+	    if (isnan(*vPix)) continue;
+
+	    float xDiff = ix - xCM;
+	    if (fabs(xDiff) > Rmax) continue;
+
+	    // radius is just a function of (xDiff, yDiff)
+	    float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+	    if (r2 > Rmax2) continue;
+
+	    psVectorAppend(pixRadius2, r2);
+	    psVectorAppend(pixFlux, *vPix);
+	    psVectorAppend(pixVar, *vWgt);
+	}
+    }
+
+    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 *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
+
+    psVectorInit (flux,    0.0);
+    psVectorInit (fluxErr, 0.0);
+    psVectorInit (fill,    0.0);
+
+    float *rPix2 = pixRadius2->data.F32;
+    for (int i = 0; i < pixRadius2->n; i++, rPix2++) {
+
+	float *aRad2 = aperRadii2->data.F32;
+	for (int j = 0; (*rPix2 < *aRad2) && (j < aperRadii2->n); j++, aRad2++) {
+	    flux->data.F32[j]    += pixFlux->data.F32[i];
+	    fluxErr->data.F32[j] += pixVar->data.F32[i];
+	    fill->data.F32[j]    += 1.0;
+	}
+    }
+
+    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);
     }
     
-    psVector *radius  = psVectorAllocEmpty(100, PS_TYPE_F32);
+    radialAper->flux = flux;
+    radialAper->fluxErr = fluxErr;
+    radialAper->fill = fill;
+
+    psFree (aperRadii2);
+    psFree (pixRadius2);
+    psFree (pixFlux);
+    psFree (pixVar);
+
+    return true;
+}
+
+static int nCalls = 0;
+static int nPass = 0;
+static int nPix = 0;
+
+bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry) {
+
+    psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?");
+
+    pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc ();
+    source->radialAper->data[entry] = radialAper;
+
+    psVector *pixRadius  = psVectorAllocEmpty(100, PS_TYPE_F32);
     psVector *pixFlux = psVectorAllocEmpty(100, PS_TYPE_F32);
     psVector *pixVar  = psVectorAllocEmpty(100, PS_TYPE_F32);
@@ -144,5 +259,5 @@
 	for (int ix = 0; ix < source->pixels->numCols; ix++) {
 
-	    // 0.5 PIX: get radius as a function of pixel coord
+	    // 0.5 PIX: get pixRadius 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;
@@ -150,10 +265,12 @@
 	    float r = hypot(x, y);
 
-	    psVectorAppend(radius, r);
+	    psVectorAppend(pixRadius, r);
 	    psVectorAppend(pixFlux, source->pixels->data.F32[iy][ix]);
 	    psVectorAppend(pixVar, source->variance->data.F32[iy][ix]);
-	}
-    }
-    psphotRadialAperturesSortFlux(radius, pixFlux, pixVar);
+	    nPix ++;
+	    // if (nPix % 10000 == 0) {fprintf (stderr, "?");}
+	}
+    }
+    psphotRadialAperturesSortFlux(pixRadius, pixFlux, pixVar);
 
     psVector *flux    = psVectorAllocEmpty(radMax->n, PS_TYPE_F32); // surface brightness of radial bin
@@ -174,6 +291,6 @@
 
     // 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) {
+    for (int i = 0; !done && (i < pixRadius->n); i++) {
+	if (pixRadius->data.F32[i] > Rmax) {
 	    // calculate the total flux for bin 'nOut'
 	    float Area = M_PI*PS_SQR(Rmax);
@@ -185,4 +302,7 @@
 		     nOut, radMax->data.F32[nOut], flux->data.F32[nOut], fluxErr->data.F32[nOut], fill->data.F32[nOut], Area);
 
+	    nPass ++;
+	    // if (nPass % 1000 == 0) {fprintf (stderr, "!");}
+
 	    nOut ++;
 	    if (nOut >= radMax->n) break;
@@ -195,16 +315,14 @@
     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);
+    radialAper->flux = flux;
+    radialAper->fluxErr = fluxErr;
+    radialAper->fill = fill;
+
+    psFree (pixRadius);
     psFree (pixFlux);
     psFree (pixVar);
 
+    nCalls ++;
+    // if (nCalls % 100 == 0) {fprintf (stderr, "*");}
     return true;
 }