Index: /branches/eam_branches/ipp-20110906/psphot/src/Makefile.am
===================================================================
--- /branches/eam_branches/ipp-20110906/psphot/src/Makefile.am	(revision 32615)
+++ /branches/eam_branches/ipp-20110906/psphot/src/Makefile.am	(revision 32616)
@@ -202,5 +202,4 @@
 	psphotRadialBins.c	       \
 	psphotRadialApertures.c	       \
-	psphotRadialAperturesByObject.c \
 	psphotPetrosian.c	       \
         psphotPetrosianRadialBins.c    \
Index: /branches/eam_branches/ipp-20110906/psphot/src/psphot.h
===================================================================
--- /branches/eam_branches/ipp-20110906/psphot/src/psphot.h	(revision 32615)
+++ /branches/eam_branches/ipp-20110906/psphot/src/psphot.h	(revision 32616)
@@ -428,5 +428,6 @@
 bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, int entry);
 bool psphotRadialApertures_Threaded (psThreadJob *job);
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry);
+bool psphotRadialApertureSource (pmSource *source, pmReadout *readout, int entry, psVector *pixRadius2, psVector *pixFlux, psVector *pixVer);
+// bool psphotRadialApertureSource (pmSource *source, int entry);
 
 bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
Index: /branches/eam_branches/ipp-20110906/psphot/src/psphotKronIterate.c
===================================================================
--- /branches/eam_branches/ipp-20110906/psphot/src/psphotKronIterate.c	(revision 32615)
+++ /branches/eam_branches/ipp-20110906/psphot/src/psphotKronIterate.c	(revision 32616)
@@ -192,5 +192,6 @@
     float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
 
-    for (int j = 0; j < 5; j++) {
+    // XXX TEST : set iteration to 1
+    for (int j = 0; j < 1; j++) {
 	for (int i = 0; i < sources->n; i++) {
 
@@ -209,13 +210,17 @@
 
 	    // use S/N to control max window size
-	    float kronSN = source->moments->KronFlux / source->moments->KronFluxErr;
+	    // float kronSN = source->moments->KronFlux / source->moments->KronFluxErr;
 
 	    // maxWindow -> 1.5*RADIUS for kronSN = 5.0, keeping S.B. constant (kronSN ~ flux)
 	    // (kronSN / maxWindow^2) = (5.0 / (1.5 RADIUS)^2)
 	    // maxWindow = 1.5 * RADIUS * sqrt(kronSN / 5.0)
-	    float maxWindow = (isfinite(kronSN) && (kronSN > 5.0)) ? 1.5 * RADIUS * sqrt(kronSN / 5.0) : 1.5*RADIUS;
+	    // XXX float maxWindow = (isfinite(kronSN) && (kronSN > 5.0)) ? 1.5 * RADIUS * sqrt(kronSN / 5.0) : 1.5*RADIUS;
 
 	    // iterate to the window radius
-	    float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow);
+	    // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow);
+
+	    // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS
+	    float maxWindow = source->skyRadius;
+	    float windowRadius = PS_MAX(RADIUS, maxWindow);
 
 	    // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
Index: /branches/eam_branches/ipp-20110906/psphot/src/psphotRadialApertures.c
===================================================================
--- /branches/eam_branches/ipp-20110906/psphot/src/psphotRadialApertures.c	(revision 32615)
+++ /branches/eam_branches/ipp-20110906/psphot/src/psphotRadialApertures.c	(revision 32616)
@@ -36,4 +36,11 @@
 }
 
+// these values are used by all threads repeatedly (and are not modified)
+static psVector *aperRadii = NULL;
+static psVector *aperRadii2 = NULL;
+static float outerRadius = NAN;
+static float SN_LIM = NAN;
+static psImageMaskType maskVal = 0;
+
 // aperture-like measurements for extended sources
 // flux in simple, circular apertures
@@ -75,4 +82,25 @@
         nThreads = 0;
     }
+
+    // aperRadii stores the upper bounds of the annuli
+    // XXX keep the same name here as for the petrosian / elliptical apertures?
+    aperRadii = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    psAssert (aperRadii, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
+    psAssert (aperRadii->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
+
+    outerRadius = aperRadii->data.F32[aperRadii->n - 1];
+
+    // save the R^2 values as well for quicker comparison
+    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]);
+    }
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
+
+    // S/N limit to perform full non-linear fits
+    SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
     // source analysis is done in S/N order (brightest first)
@@ -117,12 +145,9 @@
             psArrayAdd(job->args, 1, cells->data[j]); // sources
             psArrayAdd(job->args, 1, AnalysisRegion);
-            psArrayAdd(job->args, 1, recipe);
-
             PS_ARRAY_ADD_SCALAR(job->args, entry,  PS_TYPE_S32);
             PS_ARRAY_ADD_SCALAR(job->args, nEntry, PS_TYPE_S32);
-
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial
 
-// set this to 0 to run without threading
+	    // set this to 0 to run without threading
 # if (1)	    
             if (!psThreadJobAddPending(job)) {
@@ -138,5 +163,5 @@
 	    }
 	    psScalar *scalar = NULL;
-	    scalar = job->args->data[6];
+	    scalar = job->args->data[5];
 	    Nradial += scalar->data.S32;
 	    psFree(job);
@@ -158,5 +183,5 @@
             } else {
 		psScalar *scalar = NULL;
-		scalar = job->args->data[6];
+		scalar = job->args->data[5];
 		Nradial += scalar->data.S32;
             }
@@ -166,4 +191,5 @@
     psFree (cellGroups);
     psFree(AnalysisRegion);
+    psFree (aperRadii2);
 
     psLogMsg ("psphot", PS_LOG_WARN, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
@@ -173,5 +199,4 @@
 bool psphotRadialApertures_Threaded (psThreadJob *job) {
 
-    bool status;
     int Nradial = 0;
 
@@ -180,23 +205,11 @@
     psArray *sources        = job->args->data[1];
     psRegion *region        = job->args->data[2];
-    psMetadata *recipe      = job->args->data[3];
-
-    int entry               = PS_SCALAR_VALUE(job->args->data[4],S32); // which psf-matched image are we working on? (0 == unmatched)
-    int nEntry              = PS_SCALAR_VALUE(job->args->data[5],S32); // total number of psf-matched images + 1 unmatched
-
-    // 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);
-
-    // S/N limit to perform full non-linear fits
-    float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
-
-    float outerRadius = radMax->data.F32[radMax->n - 1];
+    int entry               = PS_SCALAR_VALUE(job->args->data[3],S32); // which psf-matched image are we working on? (0 == unmatched)
+    int nEntry              = PS_SCALAR_VALUE(job->args->data[4],S32); // total number of psf-matched images + 1 unmatched
+
+    // storage for the derived pixel values (these are passed into psphotRadialApertureSource)
+    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
 
     // choose the sources of interest
@@ -240,17 +253,7 @@
 	}
 
-	// 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...
-	pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
-
-	if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, entry)) {
+	if (!psphotRadialApertureSource (source, readout, entry, pixRadius2, pixFlux, pixVar)) {
 	    psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
 	} else {
@@ -258,20 +261,19 @@
 	}
 
-	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);
     }
-    psScalar *scalar = job->args->data[6];
+    psScalar *scalar = job->args->data[5];
     scalar->data.S32 = Nradial;
+
+    psFree (pixRadius2);
+    psFree (pixFlux);
+    psFree (pixVar);
 
     return true;
 }
 
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
-
+bool psphotRadialApertureSource (pmSource *source, pmReadout *readout, int entry, psVector *pixRadius2, psVector *pixFlux, psVector *pixVar) {
+					    
     // if we are a child source, save the results to the parent source radial aperture array
     psArray *radialAperSet = source->radialAper;
@@ -285,59 +287,68 @@
     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);
+    // find the largest aperture of interest (use only apertures with inner radii <=
+    // source->skyRadius)
+    int lastAp = aperRadii->n;
+    for (int i = 0; i < aperRadii->n; i++) {
+	if (aperRadii->data.F32[i] < source->skyRadius) continue;
+	lastAp = i + 1;
+	break;
+    }
 
     // outer-most radius for initial truncation
-    float Rmax  = aperRadii->data.F32[aperRadii->n - 1];
+    float Rmax  = aperRadii->data.F32[lastAp - 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]);
-    }
+    // in this function, the operatins are relative to the full image (readout->image, etc)
 
     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
+	xCM = source->moments->Mx; // index coord of peak in readout
+	yCM = source->moments->My; // index coord of peak in readout
     } 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
-    }
+	xCM = source->peak->xf; // index coord of peak in readout
+	yCM = source->peak->yf; // index coord of peak in readout
+    }
+
+    int Nx = readout->image->numCols;
+    int Ny = readout->image->numRows;
+
+    pixRadius2->n = 0;
+    pixFlux->n = 0;
+    pixVar->n = 0;
 
     // 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;
+    for (int iy = -Rmax; iy < Rmax + 1; iy++) {
+
+	float yDiff = iy + 0.5 + yCM;  // y-coordinate at this offse
+	int yPix = (int) yDiff;
+
+	if (yPix < 0) continue;
+	if (yPix > Ny - 1) continue;
+	if (fabs(iy) > Rmax) continue;
+
+	float *vPix = readout->image->data.F32[yPix];
+	float *vWgt = readout->variance->data.F32[yPix];
+	psImageMaskType  *vMsk = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[yPix];
+
+	for (int ix = -Rmax; ix < Rmax + 1; ix++) {
+
+	    float xDiff = ix + 0.5 + xCM;  // x-coordinate at this offse
+	    int xPix = (int) xDiff;
+	    
+	    if (xPix < 0) continue;
+	    if (xPix > Nx - 1) continue;
+	    if (fabs(ix) > Rmax) continue;
+	    
+	    if (vMsk[xPix] & maskVal) continue;
+	    if (isnan(vPix[xPix])) continue;
 
 	    // radius is just a function of (xDiff, yDiff)
-	    float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+	    float r2  = PS_SQR(ix) + PS_SQR(iy);
 	    if (r2 > Rmax2) continue;
 
 	    psVectorAppend(pixRadius2, r2);
-	    psVectorAppend(pixFlux, *vPix);
-	    psVectorAppend(pixVar, *vWgt);
+	    psVectorAppend(pixFlux, vPix[xPix]);
+	    psVectorAppend(pixVar, vWgt[xPix]);
 	}
     }
@@ -348,8 +359,15 @@
     psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
 
+    // init the apertures of interest to 0.0, the rest go to NAN
     psVectorInit (flux,    0.0);
     psVectorInit (fluxStd, 0.0);
     psVectorInit (fluxErr, 0.0);
     psVectorInit (fill,    0.0);
+    for (int i = lastAp; i < flux->n; i++) {
+	flux->data.F32[i] = NAN;
+	fluxStd->data.F32[i] = NAN;
+	fluxErr->data.F32[i] = NAN;
+	fill->data.F32[i] = NAN;
+    }
 
     float *rPix2 = pixRadius2->data.F32;
@@ -358,6 +376,8 @@
 	int j = 0;
 	float *aRad2 = aperRadii2->data.F32;
-	for (; (*aRad2 < *rPix2) && (j < aperRadii2->n); j++, aRad2++);
-	for (; j < aperRadii2->n; j++, aRad2++) {
+	for (; (*aRad2 < *rPix2) && (j < lastAp); j++, aRad2++);
+
+	// XXX I can speed this up by only saving this single aperture
+	for (; j < lastAp; j++, aRad2++) {
 	    flux->data.F32[j]    += pixFlux->data.F32[i];
 	    fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]);
@@ -371,7 +391,7 @@
        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++) {
+    */
+
+    for (int i = 0; i < lastAp; i++) {
 	// calculate the total flux for bin 'nOut'
 	float Area = M_PI*aperRadii2->data.F32[i];
@@ -383,5 +403,5 @@
 	// XXX report the total flux or the mask-corrected flux?
 	// flux->data.F32[i]    = SBmean * Area;
-	// fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
+	// fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / Pinx;
 
 	fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
@@ -393,16 +413,21 @@
     }
     
+# if (1)
     radialAper->flux = flux;
     radialAper->fluxStdev = fluxStd;
     radialAper->fluxErr = fluxErr;
     radialAper->fill = fill;
-
-    psFree (aperRadii2);
-    psFree (pixRadius2);
-    psFree (pixFlux);
-    psFree (pixVar);
+# else
+    // XXX TEST
+    psFree(flux);
+    psFree(fluxStd);
+    psFree(fluxErr);
+    psFree(fill);
+# endif
 
     return true;
 }
+
+/*** below is a test to use a sort to speed this up, not very successfully ***/
 
 static int nCalls = 0;
@@ -495,18 +520,18 @@
 // *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
 # define COMPARE_VECT(A,B) (radius->data.F32[A] < radius->data.F32[B])
-# define SWAP_VECT(TYPE,A,B) { \
-  float tmp; \
-  if (A != B) { \
-    tmp = radius->data.F32[A]; \
-    radius->data.F32[A] = radius->data.F32[B];	\
-    radius->data.F32[B] = tmp; \
-    tmp = pixFlux->data.F32[A]; \
-    pixFlux->data.F32[A] = pixFlux->data.F32[B]; \
-    pixFlux->data.F32[B] = tmp; \
-    tmp = pixVar->data.F32[A]; \
-    pixVar->data.F32[A] = pixVar->data.F32[B]; \
-    pixVar->data.F32[B] = tmp; \
-  } \
-}
+# define SWAP_VECT(TYPE,A,B) {					\
+	float tmp;						\
+	if (A != B) {						\
+	    tmp = radius->data.F32[A];				\
+	    radius->data.F32[A] = radius->data.F32[B];		\
+	    radius->data.F32[B] = tmp;				\
+	    tmp = pixFlux->data.F32[A];				\
+	    pixFlux->data.F32[A] = pixFlux->data.F32[B];	\
+	    pixFlux->data.F32[B] = tmp;				\
+	    tmp = pixVar->data.F32[A];				\
+	    pixVar->data.F32[A] = pixVar->data.F32[B];		\
+	    pixVar->data.F32[B] = tmp;				\
+	}							\
+    }
 
 bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar) {
Index: /branches/eam_branches/ipp-20110906/psphot/src/psphotRadialProfileWings.c
===================================================================
--- /branches/eam_branches/ipp-20110906/psphot/src/psphotRadialProfileWings.c	(revision 32615)
+++ /branches/eam_branches/ipp-20110906/psphot/src/psphotRadialProfileWings.c	(revision 32616)
@@ -132,5 +132,5 @@
 
 // set this to 0 to run without threading
-# if (0)
+# if (1)
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
@@ -217,4 +217,9 @@
 
     // psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+
+    // XXX unclear if I should run this analysis on both 1st and 2nd pass for 1st pass objects,
+    // or just use the 1st pass value.  I think I should just use the 1st pass value, so skip
+    // any that have already been assigned
+    if (isfinite(source->skyRadius)) return true;
 
     // radii will be MIN_RADIUS to MAX_RADIUS in NN log steps:
Index: /branches/eam_branches/ipp-20110906/psphot/src/psphotReadout.c
===================================================================
--- /branches/eam_branches/ipp-20110906/psphot/src/psphotReadout.c	(revision 32615)
+++ /branches/eam_branches/ipp-20110906/psphot/src/psphotReadout.c	(revision 32616)
@@ -308,4 +308,7 @@
 pass1finish:
 
+    // measure the radial profiles to the sky (only measures new objects)
+    psphotRadialProfileWings (config, view, filerule);
+
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, filerule);
Index: /branches/eam_branches/ipp-20110906/psphot/src/psphotSetThreads.c
===================================================================
--- /branches/eam_branches/ipp-20110906/psphot/src/psphotSetThreads.c	(revision 32615)
+++ /branches/eam_branches/ipp-20110906/psphot/src/psphotSetThreads.c	(revision 32616)
@@ -50,5 +50,5 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 7);
+    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 6);
     task->function = &psphotRadialApertures_Threaded;
     psThreadTaskAdd(task);
Index: /branches/eam_branches/ipp-20110906/psphot/src/psphotSourceSize.c
===================================================================
--- /branches/eam_branches/ipp-20110906/psphot/src/psphotSourceSize.c	(revision 32615)
+++ /branches/eam_branches/ipp-20110906/psphot/src/psphotSourceSize.c	(revision 32616)
@@ -156,5 +156,5 @@
 
     // XXX fix this (was source->n  - first)
-    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size"));
+    psLogMsg ("psphot.size", PS_LOG_WARN, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size"));
 
     psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
Index: /branches/eam_branches/ipp-20110906/psphot/src/psphotStackReadout.c
===================================================================
--- /branches/eam_branches/ipp-20110906/psphot/src/psphotStackReadout.c	(revision 32615)
+++ /branches/eam_branches/ipp-20110906/psphot/src/psphotStackReadout.c	(revision 32616)
@@ -168,4 +168,7 @@
     psphotFitSourcesLinear (config, view, STACK_SRC, false);
     psphotStackVisualFilerule(config, view, STACK_SRC);
+
+    // measure the radial profiles to the sky
+    psphotRadialProfileWings (config, view, STACK_SRC);
 
     // re-measure the kron mags with models subtracted.  this pass starts with a circular
@@ -245,4 +248,7 @@
 pass1finish:
 
+    // measure the radial profiles to the sky (only measures new objects)
+    psphotRadialProfileWings (config, view, STACK_SRC);
+
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, STACK_SRC);
@@ -254,4 +260,7 @@
 
     psMemDump("psfstats");
+
+    // XXX *** NOTE : this should happen above before the last psphotFitSourcesLinear in order
+    // to generate psf forced photometry for the undetected sources
 
     // generate the objects (objects unify the sources from the different images)
@@ -281,7 +290,5 @@
 
     // measure circular, radial apertures (objects sorted by S/N)
-    // XXX can we just use psphotRadialApertures
-    // XXX make sure the headers are consistent with this (which PSF convolutions, ie mark 'none')
-    // psphotRadialAperturesByObject (config, objectsRadial, view, STACK_SRC, nMatchedPSF); 
+    // this forces photometry on the undetected sources from other images
     psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
     psMemDump("extmeas");