Index: /branches/eam_branches/ipp-20101205/psphot/src/Makefile.am
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/Makefile.am	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/Makefile.am	(revision 30101)
@@ -98,4 +98,5 @@
 	psphotStackMatchPSFsUtils.c   \
 	psphotStackMatchPSFsPrepare.c \
+	psphotStackMatchPSFsNext.c    \
 	psphotStackOptions.c          \
 	psphotStackObjects.c          \
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphot.h
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphot.h	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphot.h	(revision 30101)
@@ -312,4 +312,5 @@
 
 int psphotFileruleCount(const pmConfig *config, const char *filerule);
+bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num);
 
 bool psphotAddKnownSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *inSources);
@@ -403,4 +404,6 @@
 bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
 bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
+bool psphotStackMatchPSFsNext (bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int nextSize);
+bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, float currentFWHM, float targetFWHM);
 
 // psphotStackMatchPSFsUtils
@@ -426,8 +429,8 @@
 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);
+bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry);
 
 bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
-bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);
+bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule, int nMatchedPSF);
 
 bool psphotStackObjectsUnifyPosition (psArray *objects);
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotExtendedSourceAnalysisByObject.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphotExtendedSourceAnalysisByObject.c	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphotExtendedSourceAnalysisByObject.c	(revision 30101)
@@ -53,4 +53,7 @@
 	pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
 	psAssert (readout, "missing readout?");
+
+	psLogMsg("psphot", PS_LOG_INFO, "petrosians for image %d", i);
+	psphotVisualShowImage(readout);
 
 	readouts->data[i] = psMemIncrRefCounter(readout);
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotExtendedSourceFits.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphotExtendedSourceFits.c	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphotExtendedSourceFits.c	(revision 30101)
@@ -18,6 +18,11 @@
     int num = psphotFileruleCount(config, filerule);
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+	if (i == chisqNum) continue; // skip chisq image
 	if (!psphotExtendedSourceFitsReadout (config, view, filerule, i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on to fit extended sources for %s entry %d", filerule, i);
@@ -37,4 +42,5 @@
     int Nplain = 0;
     int NplainPass = 0;
+    int Nfaint = 0;
 
     psTimerStart ("psphot.extended");
@@ -46,4 +52,7 @@
     pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     psAssert (readout, "missing readout?");
+
+    psLogMsg("psphot", PS_LOG_INFO, "extended source fits for image %d", index);
+    psphotVisualShowImage(readout);
 
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
@@ -167,4 +176,5 @@
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nplain
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NplainPass
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfain
 
             if (false && !psThreadJobAddPending(job)) {
@@ -189,4 +199,6 @@
                 scalar = job->args->data[11];
                 NplainPass += scalar->data.S32;
+                scalar = job->args->data[12];
+                Nfaint += scalar->data.S32;
 		psFree(job);
 	    }
@@ -217,4 +229,6 @@
                 scalar = job->args->data[11];
                 NplainPass += scalar->data.S32;
+                scalar = job->args->data[12];
+                Nfaint += scalar->data.S32;
             }
             psFree(job);
@@ -227,4 +241,5 @@
     psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit\n", Nfaint);
     return true;
 }
@@ -238,4 +253,5 @@
     int NconvolvePass = 0;
     int Nplain = 0;
+    int Nfaint = 0;
     int NplainPass = 0;
     bool savePics = false;
@@ -350,5 +366,8 @@
           // limit selection to some SN limit
           assert (source->peak); // how can a source not have a peak?
-          if (source->peak->SN < SNlim) continue;
+          if (source->peak->SN < SNlim) {
+	      Nfaint ++;
+	      continue;
+	  }
 
           // check on the model type
@@ -497,4 +516,7 @@
     scalar->data.S32 = NplainPass;
 
+    scalar = job->args->data[12];
+    scalar->data.S32 = Nfaint;
+
     return true;
 }
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotRadialApertures.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphotRadialApertures.c	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphotRadialApertures.c	(revision 30101)
@@ -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,116 @@
 }
 
-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) {
+
+    psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?");
+
+    pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc ();
+    source->radialAper->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 +253,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 +259,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 +285,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 +296,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 +309,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;
 }
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotRadialAperturesByObject.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphotRadialAperturesByObject.c	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphotRadialAperturesByObject.c	(revision 30101)
@@ -3,5 +3,14 @@
 // aperture-like measurements for extended sources
 // flux in simple, circular apertures
-bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule) {
+
+// **** it looks like this function will re-point the source pixels at the specified FILERULE
+// **** I need to distinguish PSF-matched images from raw
+// **** save (somewhere) the PSF-matched PSF values
+
+// this function measures the radial aperture fluxes for the set of readouts.  this function
+// may be called multiple times (presumably with different matched PSF sizes).  we must have
+// already added an entry to the readout->analysis identifying the FWHM of this version.
+
+bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule, int nMatchedPSF) {
 
     bool status;
@@ -28,4 +37,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)
@@ -39,4 +49,17 @@
     float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
+    // how many target PSFs do we want?
+    int nPSFsizes = 0;
+    { 
+	psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
+	if (status) {
+	    nPSFsizes = 1;
+	} else {
+	    psVector *fwhmValues = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
+	    psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+	    nPSFsizes = fwhmValues->n;
+	}
+    }
+    
     // source analysis is done in S/N order (brightest first)
     objects = psArraySort (objects, pmPhotObjSortBySN);
@@ -53,9 +76,21 @@
 	psAssert (readout, "missing readout?");
 
+	psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+	if (!fwhmValues) {
+	    psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout");
+	    return false;
+	}
+	if (fwhmValues->n != nMatchedPSF + 1) {
+	    psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM sequence is inconsistent this readout");
+	    return false;
+	}
+	psLogMsg ("psphot", PS_LOG_DETAIL, "PSF FWHM of %s : %f pixels\n", file->name, fwhmValues->data.F32[nMatchedPSF]);
+
 	readouts->data[i] = psMemIncrRefCounter(readout);
     }
 
     // process the objects in order.  
-    for (int i = 0; i < objects->n; i++) {
+    // XXX TEST: fix this!
+    for (int i = 0; (i < 50) && (i < objects->n); i++) {
         pmPhotObj *object = objects->data[i];
 	if (!object) continue;
@@ -77,4 +112,16 @@
 	    if (source->peak->SN < SN_LIM) continue;
 
+	    int index = source->imageID;
+	    if (index >= readouts->n) continue; // skip the sources generated by the chisq image
+	    pmReadout *readout = readouts->data[index];
+
+	    // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
+	    // psphotVisualShowImage(readout);
+
+	    // allocate pmSourceExtendedParameters, if not already defined
+	    if (!source->radialAper) {
+		source->radialAper = psArrayAlloc(nPSFsizes);
+	    }
+
 	    // replace object in image
 	    if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
@@ -83,15 +130,15 @@
 	    Nradial ++;
 
-	    int index = source->imageID;
-	    pmReadout *readout = readouts->data[index];
+	    // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
+	    // psphotVisualShowImage(readout);
 
 	    // 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);
+	    // 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, outerRadius + 2);
 
-	    if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) {
+	    if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, nMatchedPSF)) {
 		psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
 	    } else {
@@ -101,4 +148,7 @@
 	    // re-subtract the object, leave local sky
 	    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+
+	    // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
+	    // psphotVisualShowImage(readout);
 	}
     }
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotSetThreads.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphotSetThreads.c	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphotSetThreads.c	(revision 30101)
@@ -35,5 +35,5 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 12);
+    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 13);
     task->function = &psphotExtendedSourceFits_Threaded;
     psThreadTaskAdd(task);
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotStackChisqImage.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphotStackChisqImage.c	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphotStackChisqImage.c	(revision 30101)
@@ -6,5 +6,5 @@
 
 // XXX supply filename or keep PSPHOT.INPUT fixed?
-bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleCnv)
+bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleSrc)
 {
     psTimerStart ("psphot.chisq.image");
@@ -27,6 +27,7 @@
 
     psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.CHISQ.NUM", PS_META_REPLACE, "", num);
+
+    // we need to increment the counter for ruleDet and ruleSrc:
     num++;
-    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", num);
 
     // save the resulting image in the 'detection' set
@@ -35,10 +36,13 @@
         return false;
     }
+    psphotFileruleCountSet(config, ruleDet, num);
 
-    // save the resulting image in the 'convolved' set
-    // XXX whil am I doing this as well?
-    if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleCnv, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) {
-        psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files");
-        return false;
+    // also save the resulting image in the 'source' set (analysis set)
+    if (strcmp(ruleDet, ruleSrc)) {
+	if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleSrc, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) {
+	    psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files");
+	    return false;
+	}
+	psphotFileruleCountSet(config, ruleSrc, num);
     }
 
@@ -120,10 +124,10 @@
     psAssert (status, "programming error: must define PSPHOT.CHISQ.NUM");
 
-    int inputNum = psphotFileruleCount(config, "PSPHOT.INPUT");
+    int inputNum = psphotFileruleCount(config, filerule);
 
     pmFPAfileRemoveSingle (config->files, filerule, chisqNum);
 
     inputNum --;
-    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", inputNum);
+    psphotFileruleCountSet(config, filerule, inputNum);
 
     return true;
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFs.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFs.c	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFs.c	(revision 30101)
@@ -106,5 +106,14 @@
     rescaleData(readoutOut, config, options, index);
 
-    // dumpImage(readoutOut, readoutSrc, index, "convolved");
+    // save the output fwhm values in the readout->analysis.  we may have / will have multiple output PSF sizes,
+    // so we save this in a vector.  if the vector is not yet defined, create it
+    bool mdok = false;
+    psVector *fwhmValues = psMetadataLookupVector(&mdok, readoutOut->analysis, "STACK.PSF.FWHM.VALUES");
+    if (!fwhmValues) {
+	fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
+	psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
+	psFree(fwhmValues); // drops the extra copy
+    }
+    psVectorAppend(fwhmValues, options->targetSeeing);
 
     return true;
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFsNext.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFsNext.c	(revision 30101)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphotStackMatchPSFsNext.c	(revision 30101)
@@ -0,0 +1,138 @@
+# include "psphotInternal.h"
+
+// smooth the input image to match the next target PSF
+// this function assumes the image has already been smoothed to match the first value (array element 0),
+// and that the smoothing can use a 1D Gaussian kernel of width sqrt(TARGET^2 - CURRENT^2)
+// each subsequent call
+bool psphotStackMatchPSFsNext(bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    psVector *fwhmValues = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
+    if (!status) {
+	// must not be a vector, only one value requested
+	*smoothAgain = false;
+	return true;
+    }
+
+    if (lastSize + 1 >= fwhmValues->n) {
+	// all done with target FWHM values 
+	*smoothAgain = false;
+	return true;
+    }
+
+    float targetFWHM = fwhmValues->data.F32[lastSize + 1];
+    float currentFWHM = fwhmValues->data.F32[lastSize];
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // smooth the image and variance map
+    bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading in psImageConvolve
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+	if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, recipe, currentFWHM, targetFWHM)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF %f", filerule, i, targetFWHM);
+	    psImageConvolveSetThreads(oldThreads);
+	    return false;
+	}
+    }
+
+    psImageConvolveSetThreads(oldThreads);
+    return true;
+}
+
+bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, float currentFWHM, float targetFWHM) {
+
+    bool status = false;
+
+    // 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?");
+
+    psLogMsg("psphot", PS_LOG_INFO, "smooth %d to next psf", index);
+    psphotVisualShowImage(readout);
+
+    // 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
+    psAssert (maskVal, "missing mask value?");
+
+    float minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
+    if (!isfinite(minGauss)) {
+        psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value");
+        minGauss = 0.5;
+    }
+
+    if (targetFWHM <= currentFWHM) {
+	psError (PSPHOT_ERR_CONFIG, true, "target FWHM cannot be smaller than current FWHM");
+	return false;
+    }
+
+    float smoothFWHM  = sqrt(PS_SQR(targetFWHM) - PS_SQR(currentFWHM));
+    float SIGMA_SMTH  = smoothFWHM / (2.0*sqrt(2.0*log(2.0)));
+    float NSIGMA_SMTH = 3.0;
+    
+    // record the actual smoothing sigma
+    // psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGMA_SMOOTH", PS_META_REPLACE, "Smoothing sigma for detections", SIGMA_SMTH);
+
+    // smooth the image, applying the mask as we go
+    psImageSmoothMask_Threaded(readout->image, readout->image, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH, minGauss);
+    psLogMsg("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark("psphot.smooth"));
+
+    // Smooth the variance, applying the mask as we go.  The variance is smoothed by the PSF^2,
+    // renomalized to maintain the input level of the variance.  We achieve this by smoothing
+    // with a Gaussian with sigma = SIGMA_SMTH/sqrt(2) with unity normalization.  Note that
+    // this process yields a smoothed image with correlated errors.  The pixel-to-pixel
+    // variations in smooth_im will be decreased by a factor of 4*pi*SIGMA_SMTH^2, but for
+    // measurements based on apertures comparable to or larger than the smoothing kernel, the
+    // effective per-pixel variance is maintained.
+    psImageSmoothMask_Threaded(readout->variance, readout->variance, readout->mask, maskVal, SIGMA_SMTH * M_SQRT1_2, NSIGMA_SMTH, minGauss);
+    psLogMsg("psphot", PS_LOG_MINUTIA, "smooth variance: %f sec\n", psTimerMark("psphot.smooth"));
+
+    psLogMsg("psphot", PS_LOG_INFO, "smoothed");
+    psphotVisualShowImage(readout);
+
+    // optionally save example images under trace
+    if (psTraceGetLevel("psphot") > 5) {
+	static int pass = 0;
+        char name[64];
+        sprintf (name, "stksm.v%d.fits", pass);
+        psphotSaveImage(NULL, readout->image, name);
+        sprintf (name, "stkwt.v%d.fits", pass);
+        psphotSaveImage(NULL, readout->variance, name);
+	pass ++;
+    }
+
+    // XXX need to apply this to the radial apertures somehow.
+    // Calculate correction factor for the covariance produced by the (potentially multiple) smoothing
+    // psKernel *kernel = psImageSmoothKernel(SIGMA_SMTH, NSIGMA_SMTH); // Kernel used for smoothing
+    // double sum2 = 0.0;                                               // Sum of kernel squared
+    // for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+    //     for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+    //         sum2 += PS_SQR(kernel->kernel[y][x]);
+    //     }
+    // }
+    // float factor = 1.0 / (sum2 * psImageCovarianceCalculateFactor(kernel, readout->covariance));
+    // psFree(kernel);
+
+    // record the effective area and significance scaling factor
+    float effArea = 8.0 * M_PI * PS_SQR(SIGMA_SMTH);
+    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea);
+    // psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor);
+
+    // save the output fwhm values in the readout->analysis.  we may have / will have multiple output PSF sizes,
+    // so we save this in a vector.  if the vector is not yet defined, create it
+
+    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+    psAssert(fwhmValues, "should already exist..");
+    psVectorAppend(fwhmValues, targetFWHM);
+
+    return true;
+}
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotStackPSF.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphotStackPSF.c	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphotStackPSF.c	(revision 30101)
@@ -12,12 +12,13 @@
     bool autoPSF = psMetadataLookupBool (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.AUTO");
 
+    // Get the recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+	
+    char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF
+
     if (autoPSF) {
-	// Get the recipe values
-	psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
-	psAssert(recipe, "We've thrown an error on this before.");
-
 	int psfInstances = psMetadataLookupS32(NULL, recipe, "PSF.INSTANCES"); // Number of instances for PSF
 	float psfRadius = psMetadataLookupF32(NULL, recipe, "PSF.RADIUS"); // Radius for PSF
-	const char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF
 	int psfOrder = psMetadataLookupS32(NULL, recipe, "PSF.ORDER"); // Spatial order for PSF
 
@@ -49,10 +50,14 @@
 
 	float targetFWHM = psMetadataLookupF32 (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
-	psAssert (isfinite(targetFWHM), "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+	if (!mdok) {
+	    psVector *fwhmValues = psMetadataLookupVector(&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
+	    psAssert (mdok, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+	    targetFWHM = fwhmValues->data.F32[0];
+	}
 
 	float Sxx = sqrt(2.0)*targetFWHM / 2.35;
 
 	// XXX probably should make the model type (and par 7) optional from recipe
-	psf = pmPSFBuildSimple("PS_MODEL_PS1_V1", Sxx, Sxx, 0.0, 1.0);
+	psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0);
 	if (!psf) {
 	    psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF.");
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotStackParseCamera.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphotStackParseCamera.c	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphotStackParseCamera.c	(revision 30101)
@@ -15,4 +15,6 @@
     }
 
+    int nRaw = 0;
+    int nCnv = 0;
     int nInputs = inputs->list->n;
     for (int i = 0; i < nInputs; i++) {
@@ -57,4 +59,5 @@
 		}
 	    }
+	    nRaw ++;
 	}
 
@@ -88,4 +91,5 @@
 		}
 	    }
+	    nCnv ++;
 	}
 
@@ -93,4 +97,12 @@
 	    psError(PSPHOT_ERR_CONFIG, true, "Component %s (%d) lacks both RAW:IMAGE and CNV:IMAGE of type STR", item->name, i);
 	    return false;
+	}
+
+	// XXX what if they do not match in length
+	if (nCnv && nRaw) {
+	    if (nCnv != nRaw) {
+		psError (PSPHOT_ERR_CONFIG, true, "if both RAW and CNV images are supplied, the number must match");
+		return false;
+	    }
 	}
 
@@ -150,4 +162,7 @@
     }
     psMetadataRemoveKey(config->arguments, "FILENAMES");
+    psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.INPUT.RAW.NUM", PS_META_REPLACE, "number of inputs", nRaw);
+    psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.INPUT.CNV.NUM", PS_META_REPLACE, "number of inputs", nCnv);
+    psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.OUTPUT.IMAGE.NUM", PS_META_REPLACE, "number of inputs", nInputs);
     psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", nInputs);
 
Index: /branches/eam_branches/ipp-20101205/psphot/src/psphotStackReadout.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psphot/src/psphotStackReadout.c	(revision 30100)
+++ /branches/eam_branches/ipp-20101205/psphot/src/psphotStackReadout.c	(revision 30101)
@@ -1,7 +1,13 @@
 # include "psphotInternal.h"
 
+// we have 3 possible real filesets:
 # define STACK_RAW "PSPHOT.STACK.INPUT.RAW"
 # define STACK_CNV "PSPHOT.STACK.INPUT.CNV"
-# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"
+# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"  /* the psf-matched image */
+
+// we have 3 files on which we operate:
+// DET (detection image)       : nominally RAW (optionally CNV?)
+// SRC (source analysis image) : nominally CNV (optionally RAW)
+// OUT (psf-matched images)    : always OUT
 
 bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
@@ -21,10 +27,13 @@
     PS_ASSERT_PTR_NON_NULL (breakPt, false);
 
-    BST = RAW : CNV;
+    // XXX or do I set OUT to be a pmFPAfile pointing at the input of interest?
+    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
+    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
+    char *STACK_DET = STACK_RAW; // XXX optionally allow this to be CNV?
 
     // we have 3 relevant files: RAW, CNV, OUT 
 
     // set the photcode for each image
-    if (!psphotAddPhotcode (config, view, STACK_BST)) {
+    if (!psphotAddPhotcode (config, view, STACK_SRC)) {
         psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
         return false;
@@ -33,9 +42,9 @@
     // Generate the mask and weight images
     // XXX this should be done before we perform the convolutions
-    if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) {
-	return psphotReadoutCleanup (config, view, STACK_OUT);
+    if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
+	return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "NOTHING")) {
-	return psphotReadoutCleanup (config, view, STACK_OUT);
+	return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
@@ -43,57 +52,53 @@
     // XXX I think this is not defined correctly for an array of images.
     // XXX probably need to subtract the model (same model?) for both RAW and OUT
-    if (!psphotModelBackground (config, view, STACK_RAW)) {
-	return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
-	return psphotReadoutCleanup (config, view, STACK_OUT);
+    if (!psphotModelBackground (config, view, STACK_DET)) {
+	return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_DET)) {
+	return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "BACKMDL")) {
-	return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-
-    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are determined and saved on
-    // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT
-    if (!psphotLoadPSF (config, view, STACK_RAW)) {
-        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-
-    if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_OUT)) {
+	return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "CHISQ")) {
-	return psphotReadoutCleanup (config, view, STACK_OUT);
+	return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // find the detections (by peak and/or footprint) in the image.
     // This finds the detections on Chisq image as well as the individuals
-    if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1
+    if (!psphotFindDetections (config, view, STACK_DET, true)) { // pass 1
         // this only happens if we had an error in psphotFindDetections
         psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
-    }
-
-    // copy the detections from RAW to OUT
-    if (!psphotCopySources (config, view, STACK_OUT, STACK_RAW)) {
-        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    // copy the detections from DET to SRC 
+    if (strcmp(STACK_SRC, STACK_DET)) {
+	if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
+	    psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
+	    return psphotReadoutCleanup (config, view, STACK_SRC);
+	}
     }
 
     // construct sources and measure basic stats (saved on detections->newSources)
-    // only run this on detections from the input images, not chisq image
-    if (!psphotSourceStats (config, view, STACK_OUT, true)) { // pass 1
+    if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // generate the objects (object unify the sources from the different images)
-    psArray *objects = psphotMatchSources (config, view, STACK_OUT);
+    // XXX this could just match the detections for the chisq image, and not bother measuring the
+    // source stats in that case...
+    psArray *objects = psphotMatchSources (config, view, STACK_SRC);
 
     // construct sources for the newly-generated sources (from other images)
-    if (!psphotSourceStats (config, view, STACK_OUT, false)) { // pass 1
+    if (!psphotSourceStats (config, view, STACK_SRC, false)) { // pass 1
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
@@ -101,5 +106,5 @@
     // if (!psphotDeblendSatstars (config, view)) {
     //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-    //     return psphotReadoutCleanup (config, view, STACK_OUT);
+    //     return psphotReadoutCleanup (config, view, STACK_SRC);
     // }
 
@@ -107,39 +112,39 @@
     // if (!psphotBasicDeblend (config, view)) {
     //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
-    //     return psphotReadoutCleanup (config, view, STACK_OUT);
+    //     return psphotReadoutCleanup (config, view, STACK_SRC);
     // }
 
     // classify sources based on moments, brightness
     // only run this on detections from the input images, not chisq image
-    if (!psphotRoughClass (config, view, STACK_OUT)) {
+    if (!psphotRoughClass (config, view, STACK_SRC)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
-	return psphotReadoutCleanup (config, view, STACK_OUT);
+	return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     // only run this on detections from the input images, not chisq image
-    if (!psphotImageQuality (config, view, STACK_OUT)) { // pass 1
+    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "MOMENTS")) {
-	return psphotReadoutCleanup (config, view, STACK_OUT);
+	return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
     // this step is skipped
-    if (!psphotChoosePSF (config, view, STACK_OUT)) { // pass 1
+    if (!psphotChoosePSF (config, view, STACK_SRC)) { // pass 1
         psLogMsg ("psphot", 3, "failure to construct a psf model");
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
     if (!strcasecmp (breakPt, "PSFMODEL")) {
-        return psphotReadoutCleanup (config, view, STACK_OUT);
+        return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
-    psphotGuessModels (config, view, STACK_OUT);
+    psphotGuessModels (config, view, STACK_SRC);
 
     // merge the newly selected sources into the existing list
     // NOTE: merge OLD and NEW
-    psphotMergeSources (config, view, STACK_OUT);
+    psphotMergeSources (config, view, STACK_SRC);
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
@@ -147,26 +152,37 @@
 
     // identify CRs and extended sources
-    psphotSourceSize (config, view, STACK_OUT, TRUE);
+    psphotSourceSize (config, view, STACK_SRC, TRUE);
 
     // measure aperture photometry corrections
-    if (!psphotApResid (config, view, STACK_OUT)) {
+    if (!psphotApResid (config, view, STACK_SRC)) {
 	psFree (objects);
         psLogMsg ("psphot", 3, "failed on psphotApResid");
-	return psphotReadoutCleanup (config, view, STACK_OUT);
+	return psphotReadoutCleanup (config, view, STACK_SRC);
     }
 
     psphotStackObjectsUnifyPosition (objects);
 
-    // measure circular, radial apertures (objects sorted by S/N)
-    psphotRadialAperturesByObject (config, objects, view, STACK_OUT); 
-
     // measure elliptical apertures, petrosians (objects sorted by S/N)
-    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_OUT); // pass 1 (detections->allSources)
+    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
 
     // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
-    psphotExtendedSourceFits (config, view, STACK_OUT); // pass 1 (detections->allSources)
+    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
 
     // calculate source magnitudes
-    psphotMagnitudes(config, view, STACK_OUT);
+    psphotMagnitudes(config, view, STACK_SRC);
+
+    // copy the detections from SRC to OUT (for radial aperture photometry) 
+    if (!psphotCopySources (config, view, STACK_OUT, STACK_SRC)) {
+	psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
+	return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    // XXX need to do something to reassign the source pixels here
+
+    bool smoothAgain = true;
+    for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) {
+	// measure circular, radial apertures (objects sorted by S/N)
+	psphotRadialAperturesByObject (config, objects, view, STACK_OUT, nMatchedPSF); 
+	psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF);
+    }
 
     if (0 && !psphotEfficiency(config, view, STACK_OUT)) {
@@ -179,16 +195,20 @@
 
     // replace background in residual image
-    psphotSkyReplace (config, view, STACK_RAW);
+    psphotSkyReplace (config, view, STACK_DET);
 
     // drop the references to the image pixels held by each source
-    psphotSourceFreePixels (config, view, STACK_OUT);
-
-    // remove chisq image from config->file:PSPHOT.INPUT (why?)
-    psphotStackRemoveChisqFromInputs(config, STACK_RAW);
+    // psphotSourceFreePixels (config, view, STACK_OUT);
+    psphotSourceFreePixels (config, view, STACK_SRC);
+
+    // remove chisq image from config->file:PSPHOT.INPUT
+    psphotStackRemoveChisqFromInputs(config, STACK_DET);
+    if (strcmp(STACK_SRC, STACK_DET)) {
+	psphotStackRemoveChisqFromInputs(config, STACK_SRC);
+    }
 
     psFree (objects);
 
     // create the exported-metadata and free local data
-    return psphotReadoutCleanup (config, view, STACK_OUT);
+    return psphotReadoutCleanup (config, view, STACK_SRC);
 }
 
@@ -258,5 +278,5 @@
  ******
 
- the above is all wrong:  first, we sohould be doing the full
+ the above is all wrong:  first, we should be doing the full
  morphology analysis (ExtendedAnalysis & ExtendedFits) on the CNV or
  RAW image (as desired optionally), etc.