Index: /trunk/psphot/src/psphotReadout.c
===================================================================
--- /trunk/psphot/src/psphotReadout.c	(revision 11177)
+++ /trunk/psphot/src/psphotReadout.c	(revision 11178)
@@ -12,6 +12,5 @@
     char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
 
-    // XXX does this need to invoke I/O?
-    // XXX move this input the psphotImageLoop level?
+    // generate mask & weight images if they don't already exit
     if (!pmReadoutGenerateMaskWeight (readout, true)) {
         return false;
@@ -99,28 +98,39 @@
     psphotBlendFit (readout, sources, recipe, psf);
 
-    // next steps:
-    // remeasure the sky after objects have been subrtracted
-    // re-find the peaks.
-    // re-run ensemble fit
-
-    // XXX this was an attempt to measure bias due to the weighting.
-    // psphotWeightBias (readout, sources, recipe, psf);
-
-    // replace fitted sources
-    // XXX we need to replace failed objects, not all objects
-    // psphotReplaceUnfit (sources);
-
-    // XXX find remaining peaks after first source subtraction pass
-    // psphotFindPeaks ();
-
-    // replace subtracted sources
+    // replace all sources
     psphotReplaceAll (sources);
-    // psphotSaveImage (NULL, readout->image,  "replace.fits");
 
     // linear PSF fit to remaining peaks
     psphotEnsemblePSF (readout, sources, recipe, psf, TRUE);
+    if (!strcasecmp (breakPt, "PASS1")) {
+        goto finish;
+    }
 
-    // XXX test the plotting function
-    // psphotPlotPSFModel (config, view, sources);
+    // replace background in residual image
+    psphotSkyReplace (config, view);
+
+    // re-measure background model (median, smoothed image)
+    psphotImageMedian (config, view);
+
+    // find the peaks in the image
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", 5.0);
+    psArray *newPeaks = psphotFindPeaks (readout, recipe);
+
+    // define new sources based on the new peaks
+    psArray *newSources = psphotSourceStats (readout, recipe, newPeaks);
+    psFree (newPeaks);
+
+    // set source type
+    psphotRoughClass (newSources, recipe);
+
+    // create full input models
+    psphotGuessModels (readout, newSources, recipe, psf);
+
+    // linear PSF fit to remaining peaks
+    psphotEnsemblePSF (readout, newSources, recipe, psf, TRUE);
+
+    // merge the newly selected peaks into the existing list
+    // psphotIncludeNewPeaks (sources, newSources);
+    psFree (newSources);
 
 finish:
@@ -135,4 +145,8 @@
     // calculate source magnitudes
     psphotMagnitudes(sources, recipe, psf, background);
+
+    // replace fitted sources
+    // XXX we need to replace failed objects, not all objects
+    // psphotReplaceUnfit (sources);
 
     // replace background in residual image