Index: branches/eam_branches/ipp-20130904/psphot/src/psphot.h
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphot.h	(revision 36197)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphot.h	(revision 36198)
@@ -479,5 +479,5 @@
 bool psphotSourceChildren (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);
 bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
-psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc);
+psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc, bool sourcesSubtracted);
 
 bool psphotSourceParents (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);
Index: branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceAnalysis.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceAnalysis.c	(revision 36197)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceAnalysis.c	(revision 36198)
@@ -1,3 +1,4 @@
 # include "psphotInternal.h"
+void psphotRadialProfileShowSkips ();
 
 // measure the elliptical radial profile and use this to measure the petrosian parameters for the sources
@@ -33,8 +34,30 @@
 }
 
+/*** for the moment, this test code : it is not thread safe ***/
+int    Nall = 0;
+int  Nskip1 = 0;
+int  Nskip2 = 0;
+int  Nskip3 = 0;
+int  Nskip4 = 0;
+int  Nskip5 = 0;
+int  Nskip6 = 0;
+int  Nskip7 = 0;
+int  Nskip8 = 0;
+int  Nskip9 = 0;
+int Nskip10 = 0;
+int Nskip11 = 0;
+int Nskip12 = 0;
+int Nskip13 = 0;
+int Nskip14 = 0;
+
+# define SKIP(VALUE) { VALUE++; continue; }
+
 // aperture-like measurements for extended sources
 bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
 
     bool status;
+    int NfaintEXT = 0;
+    int NfaintPSF = 0;
+
     int Next = 0;
     int Npetro = 0;
@@ -108,7 +131,9 @@
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NfaintEXT
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NfaintPSF
 
 // set this to 0 to run without threading
-# if (1)	    
+# if (0)	    
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
@@ -129,4 +154,8 @@
 	    scalar = job->args->data[7];
 	    Nannuli += scalar->data.S32;
+	    scalar = job->args->data[8];
+	    NfaintEXT += scalar->data.S32;
+	    scalar = job->args->data[9];
+	    NfaintPSF += scalar->data.S32;
 	    psFree(job);
 # endif
@@ -153,4 +182,8 @@
 		scalar = job->args->data[7];
 		Nannuli += scalar->data.S32;
+		scalar = job->args->data[8];
+		NfaintEXT += scalar->data.S32;
+		scalar = job->args->data[9];
+		NfaintPSF += scalar->data.S32;
             }
             psFree(job);
@@ -163,4 +196,22 @@
     psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
+    psLogMsg ("psphot", PS_LOG_INFO, "  skipped: %d EXT, %d PSF\n", NfaintEXT, NfaintPSF);
+
+    fprintf (stderr, "ext analysis skipped @ 1  : %d\n", Nskip1);
+    fprintf (stderr, "ext analysis skipped @ 2  : %d\n", Nskip2);
+    fprintf (stderr, "ext analysis skipped @ 3  : %d\n", Nskip3);
+    fprintf (stderr, "ext analysis skipped @ 4  : %d\n", Nskip4);
+    fprintf (stderr, "ext analysis skipped @ 5  : %d\n", Nskip5);
+    fprintf (stderr, "ext analysis skipped @ 6  : %d\n", Nskip6);
+    fprintf (stderr, "ext analysis skipped @ 7  : %d\n", Nskip7);
+    fprintf (stderr, "ext analysis skipped @ 8  : %d\n", Nskip8);
+    fprintf (stderr, "ext analysis skipped @ 9  : %d\n", Nskip9);
+    fprintf (stderr, "ext analysis skipped @ 10 : %d\n", Nskip10);
+    fprintf (stderr, "ext analysis skipped @ 11 : %d\n", Nskip11);
+    fprintf (stderr, "ext analysis skipped @ 12 : %d\n", Nskip12);
+    fprintf (stderr, "ext analysis skipped @ 13 : %d\n", Nskip13);
+    fprintf (stderr, "ext analysis skipped @ 14 : %d\n", Nskip14);
+
+    psphotRadialProfileShowSkips ();
 
     psphotVisualShowResidualImage (readout, false);
@@ -177,4 +228,7 @@
 
     bool status;
+
+    int NfaintEXT = 0;
+    int NfaintPSF = 0;
 
     int Next = 0;
@@ -207,22 +261,24 @@
 	pmSource *source = sources->data[i];
 
+	Nall ++;
+
 	// if we have checked the source validity on the basis of the object set, then 
 	// we either skip these tests below or we skip the source completely
-	if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) continue;
+	if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) SKIP (Nskip1);
 	if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource;
 
 	// skip PSF-like and non-astronomical objects
-	if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-	if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-	if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
-	if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+	if (source->type == PM_SOURCE_TYPE_DEFECT) SKIP (Nskip2);
+	if (source->type == PM_SOURCE_TYPE_SATURATED) SKIP (Nskip3);
+	if (source->mode & PM_SOURCE_MODE_DEFECT) SKIP (Nskip4);
+	if (source->mode & PM_SOURCE_MODE_SATSTAR) SKIP (Nskip5);
 
 	// skip saturated stars modeled with a radial profile 
-        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) SKIP (Nskip6);
 
 	// optionally allow non-extended objects to get petrosians as well
 	if (!doPetroStars) {
-	    if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
-	    if (source->type == PM_SOURCE_TYPE_STAR) continue;
+	    if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) SKIP (Nskip7);
+	    // if (source->type == PM_SOURCE_TYPE_STAR) SKIP (Nskip8); -- XXX this might be set for input lists (instead of EXT_LIMIT)
 	}
 
@@ -233,14 +289,16 @@
 	if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
 	    skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
+	    NfaintEXT ++;
 	} else {
 	    skipSource = (sqrt(source->peak->detValue) < SN_LIM);
-	}
-	if (skipSource) continue;
+	    NfaintPSF ++;
+	}
+	if (skipSource) SKIP (Nskip9);
 
 	// limit selection by analysis region (this automatically apply
-	if (source->peak->x < region->x0) continue;
-	if (source->peak->y < region->y0) continue;
-	if (source->peak->x > region->x1) continue;
-	if (source->peak->y > region->y1) continue;
+	if (source->peak->x < region->x0) SKIP (Nskip10);
+	if (source->peak->y < region->y0) SKIP (Nskip11);
+	if (source->peak->x > region->x1) SKIP (Nskip12);
+	if (source->peak->y > region->y1) SKIP (Nskip13);
 
     keepSource:
@@ -266,5 +324,5 @@
 		psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
 		pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-		continue;
+		SKIP (Nskip14);
 	    }
 	    Nannuli ++;
@@ -304,4 +362,10 @@
     scalar->data.S32 = Nannuli;
 
+    scalar = job->args->data[8];
+    scalar->data.S32 = NfaintEXT;
+
+    scalar = job->args->data[9];
+    scalar->data.S32 = NfaintPSF;
+
     return true;
 }
Index: branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c	(revision 36197)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c	(revision 36198)
@@ -475,4 +475,16 @@
         }
 
+	bool testObject = false;
+	// testObject |= ((fabs(source->peak->xf -  179) < 5) && (fabs(source->peak->yf - 1138) < 5));
+	// testObject |= ((fabs(source->peak->xf - 5047) < 5) && (fabs(source->peak->yf -  151) < 5));
+	// testObject |= ((fabs(source->peak->xf - 3929) < 5) && (fabs(source->peak->yf - 4109) < 5));
+	// testObject |= ((fabs(source->peak->xf -  915) < 5) && (fabs(source->peak->yf - 5998) < 5));
+	// testObject |= ((fabs(source->peak->xf - 5406) < 5) && (fabs(source->peak->yf -  326) < 5));
+	if (testObject) {
+	    fprintf (stderr, "test object @ %f, %f\n", source->peak->xf, source->peak->yf);
+	    psTraceSetLevel ("psModules.objects.pmPCM_MinimizeChisq", 5);
+	    psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 5);
+	}
+
         // loop here over the models chosen for each source (exclude by S/N)
         // Reset the iterator
@@ -531,8 +543,4 @@
 		  psFree (source->modelFlux);
 		  source->modelFlux = NULL;
-		  // if ((fabs(source->peak->xf - 2572) < 20) && (fabs(source->peak->yf - 5874) < 20)) {
-		  //   fprintf (stderr, "test object\n");
-		  // }
-
 		  modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal);
 		  if (!modelFit) {
@@ -646,4 +654,9 @@
         psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);
         psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+
+	if (testObject) {
+	    psTraceSetLevel ("psModules.objects.pmPCM_MinimizeChisq", 0);
+	    psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 0);
+	}
     }
 
Index: branches/eam_branches/ipp-20130904/psphot/src/psphotMergeSources.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotMergeSources.c	(revision 36197)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotMergeSources.c	(revision 36198)
@@ -844,5 +844,5 @@
 // array containing the child sources.  XXX currently, this is only used by psphotStackReadout
 // (sources go on allSources so that psphotChoosePSF can be called repeatedly)
-psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc) {
+psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc, bool sourcesSubtracted) {
 
     bool status;
@@ -951,5 +951,10 @@
 	    // child sources have not been subtracted in this image, but this flag may be raised if
 	    // they were subtracted in the parent's image
-	    sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+	    // XXX NOTE : in the pre-20130914 version of psphotStack, we carried a copy of the pixels 
+	    // generated before the subtraction took place (and then we smoothed to match the desired PSF).  
+	    // in the new version, we copy the image after subtraction; we need to distinguish these cases
+	    if (!sourcesSubtracted) {
+		sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+	    }
 
 	    // set the output detections:
Index: branches/eam_branches/ipp-20130904/psphot/src/psphotPetrosianRadialBins.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotPetrosianRadialBins.c	(revision 36197)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotPetrosianRadialBins.c	(revision 36198)
@@ -36,5 +36,5 @@
     pmSourceRadialProfile *profile = source->extpars->petProfile;
 
-    float skyModelErrorSQ = PS_SQR(skynoise);
+    // float skyModelErrorSQ = PS_SQR(skynoise);
 
     int nMax = radiusMax;
@@ -133,16 +133,17 @@
 	if (radius->data.F32[i] > Rmax) {
 	    // calculate the value for the nOut bin
-	    float value, dvalue;
+	    float value; // , dvalue;
 	    if (values->n > 0) {
 		psVectorStats (stats, values, NULL, NULL, 0);
 		value = stats->robustMedian;
-		dvalue = stats->robustStdev;
+		// dvalue = stats->robustStdev;
 	    } else {
 		value = NAN;
-		dvalue = NAN;
+		// dvalue = NAN;
 	    }
 
 	    binSB->data.F32[nOut] = value;
-	    binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ);
+	    // binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ);
+	    binSBstdev->data.F32[nOut] = skynoise / sqrt(values->n);
 	    binFill->data.F32[nOut] = values->n / binArea->data.F32[nOut];
 
Index: branches/eam_branches/ipp-20130904/psphot/src/psphotPetrosianStats.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotPetrosianStats.c	(revision 36197)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotPetrosianStats.c	(revision 36198)
@@ -6,4 +6,5 @@
 // generate the Petrosian radius and flux from the mean surface brightness (r_i)
 
+float InterpolateValuesQuadratic (float *Xin, float *Yin, float X);
 float InterpolateValues     (float X0, float Y0, float X1, float Y1, float X);
 float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1);
@@ -53,5 +54,5 @@
     float dFsum2 = 0.0;
 
-    float nSigma = 3.0;
+    float nSigma = 2.0;
     int lowestSignificantRadius = 0;
     float lowestSignificantRatio = 1.0;
@@ -118,6 +119,15 @@
 		petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[nOut]);
 	    } else {
-		petRadius    = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
-		petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]);
+	      // petRadius    = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
+	      if (nOut > 1) {
+		petRadius    = InterpolateValuesQuadratic (&petRatio->data.F32[nOut-2],   &refRadius->data.F32[nOut-2],   PETROSIAN_RATIO);
+	      } else {
+		petRadius    = InterpolateValuesQuadratic (&petRatio->data.F32[nOut-3], &refRadius->data.F32[nOut-3], PETROSIAN_RATIO);
+	      }
+	      float petRadiusLinear = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
+	      if (fabs(petRadius - petRadiusLinear) > fabs(refRadius->data.F32[nOut] - refRadius->data.F32[nOut-1])) { 
+		fprintf (stderr, "big difference : %f vs %f\n", petRadius, petRadiusLinear); 
+	      }
+	      petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]);
 	    }
 	    above = false;
@@ -235,4 +245,23 @@
 }
 
+// Lagrange's form of the interpolating polynomial...
+float InterpolateValuesQuadratic (float *Xin, float *Yin, float X) {
+
+  float dx01 = Xin[0] - Xin[1];
+  float dx02 = Xin[0] - Xin[2];
+  float dx12 = Xin[1] - Xin[2];
+
+  float dx0  = X - Xin[0];
+  float dx1  = X - Xin[1];
+  float dx2  = X - Xin[2];
+
+  float y0 = Yin[0]*dx1*dx2/(dx01*dx02);
+  float y1 = Yin[1]*dx0*dx2/(dx01*dx12); // need - sign
+  float y2 = Yin[2]*dx0*dx1/(dx02*dx12); 
+
+  float Y = y0 - y1 + y2;
+  return Y;
+}
+
 float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) {
     float dydx = (Y1 - Y0) / (X1 - X0);
Index: branches/eam_branches/ipp-20130904/psphot/src/psphotRadialProfile.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotRadialProfile.c	(revision 36197)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotRadialProfile.c	(revision 36198)
@@ -1,3 +1,11 @@
 # include "psphotInternal.h"
+
+static int Nskip1 = 0;
+static int Nskip2 = 0;
+static int Nskip3 = 0;
+static int Nskip4 = 0;
+static int Nskip5 = 0;
+
+# define SKIP(VALUE) { VALUE++; return false; }
 
 bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) {
@@ -23,5 +31,5 @@
     if (!psphotRadialProfilesByAngles (source, Nsec, Rmax)) {
 	psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian");
-	return false;
+	SKIP (Nskip1);
     }
     // allocate: extpars->radFlux->radii,fluxes,theta
@@ -32,5 +40,5 @@
     if (!psphotRadiiFromProfiles (source, fluxMin, fluxMax)) {
 	psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles");
-	return false;
+	SKIP (Nskip2);
     }
     // allocate : extpars->radFlux->isophotalRadii (use profile->radii,fluxes)
@@ -40,5 +48,5 @@
     if (!psphotEllipticalContour (source)) {
 	// psLogMsg ("psphot", 3, "failed to measure elliptical contour");
-	return false;
+	SKIP (Nskip3);
     }
     // use extpars->radFlux->isophotalRadii,theta (result in extpars->axes)
@@ -48,5 +56,5 @@
     if (!psphotEllipticalProfile (source, RAW_RADIUS)) {
 	psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
-	return false;
+	SKIP (Nskip4);
     }
     // allocate extpars->ellipticalFlux->radiusElliptical,fluxElliptical (use axes to scale raw pixels)
@@ -55,5 +63,5 @@
     if (!psphotRadialBins (recipe, source, Rmax, skynoise)) {
 	psError (PS_ERR_UNKNOWN, false, "failed to generate radial bins");
-	return false;
+	SKIP (Nskip5);
     }
     // allocate extpars->radProfile->binSB, binSBstdv, binSum, binFill, radialBins, area (small lengths)
@@ -62,2 +70,11 @@
     return true;
 }
+
+void psphotRadialProfileShowSkips () {
+  
+  fprintf (stderr, "radial profile skipped @ 1  : %d\n", Nskip1);
+  fprintf (stderr, "radial profile skipped @ 2  : %d\n", Nskip2);
+  fprintf (stderr, "radial profile skipped @ 3  : %d\n", Nskip3);
+  fprintf (stderr, "radial profile skipped @ 4  : %d\n", Nskip4);
+  fprintf (stderr, "radial profile skipped @ 5  : %d\n", Nskip5);
+}
Index: branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c	(revision 36197)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c	(revision 36198)
@@ -589,11 +589,17 @@
     if (TIMING) { psTimerStart ("psphotFitPCM"); }
 
+    // if we are ever (in a given psphot implementation) going to fit a parameter, we must set the options here to include
+    // that parameter (otherwise pmPCMupdate will fail to allocate the dmodelFlux image)
+    // thus, if the sersic analysis below uses an index fit, need to use this EXT_AND_SKY mode for init
+
+    options.mode = PM_SOURCE_FIT_EXT;
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-	options.mode = PM_SOURCE_FIT_NO_INDEX; // XXX note that there may be a conflict with psphotExtendedSourceFits.c:133
 	options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
-	// if we are ever (in a given psphot implementation) going to fit a parameter, we must set the options here to include
-	// that parameter (otherwise pmPCMupdate will fail to allocate the dmodelFlux image)
-	// thus, if the sersic analysis below uses an index fit, need to use this EXT_AND_SKY mode for init
-    } else {
+    }
+    if (modelType == pmModelClassGetType("PS_MODEL_DEV")) {
+	options.mode = PM_SOURCE_FIT_SHAPE;
+	options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
+    }
+    if (modelType == pmModelClassGetType("PS_MODEL_EXP")) {
 	options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
     }
Index: branches/eam_branches/ipp-20130904/psphot/src/psphotStackImageLoop.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotStackImageLoop.c	(revision 36197)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotStackImageLoop.c	(revision 36198)
@@ -22,7 +22,5 @@
 
     pmFPAview *view = pmFPAviewAlloc (0);
-    pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
-    pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV");
-    pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+    pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
 
     if (!input) {
@@ -32,23 +30,7 @@
 
     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-    if (useRaw && inputRaw == NULL) {
-        psLogMsg ("psphot", 1, "PSPHOT.STACK.USE.RAW set but no raw input.");
-        useRaw = false;
-    }
-
-    bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES");
-    bool match_psfs = radial_apertures;
-    bool needConvolved = radial_apertures || !useRaw;
-    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");
-    }
 
     // just load the full set of images up front except for EXPNUM which we defer
     pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.RAW");
-    pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.CNV");
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
 
@@ -71,17 +53,8 @@
 		psMemDump("load");
 
-                if (match_psfs) {
-                    // Generate the 1st PSF-matched image set (larger target PSFs are generated by smoothing this image)
-                    if (!psphotStackMatchPSFs (config, view)) {
-                        psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
-                        psFree (view);
-                        return false;
-                    }
-                } else {
-                    if (!psphotStackAllocateOutput (config, view, recipe)) {
-                        psError(psErrorCodeLast(), false, "failure in psphotStackAllocateOutput for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
-                        psFree (view);
-                        return false;
-                    }
+		if (!psphotStackAllocateOutput (config, view, recipe)) {
+		    psError(psErrorCodeLast(), false, "failure in psphotStackAllocateOutput for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
+		    psFree (view);
+		    return false;
                 }
 		psMemDump("stackmatch");
@@ -126,5 +99,5 @@
 
     // Load the appropriate EXPNUM image
-    pmFPAfileActivate (config->files, true, useRaw ? "PSPHOT.STACK.EXPNUM.RAW" : "PSPHOT.STACK.EXPNUM.CNV");
+    pmFPAfileActivate (config->files, true, "PSPHOT.STACK.EXPNUM.RAW");
 
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa EXPNUM in psphot.");
@@ -146,5 +119,5 @@
                 if (! readout->data_exists) { continue; }
 
-                if (!psphotSetNFrames (config, view, useRaw ? inputRaw->name : inputCnv->name)) ESCAPE ("failed to setNFrames.");
+                if (!psphotSetNFrames (config, view, input->name)) ESCAPE ("failed to setNFrames.");
             }
         }
@@ -167,7 +140,4 @@
 bool UpdateHeadersForReadout (pmConfig *config, pmFPAview *view) {
 
-    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-
     int num = psphotFileruleCount(config, "PSPHOT.INPUT");
 
@@ -179,7 +149,5 @@
 	psAssert (output, "missing file?");
 
-	pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-	pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-	pmFPAfile *input = useRaw ? inputRaw : inputCnv;
+	pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
 	psAssert (input, "missing input file");
 
@@ -205,7 +173,5 @@
 	psAssert (output, "missing file?");
 
-	pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-	pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-	pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+	pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
 	psAssert (input, "missing input file");
 
@@ -242,7 +208,5 @@
 	psAssert (output, "missing file?");
 
-	pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-	pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-	pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+	pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
 	psAssert (input, "missing input file");
 
Index: branches/eam_branches/ipp-20130904/psphot/src/psphotStackMatchPSFsNext.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotStackMatchPSFsNext.c	(revision 36197)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotStackMatchPSFsNext.c	(revision 36198)
@@ -43,5 +43,5 @@
     for (int i = 0; i < num; i++) {
 	if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, lastSize)) {
-            psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF", filerule, i);
+	    psLogMsg ("psphot", PS_LOG_INFO, "failed to smooth image %s (%d) to target PSF", filerule, i);
 	    psImageConvolveSetThreads(oldThreads);
 	    return false;
@@ -102,5 +102,7 @@
 
     if (targetFWHM <= currentFWHM) {
-	psError (PSPHOT_ERR_CONFIG, true, "target FWHM cannot be smaller than current FWHM");
+	// psError (PSPHOT_ERR_CONFIG, true, "target FWHM cannot be smaller than current FWHM");
+	psLogMsg ("psphot", PS_LOG_INFO, "target FWHM (%f) is smaller than current FWHM (%f), not smoothing\n", targetFWHM, currentFWHM);
+	fwhmValues->data.F32[lastSize + 1] = currentFWHM;
 	return false;
     }
Index: branches/eam_branches/ipp-20130904/psphot/src/psphotStackReadout.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotStackReadout.c	(revision 36197)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotStackReadout.c	(revision 36198)
@@ -1,16 +1,43 @@
 # include "psphotInternal.h"
 
+static bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *fPSF);
+static bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *filerule, const char *fPSF, int index);
 static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule);
 static void logMemStats(const char *heading);
 
-// we have 3 possible real filesets:
+// relevant filesets:
 # define STACK_RAW "PSPHOT.STACK.INPUT.RAW"
-# define STACK_CNV "PSPHOT.STACK.INPUT.CNV"
-# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"  /* the psf-matched image */
-
-// we have 3 files on which we operate:
-// DET (detection image)       : nominally RAW (optionally CNV?)
-// SRC (source analysis image) : nominally CNV (optionally RAW)
-// OUT (psf-matched images)    : always OUT
+# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"
+
+// XXX STACK_OUT currently is a copy of STACK_RAW, but should be a pointer to is as in psphot (single)
+
+// TEST CODE, can be removed
+bool psphotDumpImages (pmConfig *config, const pmFPAview *view, const char *filerule, char *base) {
+
+    // XXX do nothing
+    return true;
+
+    int num = psphotFileruleCount(config, "PSPHOT.INPUT");
+
+    for (int i = 0; i < num; i++) {
+	// find the currently selected readout
+	pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
+	psAssert (file, "missing file?");
+
+	pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+	psAssert (readout, "missing readout?");
+
+	char line[256];
+	snprintf (line, 256, "%s.%d.im.fits", base, i);
+	psphotSaveImage (NULL, readout->image, line);
+
+	snprintf (line, 256, "%s.%d.wt.fits", base, i);
+	psphotSaveImage (NULL, readout->variance, line);
+
+	snprintf (line, 256, "%s.%d.mk.fits", base, i);
+	psphotSaveImage (NULL, readout->mask, line);
+    }
+    return true;
+}
 
 bool psphotStackVisualFilerule(pmConfig *config, const pmFPAview *view, const char *filerule) {
@@ -67,18 +94,12 @@
     psAssert (breakPt, "configuration error: set BREAK_POINT");
 
-    // we have 3 relevant files: RAW (unconvolved), CNV (convolved stack), OUT (psf-matched stack)
-    // select which image (RAW or CNV) is used for analysis (RAW always used for detection)
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
-    char *STACK_DET = STACK_RAW;
-
     // load WCS 
-    if (!psphotStackLoadWCS(config, view, STACK_SRC)) {
-        psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_SRC);
+    if (!psphotStackLoadWCS(config, view, STACK_RAW)) {
+        psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_RAW);
         return false;
     }
 
     // set the photcode for each image
-    if (!psphotAddPhotcode (config, view, STACK_SRC)) {
+    if (!psphotAddPhotcode (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
         return false;
@@ -86,77 +107,64 @@
 
     // Generate the mask and weight images (if not supplied) and set mask bits. 
-    // This also insures that all invalid pixels are masked (this is done for STACK_CNV in psphotStackMatchPSFs)
-    if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (!psphotSetMaskAndVariance (config, view, STACK_OUT)) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+    if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "NOTHING")) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // generate a background model (median, smoothed image)
-    if (!psphotModelBackground (config, view, STACK_DET)) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (!psphotSubtractBackground (config, view, STACK_DET)) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (strcmp(STACK_SRC, STACK_DET)) {
-#define MODEL_BACKGROUND_SRC 1
-#ifdef MODEL_BACKGROUND_SRC
-        // work around the fact that the background levels on the convolved
-        // and unconvolved stacks can be different
-        if (!psphotModelBackground (config, view, STACK_SRC)) {
-            return psphotReadoutCleanup (config, view, STACK_SRC);
-        }
-#endif
-        if (!psphotSubtractBackground (config, view, STACK_SRC)) {
-            return psphotReadoutCleanup (config, view, STACK_SRC);
-        }
+    if (!psphotModelBackground (config, view, STACK_RAW)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "BACKMDL")) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
+	return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+
+// XXX TEST for background:
+    if (!psphotModelBackground (config, view, STACK_RAW)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!psphotModelBackground (config, view, STACK_RAW)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+// XXX TEST END 
 
 #ifdef MAKE_CHISQ_IMAGE
     // also make the chisq detection image
-    if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
+    if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!strcasecmp (breakPt, "CHISQ")) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 #endif
-    if (!strcasecmp (breakPt, "CHISQ")) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
 
     // find the detections (by peak and/or footprint) in the image.
     // This finds the detections on Chisq image as well as the individuals
-    if (!psphotFindDetections (config, view, STACK_DET, true)) { // pass 1
-        // this only happens if we had an error in psphotFindDetections
+    if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1
         psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    // If DET and SRC are different images, copy the detections from DET to SRC.  This 'copy'
-    // is just a copy of the container pointer; the sources on both DET and SRC are the same
-    // memory objects
-    if (strcmp(STACK_SRC, STACK_DET)) {
-	if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
-	    psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-	    return psphotReadoutCleanup (config, view, STACK_SRC);
-	}
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // construct sources and measure basic stats (saved on detections->newSources)
-    if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1
+    if (!psphotSourceStats (config, view, STACK_RAW, true)) { // pass 1
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "PEAKS")) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    // psphotDumpTest (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    // psphotDumpTest (config, view, STACK_RAW);
     psMemDump("sourcestats");
     logMemStats("sourcestats");
@@ -164,60 +172,48 @@
     // classify sources based on moments, brightness
     // only run this on detections from the input images, not chisq image
-    if (!psphotRoughClass (config, view, STACK_SRC)) {
+    if (!psphotRoughClass (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    // 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)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+
+    // find and subtract radial profile models for saturated stars (XXX change name eventually)
+    if (!psphotDeblendSatstars (config, view, STACK_RAW)) {
 	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);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     // only run this on detections from the input images, not chisq image
-    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
+    if (!psphotImageQuality (config, view, STACK_RAW)) { // pass 1
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "MOMENTS")) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // use bright stellar objects to measure PSF
-    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
+    if (!psphotChoosePSF (config, view, STACK_RAW, true)) { // pass 1
         psLogMsg ("psphot", 3, "failure to construct a psf model");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "PSFMODEL")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // merge the newly selected sources into the existing list
     // NOTE: merge OLD and NEW
-    psphotMergeSources (config, view, STACK_SRC);
+    psphotMergeSources (config, view, STACK_RAW);
 
     // Construct an initial model for each object, set the radius to fitRadius, set circular
     // fit mask.  NOTE: only applied to sources without guess models
-    psphotGuessModels (config, view, STACK_SRC);
+    psphotGuessModels (config, view, STACK_RAW);
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    psphotFitSourcesLinear (config, view, STACK_SRC, false, false);
-    psphotStackVisualFilerule(config, view, STACK_SRC);
+    psphotFitSourcesLinear (config, view, STACK_RAW, false, false);
+    psphotStackVisualFilerule(config, view, STACK_RAW);
 
     // measure the radial profiles to the sky
-    psphotRadialProfileWings (config, view, STACK_SRC);
+    psphotRadialProfileWings (config, view, STACK_RAW);
 
     // re-measure the kron mags with models subtracted.  this pass starts with a circular
@@ -225,16 +221,16 @@
     // but iterates to an appropriately larger size
     logMemStats("before.kron.1");
-    psphotKronIterate(config, view, STACK_SRC, 1);
+    psphotKronIterate(config, view, STACK_RAW, 1);
     logMemStats("after.kron.1");
 	
     // identify CRs and extended sources
-    psphotSourceSize (config, view, STACK_SRC, true);
+    psphotSourceSize (config, view, STACK_RAW, true);
 
     // non-linear PSF and EXT fit to brighter sources
     // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
-    psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotBlendFit (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // replace all sources (do NOT ignore subtraction state)
-    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
+    psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 1 (detections->allSources)
 
     logMemStats("pass1");
@@ -245,6 +241,6 @@
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
-    // NOTE 2 : this function subtracts the models from the given filerule (SRC), not DET
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 2 (detections->allSources)
+    // NOTE 2 : this function subtracts the models from the given filerule
+    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 2 (detections->allSources)
 
     // NOTE: possibly re-measure background model here with objects subtracted / or masked
@@ -252,22 +248,10 @@
     // NOTE: this block performs the 2nd pass low-significance PSF detection stage
     { 
-	// if DET and SRC are different images, generate children sources for all sources in
-	// the SRC image.  This operation replaces the existing DETECTION container on DET
-	// which is currently a view to the one on SRC).  children sources go to
-	// det->allSources
-	if (strcmp(STACK_SRC, STACK_DET)) {
-	    psphotSourceChildren (config, view, STACK_DET, STACK_SRC); 
-
-	    //  subtract all sources from DET (this will subtract using the psf model for SRC, which
-	    //  will somewhat oversubtract the sources -- this is OK
-	    psphotRemoveAllSources (config, view, STACK_DET, false); // do not ignore subtraction state for sources
-	}
-
 	// add noise for subtracted objects
-	psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+	psphotAddNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
 	// find fainter sources
 	// NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
-	psphotFindDetections (config, view, STACK_DET, false); // pass 2 (detections->peaks, detections->footprints)
+	psphotFindDetections (config, view, STACK_RAW, false); // pass 2 (detections->peaks, detections->footprints)
 
 	// remove noise for subtracted objects (ie, return to normal noise level)
@@ -276,43 +260,30 @@
         bool footprintsUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
         if (!footprintsUseUnsubtracted) {
-    	    psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+    	    psphotSubNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
         }
-
-	// if DET and SRC are different images, copy the detections from DET to SRC 
-	// (this operation just ensures the metadata container has a view on SRC as well
-	if (strcmp(STACK_SRC, STACK_DET)) {
-	    // replace all sources in DET
-	    psphotReplaceAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
-
-	    // copy the newly detected peaks from DET to SRC so SourceStats below can operate on them
-	    if (!psphotCopyPeaks (config, view, STACK_SRC, STACK_DET)) {
-		psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-		return psphotReadoutCleanup (config, view, STACK_SRC);
-	    }
-	}
 
 	// define new sources based on only the new peaks & measure moments
 	// NOTE: new sources are saved on detections->newSources
-	psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)
+	psphotSourceStats (config, view, STACK_RAW, false); // pass 2 (detections->newSources)
 
 	// set source type
 	// NOTE: apply only to detections->newSources
-	if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
+	if (!psphotRoughClass (config, view, STACK_RAW)) { // pass 2 (detections->newSources)
 	    psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
-	    return psphotReadoutCleanup (config, view, STACK_SRC);
+	    return psphotReadoutCleanup (config, view, STACK_RAW);
 	}
 
 	// replace all sources so fit below applies to all at once
 	// NOTE: apply only to OLD sources (which have been subtracted)
-	psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2
+	psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 2
 
 	// merge the newly selected sources into the existing list
 	// NOTE: merge OLD and NEW
 	// XXX check on free of sources...
-	psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
+	psphotMergeSources (config, view, STACK_RAW); // (detections->newSources + detections->allSources -> detections->allSources)
 
 	// Construct an initial model for each object, set the radius to fitRadius, set circular
 	// fit mask.  NOTE: only applied to sources without guess models
-	psphotGuessModels (config, view, STACK_SRC);
+	psphotGuessModels (config, view, STACK_RAW);
     }
 
@@ -325,7 +296,7 @@
     if (splitLinearFit) {
         psLogMsg ("psphot", 3, "splitting fit of detected and matched soures\n");
-        // Fit the detected sources separately from matched that wea are about to create.
+        // Fit the detected sources separately from matched ones that we are about to create.
         // NOTE: apply to ALL sources but only include sources with postitive flux in the fit
-        psphotFitSourcesLinear (config, view, STACK_SRC, true, true); // pass 3 (detections->allSources)
+        psphotFitSourcesLinear (config, view, STACK_RAW, true, true); // pass 3 (detections->allSources)
     }
 
@@ -335,5 +306,5 @@
     // this just match the detections for the chisq image, and not bother measuring the source
     // stats in that case...?
-    objects = psphotMatchSources (config, view, STACK_SRC);
+    objects = psphotMatchSources (config, view, STACK_RAW);
     psMemDump("matchsources");
 
@@ -344,26 +315,26 @@
     // Construct an initial model for each object, set the radius to fitRadius, set circular
     // fit mask.  NOTE: only applied to sources without guess models
-    psphotGuessModels (config, view, STACK_SRC);
+    psphotGuessModels (config, view, STACK_RAW);
 
     psphotStackObjectsUnifyPosition (objects);
 
-    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
+    psphotStackObjectsSelectForAnalysis (config, view, STACK_RAW, objects);
 
     // 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)
+    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 4 (detections->allSources)
 
     // measure the radial profiles to the sky (only measures new objects)
-    psphotRadialProfileWings (config, view, STACK_SRC);
+    psphotRadialProfileWings (config, view, STACK_RAW);
 
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, STACK_SRC);
     logMemStats("before.kron.2");
-    psphotKronIterate(config, view, STACK_SRC, 2);
+    psphotKronIterate(config, view, STACK_RAW, 2);
     logMemStats("after.kron.2");
 
     // measure source size for the remaining sources
     // NOTE: applies only to NEW (unmeasured) sources
-    psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)
+    psphotSourceSize (config, view, STACK_RAW, false); // pass 2 (detections->allSources)
 
     psMemDump("psfstats");
@@ -371,22 +342,27 @@
     // drop matched sources without any useful measurements and set kron radii for the ones
     // we decide to keep
-    psphotFilterMatchedSources (config, view, STACK_SRC, objects);
+    psphotFilterMatchedSources (config, view, STACK_RAW, objects);
 
     // measure kron fluxes for the matched sources only
-    psphotKronIterate(config, view, STACK_SRC, 3);
+    psphotKronIterate(config, view, STACK_RAW, 3);
 
     // measure elliptical apertures, petrosians (objects sorted by S/N)
     // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
-    psphotExtendedSourceAnalysis (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotExtendedSourceAnalysis (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
-    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotExtendedSourceFits (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // create source children for the OUT filerule (for radial aperture photometry and output) 
-    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
+    // NOTE: The new source children have image arrays pointing to the readout associated with 
+    // STACK_OUT.  in psphotStackMatchPSFsetup, we copy the current pixel values from RAW to OUT, 
+    // but keep the pointers the same so we do not break these source image references
+    // XXX NOTE : if we use the pre-20130914 psphotStackReadout code, we need to use 'false' for the
+    // sourcesSubtracted argument
+    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects, true);
     if (!objectsOut) {
 	psFree(objects);
 	psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
@@ -396,98 +372,82 @@
         // this forces photometry on the undetected sources from other images
 
-        // NOTE: we always do the radial apertures analysis on the convolved image since
-        // those are the ones that are psf matched and are the source of STACK_OUT's pixels
-        // XXX: Actually if PSPHOT.STACK.MATCH.PSF.SOURCE were set to RAW this wouldn't be true.
-        // but in that case we don't get past the psf matching step because there is no
-        // target psf for the RAW inputs
-
-        // If useRaw copy the sources to the convolved readout
-        if (strcmp(STACK_SRC, STACK_CNV)) {
-            if (!psphotCopySources (config, view, STACK_CNV, STACK_SRC)) {
-                psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-                return psphotReadoutCleanup (config, view, STACK_SRC);
-            }
-        }
-        // 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");
-
+	// set up the FWHM vector
+	psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW);
+	psphotDumpImages (config, view, STACK_RAW, "raw.t0");
+	psphotDumpImages (config, view, STACK_OUT, "out.t0");
 
         int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
-        for (int entry = 1; entry < nRadialEntries; entry++) {
+
+        for (int entry = 0; entry < nRadialEntries; entry++) {
             // NOTE: entry 0 is the unmatched image set
 
-            // re-measure the PSF for the smoothed image (using entries in 'allSources')
-            psphotChoosePSF (config, view, STACK_OUT, false);
-
-            // this is necessary to update the models based on the new PSF
-            psphotResetModels (config, view, STACK_OUT);
-
-            // this is necessary to get the right normalization for the new models
-            psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
+	    char line[256];
 
             // measure circular, radial apertures (objects sorted by S/N)
             psphotRadialApertures (config, view, STACK_OUT, entry); 
+	    snprintf (line, 256, "%s.%d", "out.t1", entry);
+	    psphotDumpImages (config, view, STACK_OUT, line);
 
             // replace the flux in the image so it is returned to its original state
             psphotReplaceAllSources (config, view, STACK_OUT, false);
-
-            // smooth to the next FWHM, or set 'smoothAgain' to false if no more 
-            psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
-            psMemDump("matched");
+	    snprintf (line, 256, "%s.%d", "out.t2", entry);
+	    psphotDumpImages (config, view, STACK_OUT, line);
+
+	    if (entry < nRadialEntries - 1) {
+		// smooth to the next FWHM
+		// this function does nothing if the targetFWHM is smaller than the currentFWHM
+		psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
+		snprintf (line, 256, "%s.%d", "out.t3", entry);
+		psphotDumpImages (config, view, STACK_OUT, line);
+		psMemDump("matched");
+
+		// re-measure the PSF for the smoothed image (using entries in 'allSources')
+		psphotChoosePSF (config, view, STACK_OUT, false);
+
+		// this is necessary to update the models based on the new PSF
+		psphotResetModels (config, view, STACK_OUT);
+
+		// this is necessary to get the right normalization for the new models
+		// and to subtract the sources
+		psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
+		snprintf (line, 256, "%s.%d", "out.t4", entry);
+		psphotDumpImages (config, view, STACK_OUT, line);
+	    }
         }
     }
-    psphotStackSetInputsToSkip(config, view, STACK_CNV, false);
-    psphotStackSetInputsToSkip(config, view, STACK_OUT, false);
 
     // measure aperture photometry corrections
-    if (!psphotApResid (config, view, STACK_SRC)) {
+    if (!psphotApResid (config, view, STACK_RAW)) {
 	psFree (objects);
 	psFree (objectsOut);
         psLogMsg ("psphot", 3, "failed on psphotApResid");
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // calculate source magnitudes
-    psphotMagnitudes(config, view, STACK_SRC);
-
-    if (!useRaw) {
-        // psphotEfficiency wants to have the PSF of the image, but since we are measuring on
-        // the convolved images we need to generate PSFs for the DET images
-        if (!psphotChoosePSF (config, view, STACK_DET, false)) {
-            psLogMsg ("psphot", 3, "failure to construct a psf model for raw input");
-            return psphotReadoutCleanup (config, view, STACK_DET);
-        }
-    }
-    if (!psphotEfficiency(config, view, STACK_DET)) {
+    psphotMagnitudes(config, view, STACK_RAW);
+
+    if (!psphotEfficiency(config, view, STACK_RAW)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();
     }
-    psphotCopyEfficiency (config, view, STACK_OUT, STACK_DET);
+    psphotCopyEfficiency (config, view, STACK_OUT, STACK_RAW);
 
     logMemStats("final");
 #if (1)
-    psphotSourceMemory(config, view, STACK_SRC);
+    psphotSourceMemory(config, view, STACK_RAW);
     psphotSourceMemory(config, view, STACK_OUT);
 #endif
 
-    // replace failed sources?
-    // psphotReplaceUnfitSources (sources);
-
     // replace background in residual image
-    psphotSkyReplace (config, view, STACK_DET);
+    psphotSkyReplace (config, view, STACK_RAW);
 
     // drop the references to the image pixels held by each source
+    psphotSourceFreePixels (config, view, STACK_RAW);
     psphotSourceFreePixels (config, view, STACK_OUT);
-    psphotSourceFreePixels (config, view, STACK_SRC);
 
 #ifdef MAKE_CHISQ_IMAGE
     // remove chisq image from config->file:PSPHOT.INPUT
-    psphotStackRemoveChisqFromInputs(config, STACK_DET);
-    if (strcmp(STACK_SRC, STACK_DET)) {
-	psphotStackRemoveChisqFromInputs(config, STACK_SRC);
-    }
+    psphotStackRemoveChisqFromInputs(config, STACK_RAW);
 #endif
 
@@ -496,5 +456,5 @@
 
     // create the exported-metadata and free local data
-    return psphotReadoutCleanup (config, view, STACK_SRC);
+    return psphotReadoutCleanup (config, view, STACK_RAW);
 }
 
@@ -552,11 +512,8 @@
 }
 
-
-
 /* here is the process:
 
- * we have three(*) images:
- * RAW : unconvolved image stack 
- * CNV : input convolved image
+ * we have two image sets:
+ * RAW : unconvolved image stacks
 
  * OUT : psf-matched output image (there may be more than one of
@@ -649,2 +606,124 @@
 
    */
+
+
+// generate a vector fwhmValues where the first has the fwhm of the raw image and the
+// successive entries have the target values
+
+bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *filerulePSF) {
+
+    bool status;
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+
+	if (!psphotStackMatchPSFsetupReadout (config, view, filerule, filerulePSF, i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to define target PSF sizes");
+	    return false;
+	}
+    }
+
+    return true;
+}
+
+float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index);
+
+// copy the pixels from RAW to OUT (
+
+bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *fileruleOut, const char *fileruleRaw, int index) {
+
+    bool status;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+
+    // find the currently selected readout
+    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, fileruleOut, index); // File of interest
+    psAssert (fileOut, "missing file?");
+
+    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
+    psAssert (readoutOut, "missing readout?");
+
+    // find the currently selected readout
+    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, fileruleRaw, index); // File of interest
+    psAssert (fileRaw, "missing file?");
+
+    pmReadout *readoutRaw = pmFPAviewThisReadout(view, fileRaw->fpa);
+    psAssert (readoutRaw, "missing readout?");
+
+    readoutOut->image = psImageCopy(readoutOut->image, readoutRaw->image, PS_TYPE_F32);
+    if (readoutRaw->variance) {
+	readoutOut->variance = psImageCopy(readoutOut->variance, readoutRaw->variance, PS_TYPE_F32);
+    }
+    if (readoutRaw->mask) {
+	readoutOut->mask = psImageCopy(readoutOut->mask, readoutRaw->mask, PS_TYPE_IMAGE_MASK);
+    }
+
+    // pmChip *chipRaw = pmFPAviewThisChip(view, fileRaw->fpa); // The chip holds the PSF
+    // psAssert (chipRaw, "missing chip");
+
+    pmPSF *psf = psMetadataLookupPtr(&status, readoutRaw->analysis, "PSPHOT.PSF"); // PSF
+    if (!psf) {
+	// we should have a PSF by this point in psphot
+	psError(PSPHOT_ERR_PROG, true, "Unable to find PSF.");
+	return false;
+    }
+
+    float fwhmRaw = psphotPSFseeing (psf, readoutRaw, index);
+
+    psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
+    psVectorAppend(fwhmValues, fwhmRaw);
+
+    // is a single target FWHM specified, or a set of values?  set up the vector options->targetSeeing and the local 1st value
+    float targetSeeing = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
+    if (!status) {
+	psVector *targetSeeing = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
+	psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+	for (int i = 0; i < targetSeeing->n; i++) {
+	    psVectorAppend(fwhmValues, targetSeeing->data.F32[i]);
+	}	    
+    } else {
+        psVectorAppend(fwhmValues, targetSeeing);
+    }
+
+    psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
+    psFree (fwhmValues);
+
+    return true;
+}
+
+float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index) {
+
+    psImage *image = readout->image;
+
+    int Nx = image->numCols;
+    int Ny = image->numRows;
+
+    float sumFWHM = 0.0;		  // FWHM for image
+    int numFWHM = 0;			  // Number of FWHM measurements
+    for (float x = 0; x < Nx; x += 0.25*Nx) {
+	for (float y = 0; y < Ny; y += 0.25*Ny) {
+	    float fwhm = pmPSFtoFWHM(psf, x, y);
+	    if (isfinite(fwhm)) {
+		sumFWHM += fwhm;
+		numFWHM++;
+	    }
+	}
+    }
+    if (numFWHM == 0) {
+	psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", index);
+	return NAN;
+    } 
+
+    float fwhm = sumFWHM / (float) numFWHM;
+    psLogMsg ("psphotStack", PS_LOG_INFO, "Input Seeing for %d: %f\n", index, fwhm);
+
+    return fwhm;
+}