Index: branches/czw_branch/20120906/psphot/src/psphot.h
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphot.h	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphot.h	(revision 34772)
@@ -81,5 +81,5 @@
 
 bool            psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule);
-bool            psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
+bool            psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
 
 bool            psphotBasicDeblend (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -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/czw_branch/20120906/psphot/src/psphotAddNoise.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotAddNoise.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotAddNoise.c	(revision 34772)
@@ -1,3 +1,5 @@
 # include "psphotInternal.h"
+
+bool psphotMaskSource(pmSource *source, bool add, psImageMaskType maskVal);
 
 bool psphotAddNoise (pmConfig *config, const pmFPAview *view, const char *filerule) {
@@ -30,4 +32,7 @@
 }
 
+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) {
 
@@ -47,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");
@@ -56,4 +60,7 @@
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     psAssert (maskVal, "missing mask value?");
+
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
+    psAssert (markVal, "missing mask value?");
 
     // increase variance by factor*(object noise):
@@ -67,5 +74,5 @@
 
     if (SIZE <= 0) {
-       return true;
+       return false;
     }
 
@@ -94,4 +101,6 @@
 
 	pmSourceNoiseOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, FACTOR, SIZE, add, maskVal, 0, 0);
+
+	psphotMaskSource (source, add, markVal);
     }
     if (add) {
@@ -100,4 +109,5 @@
         psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
     }
+    fprintf (stderr, "masked %d objects\n", Nmasked);
 
     psphotVisualShowImage (readout);
@@ -105,2 +115,34 @@
     return true;
 }
+
+bool psphotMaskSource(pmSource *source, bool add, psImageMaskType maskVal) {
+
+    if (!source) return false;
+    if (!source->peak) return false; // XXX how can we have a peak-less source?
+    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
+    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
+
+    float Xc = source->peak->xf - source->pixels->col0 - 0.5;
+    float Yc = source->peak->yf - source->pixels->row0 - 0.5;
+
+    psImageMaskType notMaskVal = ~maskVal;
+
+    for (int iy = 0; iy < source->pixels->numRows; iy++) {
+	for (int ix = 0; ix < source->pixels->numCols; ix++) {
+
+	    float radius = hypot (ix - Xc, iy - Yc) ;
+
+	    if (radius > 4) continue;
+
+	    if (add) {
+	      source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= maskVal;
+	    } else {
+	      source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] &= notMaskVal;
+	    }
+	}
+    }
+    Nmasked ++;
+
+    return true;
+}
+
Index: branches/czw_branch/20120906/psphot/src/psphotBlendFit.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotBlendFit.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotBlendFit.c	(revision 34772)
@@ -1,4 +1,6 @@
 # include "psphotInternal.h"
 bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float sigSigma);
+bool psphotBlendSetSourceSatstar (psImage *image, pmSource *source);
+float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -256,10 +258,10 @@
 
 	int TEST_ON = false;
-# if (1)
+# 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");
-	    //  psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
+	    psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
 	    TEST_ON = true;
 	}
@@ -377,4 +379,15 @@
     if (!source) return false;
     if (!source->peak) return false; // XXX how can we have a peak-less source?
+
+# 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");
+	}
+# undef TEST_X
+# undef TEST_Y
+# endif
+
     if (!source->moments) return false;
     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
@@ -385,5 +398,6 @@
     if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) {
 	// find the radius at which we hit something like 0.1*SATURATION
-	return false;
+	psphotBlendSetSourceSatstar (IDimage, source);
+	return true;
     }
 
@@ -431,6 +445,6 @@
     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?
+    float Sxy = shape.sxy;    // factor of -1 is included to match the previous window function
+    // implementation. XXX: Is this correct?
 
     int ID = source->id;
@@ -443,4 +457,5 @@
 
 	    float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
+	    if (z > 2.311) continue;
 
 	    float f = Io*exp(-z);
@@ -453,2 +468,74 @@
     return true;
 }
+
+bool psphotBlendSetSourceSatstar (psImage *IDimage, pmSource *source) {
+
+    float Xo = source->satstar->Xo;
+    float Yo = source->satstar->Yo;
+    psVector *logRmodel = source->satstar->logRmodel;
+    psVector *logFmodel = source->satstar->logFmodel;
+
+    float lPeak = logFmodel->data.F32[0];
+    float fPeak = pow(10.0, lPeak);
+    float threshold = 0.1*fPeak;
+    float logThresh = log10(threshold);
+
+    float radius = NAN;
+
+    // what is the radius for a specific peak fraction?
+    for (int i = 1; i < logFmodel->n; i++) {
+	float logF = logFmodel->data.F32[i];
+	if (!isfinite(logF)) continue;
+
+	float flux = pow(10.0, logF);
+	if (flux > threshold) continue;
+	
+	float logF0 = logFmodel->data.F32[i - 1];
+	float logR0 = logRmodel->data.F32[i - 1];
+
+	float logF1 = logFmodel->data.F32[i];
+	float logR1 = logRmodel->data.F32[i];
+
+	float logR = InterpolateValues (logF0, logR0, logF1, logR1, logThresh);
+	radius = pow(10.0, logR);
+	break;
+    }
+
+    if (!isfinite(radius)) {
+	for (int i = logFmodel->n - 1; i >= 0; i--) {
+	    if (!isfinite(logFmodel->data.F32[i])) continue;
+	    float logR = logRmodel->data.F32[i];
+	    radius = pow(10.0, logR);
+	    break;
+	}
+    }
+
+    if (!isfinite(radius)) return false;
+
+    int Nx = IDimage->numCols;
+    int Ny = IDimage->numRows;
+
+    // region to mask
+    int minX = PS_MIN(PS_MAX(Xo - radius, 0), Nx - 1);
+    int maxX = PS_MIN(PS_MAX(Xo + radius, 0), Nx - 1);
+    int minY = PS_MIN(PS_MAX(Yo - radius, 0), Ny - 1);
+    int maxY = PS_MIN(PS_MAX(Yo + radius, 0), Ny - 1);
+
+    int ID = source->id;
+
+    fprintf (stderr, "Xo,Yo: %f,%f; radius: %f\n", Xo, Yo, radius);
+
+    for (int iy = minY; iy < maxY; iy++) {
+	for (int ix = minX; ix < maxX; ix++) {
+
+	    float dX = (ix - Xo);
+	    float dY = (iy - Yo);
+	    float R = hypot (dX, dY) ;
+	    if (R > radius) continue;
+
+            IDimage->data.S32[iy][ix] = ID;
+	}
+    }
+    return true;
+}
+
Index: branches/czw_branch/20120906/psphot/src/psphotDeblendSatstars.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotDeblendSatstars.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotDeblendSatstars.c	(revision 34772)
@@ -16,9 +16,21 @@
 bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule)
 {
+    bool status = false;
+
     int num = psphotFileruleCount(config, filerule);
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    // perform full extended source non-linear fits?
+    if (!psMetadataLookupBool (&status, recipe, "SUBTRACT_SATSTAR_PROFILE")) {
+        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source fits\n");
+        return true;
+    }
 
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-	if (!psphotDeblendSatstarsReadout (config, view, filerule, i)) {
+	if (!psphotDeblendSatstarsReadout (config, view, filerule, i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for %s entry %d", filerule, i);
 	    return false;
@@ -40,15 +52,12 @@
     4) subtract the radial profile
     
+    There is also support code for calculation of magnitudes and add/subtract the profile (a la pmSourceOp)
+
     TBD:
 
-    * function to replace the radial profile for a source
     * recenter the profile (based on cross-correlation / convolution with 1D profile)
-
-    * raise a bit somewhere for these super saturated stars (if they were so modeled)
-    * consider the subtraction bit : when to raise it?
-
  **/
 
-bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) {
+bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex, psMetadata *recipe) {
 
     int N;
@@ -56,7 +65,4 @@
     bool status;
 
-    // XXX disable this function for now
-    return true;
-
     psTimerStart ("psphot.deblend.sat");
 
@@ -78,8 +84,4 @@
 	return true;
     }
-
-    // select the appropriate recipe information
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
-    psAssert (recipe, "missing recipe?");
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -124,6 +126,6 @@
     int BIG_SIGMA = BIG_RADIUS / 4.0;
 
-    int display = psphotKapaChannel (1);
-    psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 0);
+    // int display = psphotKapaChannel (1);
+    // psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 0);
 
     // examine sources in decreasing SN order
@@ -185,5 +187,5 @@
     }
     // show the image after object have been subtracted
-    psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 1);
+    // psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 1);
 
     psFree (SN);
@@ -451,4 +453,6 @@
     // XXX do something sensible here if the profile is crap
 
+    // replace an existing profile
+    psFree (source->satstar);
     source->satstar = pmSourceSatstarAlloc();
     source->satstar->Xo = Xo;
@@ -494,5 +498,5 @@
     psStatsInit (fluxStats);
 
-    while (logR->data.F32[bin] < lRmax) {
+    while (bin < logR->n && logR->data.F32[bin] < lRmax) {
       if (isfinite(flux->data.F32[bin])) {
 	  psVectorAppend (fluxVals, flux->data.F32[bin]);
@@ -540,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++) {
@@ -1044,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");
@@ -1061,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;
     }
 
@@ -1078,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/czw_branch/20120906/psphot/src/psphotDefineFiles.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotDefineFiles.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotDefineFiles.c	(revision 34772)
@@ -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/czw_branch/20120906/psphot/src/psphotFindDetections.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotFindDetections.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotFindDetections.c	(revision 34772)
@@ -46,4 +46,10 @@
     psAssert (maskVal, "missing mask value?");
 
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
+    psAssert (markVal, "missing mark value?");
+
+    maskVal |= markVal;
+
     // Use the new pmFootprints approach?
     const bool useFootprints = psMetadataLookupBool(NULL, recipe, "USE_FOOTPRINTS");
@@ -64,4 +70,10 @@
 	psMemIncrRefCounter(detections); // so we can free the detections below
     }
+
+    // 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;
@@ -120,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/czw_branch/20120906/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotFitSourcesLinear.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotFitSourcesLinear.c	(revision 34772)
@@ -439,4 +439,9 @@
         model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
         model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
+
+	if (norm->data.F32[i] < MIN_VALID_FLUX) {
+	  fprintf (stderr, "fit out of bounds for %f,%f : %f\n", source->peak->xf, source->peak->yf, norm->data.F32[i]);
+	  model->params->data.F32[PM_PAR_I0] = MIN_VALID_FLUX;
+	}
 
 	// clear the 'mark' pixels so the subtraction covers the full window
Index: branches/czw_branch/20120906/psphot/src/psphotKronIterate.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotKronIterate.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotKronIterate.c	(revision 34772)
@@ -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,17 +296,35 @@
 	    pmSource *source = sources->data[i];
 	    if (!source->peak) continue; // XXX how can we have a peak-less source?
-
-	    // 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 (!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");
+            }
+# undef TEST_X
+# undef TEST_Y
+# endif
+
+            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
@@ -348,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
@@ -370,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
 
@@ -433,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);
@@ -442,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
@@ -455,5 +441,4 @@
     // Xn  = SUM (x - xc)^n * (z - sky)
 
-
     psF32 RF = 0.0;
     psF32 RS = 0.0;
@@ -471,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) {
@@ -480,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++) {
@@ -501,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/czw_branch/20120906/psphot/src/psphotMaskBackground.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotMaskBackground.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotMaskBackground.c	(revision 34772)
@@ -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/czw_branch/20120906/psphot/src/psphotModelBackground.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotModelBackground.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotModelBackground.c	(revision 34772)
@@ -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/czw_branch/20120906/psphot/src/psphotOutput.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotOutput.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotOutput.c	(revision 34772)
@@ -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/czw_branch/20120906/psphot/src/psphotReadout.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotReadout.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotReadout.c	(revision 34772)
@@ -185,5 +185,5 @@
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    psphotFitSourcesLinear (config, view, filerule, false, false); // pass 1 (detections->allSources)
+    psphotFitSourcesLinear (config, view, filerule, false, true); // pass 1 (detections->allSources)
 
     // measure the radial profiles to the sky
@@ -212,5 +212,5 @@
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
-    psphotFitSourcesLinear (config, view, filerule, true, false); // pass 2 (detections->allSources)
+    psphotFitSourcesLinear (config, view, filerule, true, true); // pass 2 (detections->allSources)
 
     // if we only do one pass, skip to extended source analysis
@@ -260,5 +260,5 @@
 
 	// NOTE: apply to ALL sources
-	psphotFitSourcesLinear (config, view, filerule, true, false); // pass 3 (detections->allSources)
+	psphotFitSourcesLinear (config, view, filerule, true, true); // pass 3 (detections->allSources)
     }
 
@@ -300,5 +300,5 @@
 
 	// NOTE: apply to ALL sources
-	psphotFitSourcesLinear (config, view, filerule, true, false); // pass 3 (detections->allSources)
+	psphotFitSourcesLinear (config, view, filerule, true, true); // pass 3 (detections->allSources)
     }
 
Index: branches/czw_branch/20120906/psphot/src/psphotReplaceUnfit.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotReplaceUnfit.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotReplaceUnfit.c	(revision 34772)
@@ -35,12 +35,11 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-	if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, ignoreState)) {
-	    psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
-	    return false;
-	}
-    }
-    return true;
-}
-
+        // ignore return value. False means "sources not subtracted", not a failure condition
+	(void) psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, ignoreState);
+    }
+    return true;
+}
+
+// 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 +57,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 +64,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/czw_branch/20120906/psphot/src/psphotSetThreads.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotSetThreads.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotSetThreads.c	(revision 34772)
@@ -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/czw_branch/20120906/psphot/src/psphotSkyReplace.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotSkyReplace.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotSkyReplace.c	(revision 34772)
@@ -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/czw_branch/20120906/psphot/src/psphotSourceMatch.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotSourceMatch.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotSourceMatch.c	(revision 34772)
@@ -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/czw_branch/20120906/psphot/src/psphotSourceSize.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotSourceSize.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotSourceSize.c	(revision 34772)
@@ -675,6 +675,11 @@
             return false;
         }
-        psFree(readout->mask);
-        readout->mask = newMask;
+	// Copy the new mask pixel values to the old mask array.  NOTE: we cannot replace
+	// the mask pointer because the source->maskView objects point to the old data
+	// area
+	psImage *oldMask = readout->mask;
+	readout->mask = psImageCopy (readout->mask, newMask, PS_TYPE_IMAGE_MASK);
+	psAssert (oldMask == readout->mask, "should have copied the values in-situ");
+        psFree(newMask);
     }
 
Index: branches/czw_branch/20120906/psphot/src/psphotStack.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotStack.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotStack.c	(revision 34772)
@@ -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/czw_branch/20120906/psphot/src/psphotStackImageLoop.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotStackImageLoop.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotStackImageLoop.c	(revision 34772)
@@ -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/czw_branch/20120906/psphot/src/psphotStackMatchPSFs.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotStackMatchPSFs.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotStackMatchPSFs.c	(revision 34772)
@@ -108,16 +108,13 @@
     bool notMatched = psMetadataLookupBool(NULL, readoutOut->analysis, "NOT.PSF.MATCHED");
     if (notMatched) {
-        if (!readoutOut->image) {
-            // our images got freed in the matching process. Recreate them so that they
-            // are available when sources are copied to this readout.
-            readoutOut->image = psImageCopy(NULL, readoutSrc->image, PS_TYPE_F32);
-            if (readoutSrc->variance) {
-                psFree(readoutOut->variance);
-                readoutOut->variance = psImageCopy(NULL, readoutSrc->variance, PS_TYPE_F32);
-            }
-            if (readoutSrc->mask) {
-                psFree(readoutOut->mask);
-                readoutOut->mask = psImageCopy(NULL, readoutSrc->mask, PS_TYPE_IMAGE_MASK);
-            }
+        psFree(readoutOut->image);
+        readoutOut->image = psImageCopy(NULL, readoutSrc->image, PS_TYPE_F32);
+        psFree(readoutOut->variance);
+        if (readoutSrc->variance) {
+            readoutOut->variance = psImageCopy(NULL, readoutSrc->variance, PS_TYPE_F32);
+        }
+        psFree(readoutOut->mask);
+        if (readoutSrc->mask) {
+            readoutOut->mask = psImageCopy(NULL, readoutSrc->mask, PS_TYPE_IMAGE_MASK);
         }
     }
Index: branches/czw_branch/20120906/psphot/src/psphotStackMatchPSFsUtils.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotStackMatchPSFsUtils.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotStackMatchPSFsUtils.c	(revision 34772)
@@ -303,9 +303,14 @@
 	if (!pmSubtractionMatch(NULL, readoutOut, fake, readoutSrc, footprint, stride, regionSize, spacing, threshold, stampSources, stampsName, type, size, order, widthsCopy, orders, inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
             int errorCode = psErrorCodeLast();
-            if (errorCode == PM_ERR_SMALL_AREA) {
+            if (errorCode == PM_ERR_SMALL_AREA || errorCode == PM_ERR_STAMPS) {
                 // failed to match but, don't fault. Just drop this input from measurements that need the
                 // matched readout.
+                psErrorClear();
                 rejectReadout = true;
-                psLogMsg("psphot", PS_LOG_WARN, "PSF match failed for input %d. Too few pixels.", index);
+                if (errorCode == PM_ERR_SMALL_AREA) {
+                    psLogMsg("psphot", PS_LOG_WARN, "PSF match failed for input %d. Too few pixels.", index);
+                } else if (errorCode == PM_ERR_STAMPS) {
+                    psLogMsg("psphot", PS_LOG_WARN, "PSF match failed for input %d. Failed to find stamps.", index);
+                }
             } else {
                 psError(psErrorCodeLast(), false, "Unable to match images.");
@@ -324,6 +329,8 @@
     
     if (rejectReadout) {
+        psMetadataAddBool(readoutSrc->analysis, PS_LIST_TAIL, "NOT.PSF.MATCHED", PS_META_REPLACE, "", true);
         psMetadataAddBool(readoutOut->analysis, PS_LIST_TAIL, "NOT.PSF.MATCHED", PS_META_REPLACE, "", true);
     } else {
+        psMetadataAddBool(readoutSrc->analysis, PS_LIST_TAIL, "NOT.PSF.MATCHED", PS_META_REPLACE, "", false);
         psMetadataAddBool(readoutOut->analysis, PS_LIST_TAIL, "NOT.PSF.MATCHED", PS_META_REPLACE, "", false);
     }
Index: branches/czw_branch/20120906/psphot/src/psphotStackParseCamera.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotStackParseCamera.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotStackParseCamera.c	(revision 34772)
@@ -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/czw_branch/20120906/psphot/src/psphotStackReadout.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotStackReadout.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotStackReadout.c	(revision 34772)
@@ -169,4 +169,22 @@
     }
 
+    // If DET and SRC are different images, subtract radial profiles for the convolved
+    // image first.  The profiles found on the convolved image will be replaced by those
+    // found for the unconvolved image.  Downstream, in psphotFindDetections, we will
+    // replace the the profiles (and re-subtract them) for the detection image, so we want
+    // to keep those versions of the profiles on the sources.
+    if (strcmp(STACK_SRC, STACK_DET)) {
+      // find and subtract radial profile models for saturated stars (XXX change name eventually)
+      if (!psphotDeblendSatstars (config, view, STACK_SRC)) {
+	psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
+	return psphotReadoutCleanup (config, view, STACK_SRC);
+      }
+    }
+    // find and subtract radial profile models for saturated stars (XXX change name eventually)
+    if (!psphotDeblendSatstars (config, view, STACK_DET)) {
+	psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
+	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
@@ -206,8 +224,8 @@
     // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
     // but iterates to an appropriately larger size
-        logMemStats("before.kron.1");
+    logMemStats("before.kron.1");
     psphotKronIterate(config, view, STACK_SRC, 1);
-        logMemStats("after.kron.1");
-
+    logMemStats("after.kron.1");
+	
     // identify CRs and extended sources
     psphotSourceSize (config, view, STACK_SRC, true);
@@ -332,12 +350,7 @@
     psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
 
-    if (splitLinearFit) {
-        // NOTE: apply to Matched sources. Since the sources that we fit above are subtracted, they will
-        // not be included in this fit.
-        psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
-    } else {
-        // Fit all sources together
-        psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 3 (detections->allSources)
-    }
+    // final linear fit. NOTE: if splitLinearFit is true above, this pass will only fit
+    // the unsubtracted (matched) sources (the sources that we fit above are subtracted)
+    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
 
     // measure the radial profiles to the sky (only measures new objects)
@@ -356,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)
@@ -392,9 +409,10 @@
             }
         }
+        // mark any inputs that we want to skip the matched apertures for
+        psphotStackSetInputsToSkip(config, view, STACK_CNV, true);
+        psphotStackSetInputsToSkip(config, view, STACK_OUT, true);
         psphotRadialApertures (config, view, STACK_CNV, 0); // entry 0 == unmatched
         psMemDump("extmeas");
 
-        // mark any inputs that we want to skip the matched apertures for
-        psphotStackSetInputsToSkip(config, view, STACK_OUT, true);
 
         int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
@@ -422,4 +440,5 @@
         }
     }
+    psphotStackSetInputsToSkip(config, view, STACK_CNV, false);
     psphotStackSetInputsToSkip(config, view, STACK_OUT, false);
 
Index: branches/czw_branch/20120906/psphot/src/psphotSubtractBackground.c
===================================================================
--- branches/czw_branch/20120906/psphot/src/psphotSubtractBackground.c	(revision 34410)
+++ branches/czw_branch/20120906/psphot/src/psphotSubtractBackground.c	(revision 34772)
@@ -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;