Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotMergeSources.c
===================================================================
--- branches/eam_branches/stackphot.20100406/psphot/src/psphotMergeSources.c	(revision 27649)
+++ branches/eam_branches/stackphot.20100406/psphot/src/psphotMergeSources.c	(revision 27650)
@@ -76,7 +76,8 @@
     pmDetections *extCMF = NULL;
     psArray *extSourcesTXT = NULL;
-
-    // find the currently selected readout
-    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest
+    int index = 0;
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest
     psAssert (file, "missing file?");
 
@@ -118,4 +119,5 @@
                 psFree (source->modelPSF);
                 source->modelPSF = NULL;
+		source->imageID = index;
 
                 psArrayAdd (detections->newSources, 100, source);
@@ -140,4 +142,5 @@
                 psFree (source->modelPSF);
                 source->modelPSF = NULL;
+		source->imageID = index;
 
                 psArrayAdd (detections->newSources, 100, source);
Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceMatch.c
===================================================================
--- branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceMatch.c	(revision 27649)
+++ branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceMatch.c	(revision 27650)
@@ -1,3 +1,5 @@
 # include "psphotInternal.h" 
+
+bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects);
  
 psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view) 
@@ -19,5 +21,5 @@
     }
 
-    // psphotMatchSourcesGenerate (objects, config, view, "PSPHOT.INPUT", i);
+    psphotMatchSourcesGenerate (config, view, objects);
 
     return objects;
@@ -135,7 +137,120 @@
 	pmPhotObjAddSource(obj, src);
 	psArrayAdd (objects, 100, obj);
+	psFree (obj);
     }
     psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);
 
+    psFree (found);
     return true;
 } 
+
+bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects) {
+
+    bool status = false;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    // determine properties (sky, moments) of initial sources
+    float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
+    psAssert (status, "missing SKY_OUTER_RADIUS in recipe?");
+
+    int nImages = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // generate look-up arrays for detections and readouts
+    psArray *detArrays = psArrayAlloc(nImages);
+    psArray *readouts = psArrayAlloc(nImages);
+
+    for (int i = 0; i < nImages; i++) {
+
+	// find the currently selected readout
+	pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest
+	psAssert (file, "missing file?");
+
+	pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+	psAssert (readout, "missing readout?");
+
+	pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+	psAssert (detections, "missing detections?");
+
+	detArrays->data[i] = psMemIncrRefCounter(detections);
+	readouts->data[i] = psMemIncrRefCounter(readout);
+    }
+
+    // vector to track if source for an image is found
+    psVector *found = psVectorAlloc(nImages, PS_TYPE_U8);
+
+    for (int i = 0; i < objects->n; i++) { 
+        pmPhotObj *obj = objects->data[i]; 
+
+	// mark the images for which sources have been found
+	psVectorInit (found, 0);
+	for (int j = 0; j < obj->sources->n; j++) {
+
+	    pmSource *src = obj->sources->data[j]; 
+	    int index = src->imageID;
+	    psAssert (index >= 0, "invalid index");
+	    psAssert (index < found->n, "invalid index");
+
+	    found->data.U8[index] = 1;
+	}
+
+	// generate new sources for the image that are missing
+	for (int index = 0; index < found->n; index++) {
+	    if (found->data.U8[index]) continue;
+
+	    pmDetections *detections = detArrays->data[index];
+	    pmReadout *readout = readouts->data[index];
+	    int row0 = readout->image->row0;
+	    int col0 = readout->image->col0;
+
+	    // XXX the peak type is not really used in psphot
+	    // PM_PEAK_LONE is certainly not true, but irrelevant
+	    float peakFlux = readout->image->data.F32[(int)(obj->y-row0-0.5)][(int)(obj->x-col0-0.5)];
+	    pmPeak *peak = pmPeakAlloc(obj->x, obj->y, peakFlux, PM_PEAK_LONE);
+	    peak->flux = peakFlux;
+	    peak->SN = 1.0;
+	    peak->xf = obj->x;
+	    peak->yf = obj->y;
+	    peak->dx = NAN;
+	    peak->dy = NAN;
+	    
+	    // XXX assign to a footprint?
+
+	    // create a new source
+	    pmSource *source = pmSourceAlloc();
+	    source->imageID = index;
+
+	    // add the peak
+	    source->peak = psMemIncrRefCounter(peak);
+
+	    // allocate space for moments
+	    source->moments = pmMomentsAlloc();
+
+	    // allocate image, weight, mask arrays for each peak (square of radius OUTER)
+	    pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
+
+	    peak->assigned = true;
+	    pmPhotObjAddSource(obj, source);
+	    psArrayAdd (detections->newSources, 100, source);
+	    psFree (source);
+	    psFree (peak);
+	}
+    }
+
+    // how many sources do we have now?
+    int nSources = 0;
+    for (int i = 0; i < objects->n; i++) { 
+        pmPhotObj *obj = objects->data[i]; 
+	nSources += obj->sources->n;
+    }
+    psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages);
+
+
+    psFree (found);
+    psFree (detArrays);
+    psFree (readouts);
+    return true;
+}
Index: branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceStats.c
===================================================================
--- branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceStats.c	(revision 27649)
+++ branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceStats.c	(revision 27650)
@@ -2,5 +2,7 @@
 
 // convert detections to sources and measure their basic properties (moments, local sky, sky
-// variance) Note: this function only generates sources for the new peaks (peak->assigned)
+// variance) Note: this function only generates sources for the new peaks (peak->assigned).
+// The new sources are added to any existing sources on detections->newSources.  The sources 
+// on detections->allSources are ignored.
 bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow)
 {
@@ -40,5 +42,4 @@
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     psAssert (detections, "missing detections?");
-    psAssert (!detections->newSources, "new sources already defined?");
 
     // XXX TEST:
@@ -84,10 +85,11 @@
 
     // generate the array of sources, define the associated pixel
-    sources = psArrayAllocEmpty (peaks->n);
+    if (!detections->newSources) {
+	detections->newSources = psArrayAllocEmpty (peaks->n);
+    }
+    sources = detections->newSources;
 
     // if there are no peaks, we save the empty source array and return
     if (!peaks->n) {
-        // save the new sources on the detection structure:
-        detections->newSources = sources;
         return true;
     }
@@ -100,4 +102,5 @@
         // create a new source
         pmSource *source = pmSourceAlloc();
+	source->imageID = index;
 
         // add the peak
@@ -119,5 +122,4 @@
         psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
         psphotVisualShowMoments (sources);
-        detections->newSources = sources;
         return true;
     }
@@ -126,5 +128,5 @@
         if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
             psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
-            psFree(sources);
+            psFree(detections->newSources);
             return false;
         }
@@ -178,5 +180,5 @@
                 psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
                 psFree (job);
-                psFree(sources);
+                psFree(detections->newSources);
                 return false;
             }
@@ -187,5 +189,5 @@
         if (!psThreadPoolWait (false)) {
             psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
-            psFree(sources);
+	    psFree(detections->newSources);
             return false;
         }
@@ -215,8 +217,4 @@
     psphotVisualShowMoments (sources);
 
-    // save the new sources on the detection structure:
-    detections->newSources = sources;
-
-
     if (detections->allSources) {
         psphotMaskCosmicRayFootprintCheck(detections->allSources);
@@ -377,4 +375,7 @@
         pmSource *source = sources->data[i];
 
+	if (source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue;
+	source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
+
         // skip faint sources for moments measurement
         if (source->peak->SN < MIN_SN) {