Index: /branches/eam_branches/ipp-20120905/psphot/src/psphot.h
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphot.h	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphot.h	(revision 34570)
@@ -372,5 +372,5 @@
 bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index);
 bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS);
-bool psphotDropBadMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
+bool psphotFilterMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
 
 bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final);
@@ -525,3 +525,6 @@
 bool psphotSourceMemoryReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index);
 
+const char * psphotGetFilerule(const char *baseRule);
+extern bool psphotINpsphotStack;
+
 #endif
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotAddNoise.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotAddNoise.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotAddNoise.c	(revision 34570)
@@ -34,4 +34,5 @@
 static int Nmasked = 0;
 
+// the return state indicates if any sources were actually replaced
 bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add) {
 
@@ -51,7 +52,6 @@
 
     psArray *sources = detections->allSources;
-    //psAssert (sources, "missing sources?");
-    // if no work, should just return true
-    if (!sources) return true;
+    // if no work to do, should just return true
+    if (!sources) return false;
 
     psTimerStart ("psphot.noise");
@@ -74,5 +74,5 @@
 
     if (SIZE <= 0) {
-       return true;
+       return false;
     }
 
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotDeblendSatstars.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotDeblendSatstars.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotDeblendSatstars.c	(revision 34570)
@@ -544,8 +544,10 @@
     float logRdel = 0.1;
 
-    float Xc = source->satstar->Xo - source->pixels->col0 - 0.5;
-    float Yc = source->satstar->Yo - source->pixels->row0 - 0.5;
-    psVector *logRmodel = source->satstar->logRmodel;
-    psVector *logFmodel = source->satstar->logFmodel;
+    pmSourceSatstar *satstar = source->satstar ? source->satstar : (source->parent ? source->parent->satstar : NULL);
+    psAssert (satstar, "null satstar");
+    float Xc = satstar->Xo - source->pixels->col0 - 0.5;
+    float Yc = satstar->Yo - source->pixels->row0 - 0.5;
+    psVector *logRmodel = satstar->logRmodel;
+    psVector *logFmodel = satstar->logFmodel;
 
     for (int iy = 0; iy < source->pixels->numRows; iy++) {
@@ -1048,7 +1050,9 @@
 }
 
+// the return state indicates if any sources were actually subtracted
 bool psphotAddOrSubSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex, psMetadata *recipe, bool add) {
 
     bool status;
+    bool modified = false;
 
     psTimerStart ("psphot.deblend.sat");
@@ -1065,9 +1069,10 @@
 
     psArray *sources = detections->allSources;
-    psAssert (sources, "missing sources?");
+    // if no work to do, should just return false
+    if (!sources) return false;
 
     if (!sources->n) {
 	psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend");
-	return true;
+	return false;
     }
 
@@ -1082,9 +1087,12 @@
 	if (!(source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE)) continue;
 
+	// tell the calling function that we modified the image
+	modified = true;
+
 	if (!psphotSatstarProfileOp (source, maskVal, 1.0, 0, add)) continue;
     }
 
     psLogMsg ("psphot", PS_LOG_DETAIL, "satstar op: %f sec\n", psTimerMark ("psphot.deblend.sat"));
-    return true;
+    return modified;
 }
 
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotDefineFiles.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotDefineFiles.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotDefineFiles.c	(revision 34570)
@@ -1,3 +1,5 @@
 # include "psphotInternal.h"
+
+bool psphotINpsphotStack = false;
 
 // List of output files
@@ -181,2 +183,26 @@
     return;
 }
+
+// psphotGetFilerule
+// Since psphotStack processes multipe FPAs at a time it has a different file rule structure than regular psphot. 
+// For the background output files we define a function psphotGetFilerule which given a base psphot file rule
+// returns the corresponding psphotStack rule *if* the program is psphotStack. That is indicated by a global
+// boolean which defaults to false, and psphotStack only sets to true
+
+const char *psphotGetFilerule(const char *psphotRule) {
+    const char *rule = psphotRule;
+    if (psphotINpsphotStack) {
+        if (!strcmp(psphotRule, "PSPHOT.BACKMDL")) {
+            rule =  "PSPHOT.STACK.BACKMDL";
+        } else if (!strcmp(psphotRule, "PSPHOT.BACKMDL.STDEV")) {
+            rule = "PSPHOT.STACK.BACKMDL.STDEV";
+        } else if (!strcmp(psphotRule, "PSPHOT.BACKSUB")) {
+            rule = "PSPHOT.STACK.BACKSUB";
+        } else if (!strcmp(psphotRule, "PSPHOT.BACKGND")) {
+            rule = "PSPHOT.STACK.BACKGND";
+        } else {
+            psAssert(0, "unsupported file rule %s", psphotRule);
+        }
+    }
+    return rule;
+}
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotFindDetections.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotFindDetections.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotFindDetections.c	(revision 34570)
@@ -71,4 +71,10 @@
     }
 
+    // first pass vs other: if this is the first pass, the code will use PEAKS_NSIGMA_LIMIT and
+    // only attempt to detect PEAKS_NMAX entries.  If 'firstPass' is false, the code will
+    // attempt to replace the subtracted sources in order to measure the footprints.  After
+    // replacement, it is necessary to regenerate the significance image.  If no sources are
+    // available, the code will skip the significance image regeneration step.
+
     bool replaceSourcesForFootprints = false;
     if (firstPass) {
@@ -126,13 +132,16 @@
     if (useFootprints) {
         if (replaceSourcesForFootprints) {
+	    bool modified = false;
 	    // subtract the noise for all sources including satstars
-            psphotAddOrSubNoiseReadout(config, view, filerule, index, recipe, false);
-            psphotReplaceAllSourcesReadout (config, view, filerule, index, recipe, false);
+	    modified |= psphotAddOrSubNoiseReadout(config, view, filerule, index, recipe, false);
+            modified |= psphotReplaceAllSourcesReadout (config, view, filerule, index, recipe, false);
 
 	    // add in the satstars
-	    psphotAddOrSubSatstarsReadout (config, view, filerule, index, recipe, true);
+	    modified |= psphotAddOrSubSatstarsReadout (config, view, filerule, index, recipe, true);
 
-            psFree (significance);
-            significance = psphotSignificanceImage (readout, recipe, maskVal);
+	    if (modified) {
+		psFree (significance);
+		significance = psphotSignificanceImage (readout, recipe, maskVal);
+	    }
 
 	    // display the significance image
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotKronIterate.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotKronIterate.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotKronIterate.c	(revision 34570)
@@ -1,10 +1,9 @@
 # include "psphotInternal.h"
 
-bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow);
-bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels);
-
-
-bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass)
-{
+bool psphotKronRadiusMeasure (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal, psImage *smoothedPixels);
+bool psphotKronFluxMeasure (pmSource *source, psImageMaskType maskVal);
+
+
+bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass) {
     bool status = true;
 
@@ -98,13 +97,4 @@
     }
 
-    bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
-    if (!status) {
-        KRON_APPLY_WEIGHT = true;
-    }
-
-    bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
-    if (!status) {
-        KRON_APPLY_WINDOW = false;
-    }
     bool KRON_SMOOTH = psMetadataLookupBool (&status, recipe, "KRON_SMOOTH");
     if (!status) {
@@ -170,17 +160,8 @@
     }
 
-    // generate the window image: multiply the flux by this to downweight neighbors
-    // XXX: we don't need this image if we aren't going to apply the window
-    psImage *kronWindow = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
-    psImageInit (kronWindow, 1.0);
-
-    // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image
-    for (int i = 0; i < sources->n; i++) {
-
-        pmSource *source = sources->data[i];
-
-	// set a window function for each source based on the moments
-	// (this skips really bad sources (no peak, no moments, DEFECT)
-	psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW);
+    // pass 3 we only measure fluxes for matched sources so we don't need smooth pixels or multiple iterations
+    if (pass == 3) {
+        KRON_SMOOTH = false;
+        KRON_ITERATIONS = 1;
     }
 
@@ -235,5 +216,4 @@
 
             psArrayAdd(job->args, 1, readout);
-            psArrayAdd(job->args, 1, kronWindow);
             psArrayAdd(job->args, 1, cells->data[j]); // sources
             psArrayAdd(job->args, 1, smoothedImage);
@@ -243,6 +223,4 @@
             PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS,    PS_TYPE_F32);
             PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS,    PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32);
             PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
             PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32);
@@ -281,5 +259,4 @@
     }
     psFree (cellGroups);
-    psFree (kronWindow);
     if (KRON_SMOOTH) {
         for (int i = 0; i < sources->n; i++) {
@@ -298,20 +275,19 @@
 
     pmReadout *readout              = job->args->data[0];
-    psImage *kronWindow             = job->args->data[1];
-    psArray *sources                = job->args->data[2];
-    psImage *smoothedImage          = job->args->data[3];
-    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
-    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
-    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[6],F32);
-    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[7],F32);
-    int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[8],S32);
-    bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[9],S32);
-    bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[10],S32);
-    float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[11],F32);
-    float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[12],F32);
-    int pass                        = PS_SCALAR_VALUE(job->args->data[13],S32);
-#ifndef REVERT_ON_BAD_MEASUREMENT
-    (void) pass;
-#endif
+    psArray *sources                = job->args->data[1];
+    psImage *smoothedImage          = job->args->data[2];
+    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
+    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[5],F32);
+    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
+    int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[7],S32);
+    float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[8],F32);
+    float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[9],F32);
+    int pass                        = PS_SCALAR_VALUE(job->args->data[10],S32);
+
+    bool measureRadius = true;
+    if (pass == 3) {
+        measureRadius = false;
+    }
 
     for (int j = 0; j < KRON_ITERATIONS; j++) {
@@ -320,27 +296,35 @@
 	    pmSource *source = sources->data[i];
 	    if (!source->peak) continue; // XXX how can we have a peak-less source?
+            if (!source->moments) continue;
+            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+	    // skip saturated stars modeled with a radial profile 
+	    if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
+
+            // in pass 3 we only measure the flux for matched sources
+            if ((pass == 3) && !(source->mode2 & PM_SOURCE_MODE2_MATCHED)) continue;
 
 # if (0)
 # define TEST_X 653
 # define TEST_Y 466
-	if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
-	    fprintf (stderr, "test object\n");
-	}
+            if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+                fprintf (stderr, "test object\n");
+            }
 # undef TEST_X
 # undef TEST_Y
 # endif
 
-	    // check status of this source's moments
-	    if (!source->moments) continue;
-	    if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
-	    if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+            if (measureRadius) {
+                // check status of this source's moments
+                // XXX: I don't think that we have to apply these restrictions since we dropped the window function
+                if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
+                if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+            }
+
             if (!isfinite(source->moments->Mrf)) {
                 // Once we save a bad Mrf measurement we give up on this source
-                // checking here allows us to avoid adding and subtracting the model
+                // XXX: is this the right thing to do?
                 continue;
             }
-
-	    // skip saturated stars modeled with a radial profile 
-	    if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
 
 	    // replace object in image
@@ -358,5 +342,10 @@
             // On first iteration set window radius to sky radius (if valid). We also use this on subsequent
             // iterations if we cannot find a better limit
-	    float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
+	    float maxWindow;
+            if (measureRadius) {
+                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
+            } else {
+                maxWindow = source->moments->Mrf;
+            }
             if (j > 0) {
                 // on subsequent iterations we use a factor times the previous radial moment value
@@ -380,52 +369,40 @@
 	    float windowRadius = PS_MAX(RADIUS, maxWindow);
 
-#ifdef REVERT_ON_BAD_MEASURMENT
-            // save previous measurements. We might revert back to them if this round fails
-            float MrfPrior = source->moments->Mrf;
-            float KronFluxPrior = source->moments->KronFlux;
-            float KronFluxErrPrior = source->moments->KronFluxErr;
-#endif
-
 	    // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
 	    bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
-	    psAssert (source->pixels, "WTF?");
+	    psAssert (source->pixels, "redefine pixels failed?");
             if (extend && smoothedPixels) {
                 psFree(source->tmpPtr);
                 smoothedPixels = psImageSubset(smoothedImage, source->region);
-                psAssert (smoothedPixels, "WTF?");
+                psAssert (smoothedPixels, "redefine smoothed pixels failed?");
                 source->tmpPtr = (psPtr) smoothedPixels ;
             }
 
-
-	    // clear the window function for this source based on the moments
-            // Note: this function also applies cuts on the source and returns false if it
-            // does not meet the requirements for measuring the Kron Radius or Magnitude.
-            // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false
-	    if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) {
-
-                // this function populates moments->Mrf,KronFlux,KronFluxErr
-                psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, smoothedPixels);
-                psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
-
-                // set a window function for each source based on the moments
-                psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW);
-            }
-
-#ifdef REVERT_ON_BAD_MEASUREMENT
-            // on pass 2 if we get an invalid measurement on a pass 1 source where we had a good one previously
-            // in pass 1 keep that measurement
-            bool reverted = false;
-            if (pass > 1 && !isfinite(source->moments->Mrf) && (source->mode2 & PM_SOURCE_MODE2_PASS1_SRC)) {
-                source->moments->Mrf = MrfPrior;    // This is finite otherwise we wouldn't have gotten here
-                source->moments->KronFlux = KronFluxPrior;
-                source->moments->KronFluxErr = KronFluxErrPrior;
-                reverted = true;
-            }
-#endif
+            bool measureFlux = true;
+            if (measureRadius) {
+                measureFlux = psphotKronRadiusMeasure (source, windowRadius, MIN_KRON_RADIUS, maskVal, smoothedPixels);
+            }
+
+            if (measureFlux) {
+                // Make sure the sources images are large enough for the measured Kron Radius
+                bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, 
+                    2.5 * source->moments->Mrf + 2);
+
+                psAssert (source->pixels, "redefine pixels failed?");
+                if (extend && smoothedPixels) {
+                    psFree(source->tmpPtr);
+                    smoothedPixels = psImageSubset(smoothedImage, source->region);
+                    psAssert (smoothedPixels, "redefine smoothed pixels failed?");
+                    source->tmpPtr = (psPtr) smoothedPixels ;
+                }
+
+                // this function populates moments->Mrf,KronFlux,KronFluxErr, KronFinner, and KronFouter
+                psphotKronFluxMeasure (source, maskVal);
+            }
+
+            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+
 #ifdef DUMP_KRS
-#ifndef REVERT_ON_BAD_MEASUREMENT
-            bool reverted = false;
-#endif
-            fprintf(dumpFile, "%7d %1d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %2d\n", source->id, reverted, source->moments->Mrf, MrfPrior, maxWindow, windowRadius, source->peak->xf, source->peak->yf, source->imageID);
+            fprintf(dumpFile, "%7d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %2d\n", source->id, source->moments->Mrf, MrfPrior, maxWindow, windowRadius, source->peak->xf, source->peak->yf, source->imageID);
 #endif
 
@@ -443,6 +420,6 @@
 }
 
-bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal,
-    bool applyWeight, psImage *smoothedPixels) {
+bool psphotKronRadiusMeasure (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal,
+    psImage *smoothedPixels) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -452,5 +429,4 @@
 
     psF32 R2 = PS_SQR(radius);
-    float rsigma2 =  applyWeight ? 0.5 / R2 : 0;
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -465,5 +441,4 @@
     // Xn  = SUM (x - xc)^n * (z - sky)
 
-
     psF32 RF = 0.0;
     psF32 RS = 0.0;
@@ -481,7 +456,4 @@
     psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
 
-    int Xwo = source->pixels->col0;
-    int Ywo = source->pixels->row0;
-
     psF32 **vPix;
     if (smoothedPixels) {
@@ -490,18 +462,142 @@
         vPix = source->pixels->data.F32;
     }
-    psF32 **vWin = kronWindow->data.F32;
+    
+    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+
+    {
+        for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+            psF32 yDiff = row - yCM;
+            if (fabs(yDiff) > radius) continue;
+
+            // coordinate of mirror pixel
+            int yFlip = yCM - yDiff;
+            if (yFlip < 0) continue;
+            if (yFlip >= source->pixels->numRows) continue;
+
+            for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+                // check mask and value for this pixel
+                if (vMsk && (vMsk[row][col] & maskVal)) continue;
+                if (isnan(vPix[row][col])) continue;
+
+                psF32 xDiff = col - xCM;
+                if (fabs(xDiff) > radius) continue;
+
+                // coordinate of mirror pixel
+                int xFlip = xCM - xDiff;
+                if (xFlip < 0) continue;
+                if (xFlip >= source->pixels->numCols) continue;
+
+                // check mask and value for mirror pixel
+                if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
+                if (isnan(vPix[yFlip][xFlip])) continue;
+
+                // radius is just a function of (xDiff, yDiff)
+                psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+                if (r2 > R2) continue;
+
+                float fDiff1 = vPix[row][col];
+                float fDiff2 = vPix[yFlip][xFlip];
+
+                float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
+
+                // Kron Flux uses the 1st radial moment
+                psF32 rf = pDiff * sqrt(r2);
+                psF32 rs = 0.5 * (fDiff1 + fDiff2);
+
+                RF  += rf;
+                RS  += rs;
+            }
+        }
+
+        float MrfTry = RF/RS;
+        if (RF <= 0. || RS <= 0 || !isfinite(MrfTry)) {
+            // We did not get a good measurement
+            source->moments->Mrf = NAN;
+            source->moments->KronFlux  = NAN;
+            source->moments->KronFluxErr  = NAN;
+            return false;
+        }
+
+        float Mrf = MAX(minKronRadius, MrfTry);
+        // Saturate the 1st radial moment
+        if (sqrt(source->peak->detValue) < 10.0) {
+            Mrf = MIN (radius, Mrf);
+        }
+        source->moments->Mrf = Mrf;
+    }
+    return true;
+}
+
+bool psphotKronFluxMeasure(pmSource *source, psImageMaskType maskVal) {
+
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+
+    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
+    // image coordinates.  the source->pixels image has an offset relative to its parent of
+    // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in
+    // this subimage.  we subtract off the peak coordinates, adjusted to this subimage, to have
+    // minimal round-off error in the sums.  since these values are subtracted just to minimize
+    // the dynamic range and are added back below, the exact value does not matter. these are
+    // (int) so they can be used in the image index below.
+
+    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
+    // Xn  = SUM (x - xc)^n * (z - sky)
+
+    // the peak position is less accurate but less subject to extreme deviations
+    float dX = source->moments->Mx - source->peak->xf;
+    float dY = source->moments->My - source->peak->yf;
+    float dR = hypot(dX, dY);
+    float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
+    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
+
+    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
+    // xCM, yCM from pixel coords to pixel index here.
+    psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
+    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
+
+    psF32 **vPix = source->pixels->data.F32;
     psF32 **vWgt = source->variance->data.F32;
     
     psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
 
+    source->moments->KronFlux  = NAN;
+    source->moments->KronFluxErr  = NAN;
+    source->moments->KronFinner = NAN;
+    source->moments->KronFouter = NAN;
+
+    // Calculate the Kron fluxes
+    float radKinner = 1.0 * source->moments->Mrf;
+    float radKron   = 2.5 * source->moments->Mrf;
+    float radKouter = 4.0 * source->moments->Mrf;
+
+    float limitRadius = MIN (radKouter, source->windowRadius);
+    if (radKouter > source->windowRadius) {
+        // This happens but the measurement isn't important enough to allocate the extra pixels
+        // psWarning ("outer kron radius: %f is larger than windowRadius: %f for %d\n", 
+        //     radKouter, source->windowRadius, source->id);
+        limitRadius = MIN (radKron, source->windowRadius);
+    }
+    if (radKron > source->windowRadius) {
+        // caller should have prevented this from happening
+        psWarning ("kron radius: %f is larger than windowRadius: %f for %d\n", 
+            radKron, source->windowRadius, source->id);
+        return false;
+    }
+
+    float Sum = 0.0;
+    float Var = 0.0;
+    float SumInner = 0.0;
+    float SumOuter = 0.0;
+
+    // set vPix to the source pixels (it may have been set to the smoothed image above)
+    vPix = source->pixels->data.F32;
+
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
 	psF32 yDiff = row - yCM;
-	if (fabs(yDiff) > radius) continue;
-
-	// coordinate of mirror pixel
-	int yFlip = yCM - yDiff;
-	if (yFlip < 0) continue;
-	if (yFlip >= source->pixels->numRows) continue;
+	if (fabs(yDiff) > limitRadius) continue;
 
 	for (psS32 col = 0; col < source->pixels->numCols ; col++) {
@@ -511,186 +607,34 @@
 
 	    psF32 xDiff = col - xCM;
-	    if (fabs(xDiff) > radius) continue;
-
-	    // coordinate of mirror pixel
-	    int xFlip = xCM - xDiff;
-	    if (xFlip < 0) continue;
-	    if (xFlip >= source->pixels->numCols) continue;
-
-	    // check mask and value for mirror pixel
-	    if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
-	    if (isnan(vPix[yFlip][xFlip])) continue;
-
-	    // radius is just a function of (xDiff, yDiff)
+	    if (fabs(xDiff) > limitRadius) continue;
+
 	    psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-	    if (r2 > R2) continue;
-
-	    // flux * window
-	    float z = r2 * rsigma2;
-	    assert (z >= 0.0);
-
-	    // weight by window image and wide Gaussian
-	    float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
-	    float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
-
-	    float fDiff1 = vPix[row][col]*weight1;
-	    float fDiff2 = vPix[yFlip][xFlip]*weight2;
-
-	    float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
-
-	    // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
-	    psF32 rf = pDiff * sqrt(r2);
-	    psF32 rs = 0.5 * (fDiff1 + fDiff2);
-
-	    RF  += rf;
-	    RS  += rs;
+            psF32 r = sqrt(r2);
+
+	    float pDiff = vPix[row][col];
+	    psF32 wDiff = vWgt[row][col];
+
+            if (r > radKinner && r < radKron) {
+                SumInner += pDiff;
+            }
+            if (r < radKron) {
+                Sum += pDiff;
+                Var += wDiff;
+            }
+            if (r > radKron && r < radKouter) {
+                SumOuter += pDiff;
+            }
 	}
     }
 
-    float MrfTry = RF/RS;
-    if (RF <= 0. || RS <= 0 || !isfinite(MrfTry)) {
-        // We did not get a good measurement
-	source->moments->Mrf = NAN;
-        source->moments->KronFlux  = NAN;
-        source->moments->KronFluxErr  = NAN;
-        return false;
-    }
-
-    float Mrf = MAX(minKronRadius, MrfTry);
-    // Saturate the 1st radial moment
-    if (sqrt(source->peak->detValue) < 10.0) {
-	Mrf = MIN (radius, Mrf);
-    }
-
-    // Calculate the Kron magnitude (make this block optional?)
-    float radKron  = 2.5*Mrf;
-    float radKron2 = radKron*radKron;
-
-    int nKronPix = 0;
-    float Sum = 0.0;
-    float Var = 0.0;
-    float Win = 0.0;
-
-    // set vPix to the source pixels (it may have been set to the
-    // smoothed image above)
-    vPix = source->pixels->data.F32;
-
-
-    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
-
-	psF32 yDiff = row - yCM;
-	if (fabs(yDiff) > radKron) continue;
-
-	for (psS32 col = 0; col < source->pixels->numCols ; col++) {
-	    // check mask and value for this pixel
-	    if (vMsk && (vMsk[row][col] & maskVal)) continue;
-	    if (isnan(vPix[row][col])) continue;
-
-	    psF32 xDiff = col - xCM;
-	    if (fabs(xDiff) > radKron) continue;
-
-	    // radKron is just a function of (xDiff, yDiff)
-	    psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-	    if (r2 > radKron2) continue;
-
-	    float weight1  = vWin[row+Ywo][col+Xwo];
-	    float fDiff1 = vPix[row][col]*weight1;
-
-	    float pDiff = fDiff1;
-	    psF32 wDiff = vWgt[row][col] * weight1;
-
-	    Sum += pDiff;
-	    Var += wDiff;
-	    Win += weight1;
-	    nKronPix ++;
-	}
-    }
-
-    source->moments->Mrf         = Mrf;
     source->moments->KronFlux    = Sum;
     source->moments->KronFluxErr = sqrt(Var);
+    source->moments->KronFinner  = SumInner;
+
+    // only save radKouter if the radius is inside the integration radius limit
+    if (radKouter <= limitRadius) {
+        source->moments->KronFouter = SumOuter;
+    }
 
     return true;
 }
-
-bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) {
-
-    if (!source) return false;
-    if (!source->peak) return false; // XXX how can we have a peak-less source?
-    if (!source->moments) return false;
-    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
-    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
-    if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) return false;
-    if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false;
-    psAssert(kronWindow, "need a window");
-
-    // XXX: If we are not applying the window then we don't need to check for valid Mrf here.
-    // We should give the this module a chance to measure a good value. 
-    // However experiments show that it hardly ever succeeds in getting a better value
-    if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
-
-    if (!applyWindow) {
-        return true;
-    }
-
-    // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has 
-    // these values (and peak of 1.0)
-
-    int Nx = kronWindow->numCols;
-    int Ny = kronWindow->numRows;
-
-    float Xo = source->moments->Mx;
-    float Yo = source->moments->My;
-
-    psEllipseMoments moments;
-    moments.x2 = source->moments->Mxx;
-    moments.y2 = source->moments->Myy;
-    moments.xy = source->moments->Mxy;
-
-    psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
-    if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
-        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
-        // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
-        return false;
-    }
-
-    // Why this factor of 0.5 ?
-    float scale = 0.5 * source->moments->Mrf / axes.major;
-    axes.major *= scale;
-    axes.minor *= scale;
-
-    psEllipseShape shape = psEllipseAxesToShape(axes);
-    if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
-        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
-        // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
-        return false;
-    }
-
-    float Smajor = axes.major;
-
-    int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1);
-    int maxX = PS_MIN(PS_MAX(Xo + 5*Smajor, 0), Nx - 1);
-    int minY = PS_MIN(PS_MAX(Yo - 5*Smajor, 0), Ny - 1);
-    int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1);
-
-    float rMxx = 0.5 / PS_SQR(shape.sx);
-    float rMyy = 0.5 / PS_SQR(shape.sy);
-    float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
-                                    // implementation. XXX: Is this correct?
-
-    for (int iy = minY; iy < maxY; iy++) {
-	for (int ix = minX; ix < maxX; ix++) {
-
-	    float dX = (ix + 0.5 - Xo);
-	    float dY = (iy + 0.5 - Yo);
-
-	    float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
-
-	    float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z));
-
-            kronWindow->data.F32[iy][ix] *= f;
-	}
-    }
-
-    return true;
-}
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotMaskBackground.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotMaskBackground.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotMaskBackground.c	(revision 34570)
@@ -17,8 +17,8 @@
 
     // find the currently selected readout (XXX note that the model is saved on PSPHOT.BACKMDL regardless of 'filename'
-    pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest
+    pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, psphotGetFilerule("PSPHOT.BACKMDL"), index); // File of interest
     assert (modelFile);
 
-    pmFPAfile *stdevFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL.STDEV", index);
+    pmFPAfile *stdevFile = pmFPAfileSelectSingle(config->files, psphotGetFilerule("PSPHOT.BACKMDL.STDEV"), index);
     assert (stdevFile);
 
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotModelBackground.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotModelBackground.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotModelBackground.c	(revision 34570)
@@ -405,6 +405,6 @@
 
     psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters
-    pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning, index);
-    pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index);
+    pmReadout *model = pmFPAGenerateReadout(config, view, psphotGetFilerule("PSPHOT.BACKMDL"), inFPA, binning, index);
+    pmReadout *modelStdev = pmFPAGenerateReadout(config, view, psphotGetFilerule("PSPHOT.BACKMDL.STDEV"), inFPA, binning, index);
 
     if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config, false)) {
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotOutput.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotOutput.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotOutput.c	(revision 34570)
@@ -34,5 +34,5 @@
     
 
-    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL");
+    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, psphotGetFilerule("PSPHOT.BACKMDL"));
     if (!file) return NULL;
 
@@ -45,5 +45,5 @@
     bool status;
 
-    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL.STDEV");
+    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, psphotGetFilerule("PSPHOT.BACKMDL.STDEV"));
     if (!file) return NULL;
 
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotReplaceUnfit.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotReplaceUnfit.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotReplaceUnfit.c	(revision 34570)
@@ -43,4 +43,5 @@
 }
 
+// the return state indicates if any sources were actually replaced
 bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) {
 
@@ -58,5 +59,5 @@
     if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) {
         psLogMsg ("psphot", PS_LOG_DETAIL, "skipping replace all sources for input file %d", index);
-        return true;
+        return false;
     }
 
@@ -65,6 +66,6 @@
 
     psArray *sources = detections->allSources;
-    //psAssert (sources, "missing sources?");
-    if (!sources) return true;
+    // if no work to do, should just return false
+    if (!sources) return false;
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotSetThreads.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotSetThreads.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotSetThreads.c	(revision 34570)
@@ -30,5 +30,5 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 14);
+    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 11);
     task->function = &psphotKronIterate_Threaded;
     psThreadTaskAdd(task);
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotSkyReplace.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotSkyReplace.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotSkyReplace.c	(revision 34570)
@@ -36,5 +36,5 @@
 
     // select background pixels, from output background file, or create
-    pmReadout *background = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKGND");
+    pmReadout *background = pmFPAfileThisReadout (config->files, view, psphotGetFilerule("PSPHOT.BACKGND"));
     if (background == NULL) psAbort("background not defined");
 
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotSourceMatch.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotSourceMatch.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotSourceMatch.c	(revision 34570)
@@ -559,9 +559,9 @@
 }
 
-bool psphotDropBadMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
+bool psphotFilterMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
 
     bool status = false;
 
-    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Drop Bad Matched Sources ---");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Filter Matched Sources ---");
 
     psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
@@ -615,4 +615,36 @@
     psFree(dropped);
 
+    // find the "best" Mrf from the detected sources. 
+    // Currently we use the smallest positive value
+    for (int i=0; i< objects->n; i++) {
+        pmPhotObj *obj = objects->data[i]; 
+
+        float minMrf = 1000.;
+        bool hasMatched = false;
+	for (int j = 0; j < obj->sources->n; j++) {
+	    pmSource *source = obj->sources->data[j]; 
+            if (source->mode2 & PM_SOURCE_MODE2_MATCHED) {
+                hasMatched = true;
+                continue;
+            }
+            float Mrf = source->moments->Mrf;
+            if (isfinite(Mrf) && Mrf < minMrf && Mrf > 0) {
+                minMrf = Mrf;
+            }
+        }
+
+        if (!hasMatched || minMrf > 120.) {
+            continue;
+        }
+
+        // set Mrf for matched sources to the value found above
+	for (int j = 0; j < obj->sources->n; j++) {
+	    pmSource *source = obj->sources->data[j]; 
+            if (source->mode2 & PM_SOURCE_MODE2_MATCHED) {
+                source->moments->Mrf = minMrf;
+            }
+        }
+    }
+
     return true;
 }
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotStack.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotStack.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotStack.c	(revision 34570)
@@ -2,4 +2,8 @@
 
 int main (int argc, char **argv) {
+
+    // Set flag to tell certain library functions that we are in the psphotStack program. 
+    // (This is defined in psphotDefineFiles.c)
+    psphotINpsphotStack = true;
 
     // uncomment to turn on memory dumps (move this to an option)
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotStackImageLoop.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotStackImageLoop.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotStackImageLoop.c	(revision 34570)
@@ -43,4 +43,7 @@
     if (!needConvolved) {
         pmFPAfileActivate (config->files, false, "PSPHOT.STACK.INPUT.CNV");
+        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.MASK.CNV");
+        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.VARIANCE.CNV");
+        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.PSF.CNV");
     }
 
@@ -97,7 +100,7 @@
 	    // drop all versions of the internal files
 	    status = true;
-	    status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
-	    status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");
-	    status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");
+	    status &= pmFPAfileDropInternal (config->files, psphotGetFilerule("PSPHOT.BACKMDL"));
+	    status &= pmFPAfileDropInternal (config->files, psphotGetFilerule("PSPHOT.BACKMDL.STDEV"));
+	    status &= pmFPAfileDropInternal (config->files, psphotGetFilerule("PSPHOT.BACKGND"));
 	    if (!status) {
 		psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files");
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotStackParseCamera.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotStackParseCamera.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotStackParseCamera.c	(revision 34570)
@@ -14,4 +14,17 @@
 	return false;
     }
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    bool savePSF = psMetadataLookupBool(&status, recipe, "SAVE.PSF");
+    bool saveBackgroundModel = psMetadataLookupBool(&status, recipe, "SAVE.BACKMDL");
+    bool saveBackground = psMetadataLookupBool(&status, recipe, "SAVE.BACKGND");
+    bool saveBackSub = psMetadataLookupBool(&status, recipe, "SAVE.BACKSUB");
+    bool saveResid = psMetadataLookupBool(&status, recipe, "SAVE.RESID");
+
+    bool saveCnv = psMetadataLookupBool(&status, recipe, "SAVE.CNV");
+    bool saveChisq = psMetadataLookupBool(&status, recipe, "SAVE.CHISQ");
+    bool useRaw = psMetadataLookupBool(&status, recipe, "PSPHOT.STACK.USE.RAW");
+
 
     int nRaw = 0;
@@ -120,4 +133,14 @@
 	    }
 	}
+        pmFPAfile *inputTemplate;
+        if (useRaw) {
+            inputTemplate = rawInputFile;
+        } else {
+            inputTemplate = cnvInputFile;
+        }
+        if (!inputTemplate) {
+            psError(PS_ERR_UNKNOWN, true, "cannot determinte inputTemplate: USE.RAW: %d\n", useRaw);
+            return false;
+        }
 
 	psString sources = psMetadataLookupStr(&status, input, "SOURCES"); // Name of mask
@@ -135,5 +158,4 @@
         }
 	// generate an pmFPAimage for the output convolved image
-	// XXX output of these files should be optional
 	{
 	    // pmFPAfile *srcInputFile = rawInputFile ? rawInputFile : cnvInputFile;
@@ -143,5 +165,5 @@
 		return false;
 	    }
-	    outputImage->save = true;
+	    outputImage->save = saveCnv;
 	    outputImage->fileID = stack_id;		// this is used to generate output names
 
@@ -155,5 +177,5 @@
 		return NULL;
 	    }
-	    outputMask->save = true;
+	    outputMask->save = saveCnv;
 	    outputMask->fileID = stack_id;		// this is used to generate output names
 
@@ -167,5 +189,5 @@
 		return NULL;
 	    }
-	    outputVariance->save = true;
+	    outputVariance->save = saveCnv;
 	    outputVariance->fileID = stack_id;		// this is used to generate output names
 
@@ -178,4 +200,53 @@
 	    outsources->save = true;
 	    outsources->fileID = stack_id;		// this is used to generate output names
+
+            if (savePSF) {
+                pmFPAfile *output = pmFPAfileDefineOutputFromFile (config, inputTemplate, "PSPHOT.STACK.PSF.SAVE");
+                if (!output) {
+                    psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.PSF.SAVE");
+                    return false;
+                }
+                output->save = true;
+                output->fileID = stack_id;
+            }
+            if (saveBackgroundModel) {
+                int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN");
+                int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN");
+                pmFPAfile *output = pmFPAfileDefineFromFile (config, inputTemplate, DX, DY, "PSPHOT.STACK.BACKMDL");
+                if (!output) {
+                    psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.BACKMDL");
+                    return false;
+                }
+
+                output->save = true;
+                output->fileID = stack_id;
+            }
+            if (saveBackground) {
+                pmFPAfile *output = pmFPAfileDefineFromFile (config, inputTemplate, 1, 1, "PSPHOT.STACK.BACKGND");
+                if (!output) {
+                    psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.STACK.BACKGND");
+                    return false;
+                }
+                output->save = true;
+                output->fileID = stack_id;
+            }
+            if (saveBackSub) {
+                pmFPAfile *output = pmFPAfileDefineFromFile (config, inputTemplate, 1, 1, "PSPHOT.STACK.BACKSUB");
+                if (!output) {
+                    psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.STACK.BACKSUB");
+                    return false;
+                }
+                output->save = true;
+                output->fileID = stack_id;
+            }
+            if (saveResid) {
+                pmFPAfile *output = pmFPAfileDefineOutputFromFile (config, inputTemplate, "PSPHOT.STACK.RESID");
+                if (!output) {
+                    psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.RESID");
+                    return false;
+                }
+                output->save = true;
+                output->fileID = stack_id;
+            }
 	}
     }
@@ -191,25 +262,5 @@
     }
 
-    // select the appropriate recipe information
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
-    bool saveCnv = psMetadataLookupBool(&status, recipe, "SAVE.CNV");
-    bool saveChisq = psMetadataLookupBool(&status, recipe, "SAVE.CHISQ");
-
-    // loop over the available readouts
-    for (int i = 0; i < nInputs; i++) {
-	pmFPAfile *file = NULL;
-
-	file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i);
-	file->save = saveCnv;
-
-	file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.MASK", i);
-	file->save = saveCnv;
-
-	file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.VARIANCE", i);
-	file->save = saveCnv;
-    }
-
     // generate an pmFPAimage for the chisqImage
-    // XXX output of these files should be optional
     {
 	pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE");
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotStackReadout.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotStackReadout.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotStackReadout.c	(revision 34570)
@@ -369,6 +369,10 @@
     psMemDump("psfstats");
 
-    // drop matched sources without any useful measurements
-    psphotDropBadMatchedSources (config, view, STACK_SRC, objects);
+    // drop matched sources without any useful measurements and set kron radii for the ones
+    // we decide to keep
+    psphotFilterMatchedSources (config, view, STACK_SRC, objects);
+
+    // measure kron fluxes for the matched sources only
+    psphotKronIterate(config, view, STACK_SRC, 3);
 
     // measure elliptical apertures, petrosians (objects sorted by S/N)
Index: /branches/eam_branches/ipp-20120905/psphot/src/psphotSubtractBackground.c
===================================================================
--- /branches/eam_branches/ipp-20120905/psphot/src/psphotSubtractBackground.c	(revision 34569)
+++ /branches/eam_branches/ipp-20120905/psphot/src/psphotSubtractBackground.c	(revision 34570)
@@ -21,5 +21,5 @@
 
     // find the currently selected readout (XXX note that the model is saved on PSPHOT.BACKMDL regardless of 'filename'
-    pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest
+    pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, psphotGetFilerule("PSPHOT.BACKMDL"), index); // File of interest
     assert (modelFile);
 
@@ -36,5 +36,5 @@
     // select background pixels, from output background file, or create
     // XXX for now, we will only allow a single background image to be generated
-    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND");
+    file = psMetadataLookupPtr (&status, config->files, psphotGetFilerule("PSPHOT.BACKGND"));
     if (file) {
         // we are using PSPHOT.BACKGND as an I/O file: select readout or create
@@ -50,5 +50,5 @@
         }
     } else {
-        background = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKGND", image->numCols, image->numRows, PS_TYPE_F32);
+        background = pmFPAfileDefineInternal (config->files, psphotGetFilerule("PSPHOT.BACKGND"), image->numCols, image->numRows, PS_TYPE_F32);
     }
     psF32 **backData = background->image->data.F32;