Index: /branches/eam_branches/ipp-20110710/psphot/src/psphot.h
===================================================================
--- /branches/eam_branches/ipp-20110710/psphot/src/psphot.h	(revision 32321)
+++ /branches/eam_branches/ipp-20110710/psphot/src/psphot.h	(revision 32322)
@@ -377,13 +377,4 @@
     bool convolve;                      // Convolve images?
     psphotStackConvolveSource convolveSource;
-    // psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images
-
-    // bool matchZPs;                      // Adjust relative fluxes based on transparency analysis?
-    // bool photometry;                    // Perform photometry?
-    // psMetadata *stats;                  // Statistics for output
-    // FILE *statsFile;                    // File to which to write statistics
-    // psArray *origImages, *origMasks, *origVariances; // Filenames of the original images
-    // psArray *origCovars;                // Original covariances matrices
-    // int quality;                        // Bad data quality flag
 
     // Prepare
@@ -392,14 +383,8 @@
     psVector *inputMask;                // Mask for inputs
 
-    float targetSeeing;                 // Target seeing FWHM
+    psVector *targetSeeing;		// Target seeing FWHMs
     psArray *sourceLists;               // Individual lists of sources for matching
     psVector *norm;                     // Normalisation for each image
     psArray *psfs;
-
-    // psVector *exposures;                // Exposure times
-    // float sumExposure;                  // Sum of exposure times
-    // float zp;                           // Zero point for output
-    // psVector *inputMask;                // Mask for inputs
-    // psArray *sources;                   // Matched sources
 
     // Convolve
@@ -407,16 +392,4 @@
     psArray *regions;                   // PSF-matching regions --- required in the stacking
     psVector *matchChi2;                // chi^2 for stamps from matching
-    psVector *weightings;               // Combination weightings for images (1/noise^2)
-    // psArray *cells;                     // Cells for convolved images --- a handle for reading again
-    // int numCols, numRows;               // Size of image
-    // psArray *convCovars;                // Convolved covariance matrices
-
-    // Combine initial
-    // pmReadout *outRO;                   // Output readout
-    // pmReadout *expRO;                   // Exposure readout
-    // psArray *inspect;                   // Array of arrays of pixels to inspect
-
-    // Rejection
-    // psArray *rejected;                  // Rejected pixels
 } psphotStackOptions;
 
@@ -425,19 +398,23 @@
 bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
 bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index);
-bool psphotStackMatchPSFsNext (bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize);
-bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize);
+
+bool psphotStackMatchPSFsNext (pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize);
+bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, int lastSize);
+int psphotStackMatchPSFsEntries (pmConfig *config, const pmFPAview *view, const char *filerule);
 
 // psphotStackMatchPSFsUtils
-psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
-pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
-bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index);
-bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index);
+// psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
+// pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
+// bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index);
+// bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index);
 bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index);
 bool matchKernel(pmConfig *config, pmReadout *cnv, pmReadout *raw, psphotStackOptions *options, int index);
-bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
-bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
-bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index);
-
-pmPSF *psphotStackPSF(const pmConfig *config, int numCols, int numRows, const psArray *psfs, const psVector *inputMask);
+// bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
+// bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
+// bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index);
+
+bool psphotStackRenormaliseVariance(const pmConfig *config, pmReadout *readout);
+
+bool psphotStackPSF(const pmConfig *config, psphotStackOptions *options);
 
 psphotStackOptions *psphotStackOptionsAlloc (int num);
@@ -448,6 +425,6 @@
 bool psphotCopySourcesReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
 
-bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule);
-bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule, int entry);
+bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, int entry);
 bool psphotRadialApertures_Threaded (psThreadJob *job);
 bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry);
@@ -489,3 +466,5 @@
 bool psphotKronIterate_Threaded (psThreadJob *job);
 
+bool psphotStackObjectsSelectForAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
+
 #endif
Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotRadialApertures.c
===================================================================
--- /branches/eam_branches/ipp-20110710/psphot/src/psphotRadialApertures.c	(revision 32321)
+++ /branches/eam_branches/ipp-20110710/psphot/src/psphotRadialApertures.c	(revision 32322)
@@ -3,8 +3,14 @@
 bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar);
 
+// this function measures the radial aperture fluxes for the set of readouts.  this function
+// may be called multiple times, presumably for different versions of PSF-matched or unmatched images.  
+
 // for now, let's store the detections on the readout->analysis for each readout
-bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule)
+bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule, int entry)
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Radial Apertures ---");
 
     // select the appropriate recipe information
@@ -22,5 +28,5 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-	if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe)) {
+	if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe, entry)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i);
 	    return false;
@@ -32,5 +38,6 @@
 // aperture-like measurements for extended sources
 // flux in simple, circular apertures
-bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+// 'entry' tells us which of the matched-PSF images we are working on (0 == unmatched image, also non-stack psphot)
+bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, int entry) {
 
     bool status;
@@ -72,4 +79,17 @@
     // XXX are we getting the objects out of order? does it matter?
     sources = psArraySort (sources, pmSourceSortByFlux);
+
+    // XXX make this consistent with entry 0 == unmatched
+    int nEntry = 1;
+    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+    if (fwhmValues) {
+	psAssert (entry < fwhmValues->n, "inconsistent matched-PSF entry");
+	nEntry = fwhmValues->n;
+    }
+    if (entry > 0) {
+	psLogMsg ("psphot", PS_LOG_DETAIL, "Radial Apertures for matched image %s : PSF FWHM = %f pixels\n", file->name, fwhmValues->data.F32[entry]);
+    } else {
+	psLogMsg ("psphot", PS_LOG_DETAIL, "Radial Apertures for unmatched image %s\n", file->name);
+    }
 
     // option to limit analysis to a specific region
@@ -98,4 +118,7 @@
             psArrayAdd(job->args, 1, AnalysisRegion);
             psArrayAdd(job->args, 1, recipe);
+
+            PS_ARRAY_ADD_SCALAR(job->args, entry,  PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, nEntry, PS_TYPE_S32);
 
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial
@@ -115,5 +138,5 @@
 	    }
 	    psScalar *scalar = NULL;
-	    scalar = job->args->data[4];
+	    scalar = job->args->data[6];
 	    Nradial += scalar->data.S32;
 	    psFree(job);
@@ -135,5 +158,5 @@
             } else {
 		psScalar *scalar = NULL;
-		scalar = job->args->data[4];
+		scalar = job->args->data[6];
 		Nradial += scalar->data.S32;
             }
@@ -159,4 +182,7 @@
     psMetadata *recipe      = job->args->data[3];
 
+    int entry               = PS_SCALAR_VALUE(job->args->data[4],S32); // which psf-matched image are we working on? (0 == unmatched)
+    int nEntry              = PS_SCALAR_VALUE(job->args->data[5],S32); // total number of psf-matched images + 1 unmatched
+
     // radMax stores the upper bounds of the annuli
     // XXX keep the same name here as for the petrosian / elliptical apertures?
@@ -198,6 +224,13 @@
 
 	// allocate pmSourceExtendedParameters, if not already defined
-	if (!source->radialAper) {
-	    source->radialAper = psArrayAlloc(1);
+	// XXX check that nPSFsizes is consistent with targets
+	if (source->parent) {
+	    if (!source->parent->radialAper) {
+		source->parent->radialAper = psArrayAlloc(nEntry);
+	    }
+	} else {
+	    if (!source->radialAper) {
+		source->radialAper = psArrayAlloc(nEntry);
+	    }
 	}
 
@@ -219,5 +252,5 @@
 	pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
 
-	if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, 0)) {
+	if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, entry)) {
 	    psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
 	} else {
@@ -233,5 +266,5 @@
 	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     }
-    psScalar *scalar = job->args->data[4];
+    psScalar *scalar = job->args->data[6];
     scalar->data.S32 = Nradial;
 
Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotReadout.c
===================================================================
--- /branches/eam_branches/ipp-20110710/psphot/src/psphotReadout.c	(revision 32321)
+++ /branches/eam_branches/ipp-20110710/psphot/src/psphotReadout.c	(revision 32322)
@@ -141,4 +141,6 @@
 
     // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
+    // XXX I've deactivated this because it was preventing galaxies close to stars from being
+    // XXX fitted as an extended source.
     if (false && !psphotBasicDeblend (config, view, filerule)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
@@ -313,5 +315,5 @@
     psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources)
     psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources)
-    psphotRadialApertures(config, view, filerule);
+    psphotRadialApertures(config, view, filerule, 0);
 
 finish:
Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotSetThreads.c
===================================================================
--- /branches/eam_branches/ipp-20110710/psphot/src/psphotSetThreads.c	(revision 32321)
+++ /branches/eam_branches/ipp-20110710/psphot/src/psphotSetThreads.c	(revision 32322)
@@ -50,5 +50,5 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 5);
+    task = psThreadTaskAlloc("PSPHOT_RADIAL_APERTURES", 7);
     task->function = &psphotRadialApertures_Threaded;
     psThreadTaskAdd(task);
Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceMatch.c
===================================================================
--- /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceMatch.c	(revision 32321)
+++ /branches/eam_branches/ipp-20110710/psphot/src/psphotSourceMatch.c	(revision 32322)
@@ -48,11 +48,7 @@
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     psAssert (detections, "missing detections?");
-    psAssert (detections->newSources, "new sources not defined?");
-    psAssert (!detections->allSources, "all sources already defined?");
-
-    // XXX TEST:
-    if (detections->newSources) {
-        psphotMatchSourcesToObjects(objects, detections->newSources, RADIUS);
-    }
+    psAssert (detections->allSources, "all sources not defined?");
+
+    psphotMatchSourcesToObjects(objects, detections->allSources, RADIUS);
 
     return true;
@@ -261,5 +257,5 @@
 	    peak->assigned = true;
 	    pmPhotObjAddSource(obj, source);
-	    psArrayAdd (detections->newSources, 100, source);
+	    psArrayAdd (detections->allSources, 100, source);
 	    psFree (source);
 	}
Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotStackMatchPSFsNext.c
===================================================================
--- /branches/eam_branches/ipp-20110710/psphot/src/psphotStackMatchPSFsNext.c	(revision 32321)
+++ /branches/eam_branches/ipp-20110710/psphot/src/psphotStackMatchPSFsNext.c	(revision 32322)
@@ -1,3 +1,26 @@
 # include "psphotInternal.h"
+
+// NOTE : element 0 of fwhmValues if the unmatched image,  
+
+int psphotStackMatchPSFsEntries (pmConfig *config, const pmFPAview *view, const char *filerule) {
+
+    int nRadialEntries = 0;
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest
+    psAssert (file, "missing file?");
+    
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+    
+    bool status = false;
+    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+    if (!fwhmValues) {
+	return 1;
+    }
+    
+    nRadialEntries = fwhmValues->n;
+    return nRadialEntries;
+}
 
 // smooth the input image to match the next target PSF
@@ -5,25 +28,6 @@
 // and that the smoothing can use a 1D Gaussian kernel of width sqrt(TARGET^2 - CURRENT^2)
 // each subsequent call
-bool psphotStackMatchPSFsNext(bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize)
+bool psphotStackMatchPSFsNext(pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize)
 {
-    bool status = true;
-
-    // select the appropriate recipe information
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
-    psAssert (recipe, "missing recipe?");
-
-    psVector *fwhmValues = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
-    if (!status) {
-	// must not be a vector, only one value requested
-	*smoothAgain = false;
-	return true;
-    }
-
-    if (lastSize + 1 >= fwhmValues->n) {
-	// all done with target FWHM values 
-	*smoothAgain = false;
-	return true;
-    }
-
     int num = psphotFileruleCount(config, filerule);
 
@@ -33,6 +37,6 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-	if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, recipe, fwhmValues, lastSize)) {
-            psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF %f", filerule, i, fwhmValues->data.F32[lastSize+1]);
+	if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, lastSize)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF", filerule, i);
 	    psImageConvolveSetThreads(oldThreads);
 	    return false;
@@ -44,5 +48,5 @@
 }
 
-bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize) {
+bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, int lastSize) {
 
     bool status = false;
@@ -58,4 +62,8 @@
     psphotVisualShowImage(readout);
 
+    // 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)
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
@@ -66,4 +74,11 @@
         psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value");
         minGauss = 0.5;
+    }
+
+    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+    psAssert (fwhmValues, "need target PSFs");
+
+    if (lastSize + 1 >= fwhmValues->n) {
+	return true;
     }
 
@@ -128,11 +143,4 @@
     // psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor);
 
-    // save the output fwhm values in the readout->analysis.  we may have / will have multiple output PSF sizes,
-    // so we save this in a vector.  if the vector is not yet defined, create it
-
-    psVector *fwhmValuesOut = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
-    psAssert(fwhmValuesOut, "should already exist..");
-    psVectorAppend(fwhmValuesOut, targetFWHM);
-
     // do not generate a PSF if we already were supplied one
     pmPSF *psfOld = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotStackReadout.c
===================================================================
--- /branches/eam_branches/ipp-20110710/psphot/src/psphotStackReadout.c	(revision 32321)
+++ /branches/eam_branches/ipp-20110710/psphot/src/psphotStackReadout.c	(revision 32322)
@@ -44,4 +44,6 @@
 bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
 
+    // measure the total elapsed time in psphotReadout.  dtime is the elapsed time used jointly
+    // by the multiple threads, not the total time used by all threads.
     psTimerStart ("psphotReadout");
 
@@ -56,12 +58,11 @@
     // optional break-point for processing
     char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
-    PS_ASSERT_PTR_NON_NULL (breakPt, false);
-
-    // XXX or do I set OUT to be a pmFPAfile pointing at the input of interest?
+    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; // XXX optionally allow this to be CNV?
-
-    // we have 3 relevant files: RAW, CNV, OUT 
+    char *STACK_DET = STACK_RAW;
 
     // set the photcode for each image
@@ -71,6 +72,5 @@
     }
 
-    // Generate the mask and weight images
-    // XXX this should be done before we perform the convolutions
+    // Generate the mask and weight images (if not supplied) and set mask bits
     if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
 	return psphotReadoutCleanup (config, view, STACK_SRC);
@@ -80,7 +80,9 @@
     }
 
+    // XXX I think this is not defined correctly for an array of images.
+    // XXX I probably need to subtract the model (same model?) for both RAW and OUT.
+    // XXX But, probably not a problem in practice since the stacks are constructed with 0.0 mean level.
+
     // generate a background model (median, smoothed image)
-    // XXX I think this is not defined correctly for an array of images.
-    // XXX probably need to subtract the model (same model?) for both RAW and OUT
     if (!psphotModelBackground (config, view, STACK_DET)) {
 	return psphotReadoutCleanup (config, view, STACK_SRC);
@@ -93,4 +95,5 @@
     }
 
+    // also make the chisq detection image
     if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
@@ -109,5 +112,5 @@
     }
 
-    // copy the detections from DET to SRC 
+    // if DET and SRC are different images, copy the detections from DET to SRC 
     if (strcmp(STACK_SRC, STACK_DET)) {
 	if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
@@ -122,111 +125,150 @@
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
-
-    if (!strcasecmp (breakPt, "TEST1")) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
+    if (!strcasecmp (breakPt, "PEAKS")) {
+	return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
     psMemDump("sourcestats");
 
-    // generate the objects (object unify the sources from the different images)
+    // classify sources based on moments, brightness
+    // only run this on detections from the input images, not chisq image
+    if (!psphotRoughClass (config, view, STACK_SRC)) {
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
+	return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
+    // only run this on detections from the input images, not chisq image
+    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    if (!strcasecmp (breakPt, "MOMENTS")) {
+	return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    // use bright stellar objects to measure PSF
+    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
+        psLogMsg ("psphot", 3, "failure to construct a psf model");
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    if (!strcasecmp (breakPt, "PSFMODEL")) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
+    psphotGuessModels (config, view, STACK_SRC);
+
+    // merge the newly selected sources into the existing list
+    // NOTE: merge OLD and NEW
+    psphotMergeSources (config, view, STACK_SRC);
+
+    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
+    // XXX why do this as a stack operation?
+    // psphotFitSourcesLinearStack (config, objects, false);
+    psphotFitSourcesLinear (config, view, STACK_SRC, false);
+    psphotStackVisualFilerule(config, view, STACK_SRC);
+
+    // re-measure the kron mags with models subtracted.  this pass starts with a circular
+    // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
+    // but iterates to an appropriately larger size
+    psphotKronIterate(config, view, STACK_SRC);
+
+    // identify CRs and extended sources
+    psphotSourceSize (config, view, STACK_SRC, 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)
+
+    // replace all sources
+    psphotReplaceAllSources (config, view, STACK_SRC); // pass 1 (detections->allSources)
+
+    // linear fit to include all sources (subtract again)
+    // NOTE : apply to ALL sources (extended + psf)
+    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)
+
+    // if we only do one pass, skip to extended source analysis
+    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
+
+    // NOTE: possibly re-measure background model here with objects subtracted / or masked
+
+    // NOTE: this block performs the 2nd pass low-significance PSF detection stage
+    { 
+	// add noise for subtracted objects
+	psphotAddNoise (config, view, STACK_DET); // 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)
+
+	// remove noise for subtracted objects (ie, return to normal noise level)
+	// NOTE: this needs to operate only on the OLD sources
+	psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+
+	// if DET and SRC are different images, copy the detections from DET to SRC 
+	if (strcmp(STACK_SRC, STACK_DET)) {
+	    // XXX how does this handle 1st vs 2nd pass sources?
+	    if (!psphotCopySources (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)
+
+	// set source type
+	// NOTE: apply only to detections->newSources
+	if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
+	    psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
+	    return psphotReadoutCleanup (config, view, STACK_SRC);
+	}
+
+	// create full input models, set the radius to fitRadius, set circular fit mask
+	// NOTE: apply only to detections->newSources
+	psphotGuessModels (config, view, STACK_SRC); // pass 2 (detections->newSources)
+
+	// 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); // 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)
+
+	// NOTE: apply to ALL sources
+	psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
+    }
+
+pass1finish:
+
+    // re-measure the kron mags with models subtracted
+    // psphotKronMasked(config, view, STACK_SRC);
+    psphotKronIterate(config, view, STACK_SRC);
+
+    // measure source size for the remaining sources
+    // NOTE: applies only to NEW (unmeasured) sources
+    psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)
+
+    psMemDump("psfstats");
+
+    // generate the objects (objects unify the sources from the different images)
     // XXX this could just match the detections for the chisq image, and not bother measuring the
     // source stats in that case...
     psArray *objects = psphotMatchSources (config, view, STACK_SRC);
-
     psMemDump("matchsources");
 
-    if (!strcasecmp (breakPt, "TEST2")) {
-	psFree(objects);
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    // construct sources for the newly-generated sources (from other images)
-    if (!psphotSourceStats (config, view, STACK_SRC, false)) { // pass 1
-	psFree(objects);
-        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    psMemDump("sourcestats");
-
-    // find blended neighbors of very saturated stars (detections->newSources)
-    // if (!psphotDeblendSatstars (config, view)) {
-    //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-    //     return psphotReadoutCleanup (config, view, STACK_SRC);
-    // }
-
-    // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
-    // if (!psphotBasicDeblend (config, view)) {
-    //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
-    //     return psphotReadoutCleanup (config, view, STACK_SRC);
-    // }
-
-    // classify sources based on moments, brightness
-    // only run this on detections from the input images, not chisq image
-    if (!psphotRoughClass (config, view, STACK_SRC)) {
-	psFree(objects);
-        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
-    // only run this on detections from the input images, not chisq image
-    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
-	psFree(objects);
-        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (!strcasecmp (breakPt, "MOMENTS")) {
-	psFree(objects);
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
-    // this step is skipped
-    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
-	psFree(objects);
-        psLogMsg ("psphot", 3, "failure to construct a psf model");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (!strcasecmp (breakPt, "PSFMODEL")) {
-	psFree(objects);
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
-    psphotGuessModels (config, view, STACK_SRC);
-
-    // merge the newly selected sources into the existing list
-    // NOTE: merge OLD and NEW
-    psphotMergeSources (config, view, STACK_SRC);
-
-    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    psphotFitSourcesLinearStack (config, objects, FALSE);
-    psphotStackVisualFilerule(config, view, STACK_SRC);
-
-    // identify CRs and extended sources
-    psphotSourceSize (config, view, STACK_SRC, TRUE);
-
-    // XXX do we want to do a preliminary (unconvolved) model fit here, and then 
-    // do a second detection pass? (like standard psphot)
-
-    // measure aperture photometry corrections
-    if (!psphotApResid (config, view, STACK_SRC)) {
-	psFree (objects);
-        psLogMsg ("psphot", 3, "failed on psphotApResid");
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    psMemDump("psfstats");
-
     psphotStackObjectsUnifyPosition (objects);
 
+    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
+
     // measure elliptical apertures, petrosians (objects sorted by S/N)
-    psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
+    // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotExtendedSourceAnalysis (config, view, STACK_SRC); // pass 1 (detections->allSources)
 
     // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
     psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
-
-    // calculate source magnitudes
-    psphotMagnitudes(config, view, STACK_SRC);
 
     // create source children for the OUT filerule (for radial aperture photometry) 
@@ -238,8 +280,15 @@
     }
 
+    // measure circular, radial apertures (objects sorted by S/N)
+    // XXX can we just use psphotRadialApertures
+    // XXX make sure the headers are consistent with this (which PSF convolutions, ie mark 'none')
+    // psphotRadialAperturesByObject (config, objectsRadial, view, STACK_SRC, nMatchedPSF); 
+    psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
     psMemDump("extmeas");
 
-    bool smoothAgain = true;
-    for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) {
+    int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
+
+    for (int entry = 1; entry < nRadialEntries; entry++) {
+	// NOTE: entry 0 is the unmatched image set
 
 	// re-measure the PSF for the smoothed image (using entries in 'allSources')
@@ -253,5 +302,6 @@
 
 	// measure circular, radial apertures (objects sorted by S/N)
-	psphotRadialAperturesByObject (config, objectsRadial, view, STACK_OUT, nMatchedPSF); 
+	// entry 0 == unmatched? pass entry + 1?
+	psphotRadialApertures (config, view, STACK_OUT, entry); 
 
 	// replace the flux in the image so it is returned to its original state
@@ -259,9 +309,19 @@
 
 	// smooth to the next FWHM, or set 'smoothAgain' to false if no more 
-	psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF);
+	psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
 	psMemDump("matched");
     }
 
-    if (0 && !psphotEfficiency(config, view, STACK_OUT)) {
+    // measure aperture photometry corrections
+    if (!psphotApResid (config, view, STACK_SRC)) {
+	psFree (objects);
+        psLogMsg ("psphot", 3, "failed on psphotApResid");
+	return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
+    // calculate source magnitudes
+    psphotMagnitudes(config, view, STACK_SRC);
+
+    if (0 && !psphotEfficiency(config, view, STACK_DET)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();