Index: /trunk/psphot/src/psphotAddNoise.c
===================================================================
--- /trunk/psphot/src/psphotAddNoise.c	(revision 12626)
+++ /trunk/psphot/src/psphotAddNoise.c	(revision 12626)
@@ -0,0 +1,54 @@
+# include "psphot.h"
+
+bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add) {
+
+    bool status = false;
+    psTimerStart ("psphot");
+
+    PS_ASSERT (readout, false);
+    PS_ASSERT (readout->parent, false);
+    PS_ASSERT (readout->parent->concepts, false);
+
+    // increase variance by factor*(object noise):
+    // weight = flux/gain + rn^2/g^2
+    // we are adding factor*flux/gain
+
+    float FACTOR = psMetadataLookupF32 (&status, recipe, "NOISE.FACTOR");
+    PS_ASSERT (status, false);
+
+    float GAIN = psMetadataLookupF32(&status, readout->parent->concepts, "CELL.GAIN"); // Cell gain
+    PS_ASSERT (status, false);
+
+    FACTOR /= GAIN;
+
+    // loop over all source
+    for (int i = 0; i < sources->n; i++) {
+	source = sources->data[i];
+
+	// skip sources which were not subtracted
+	if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
+
+	// select appropriate model
+	pmModel *model = pmSourceGetModel (NULL, source);
+	if (model == NULL) continue;  // model must be defined
+	
+	if (add) {
+	    psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+	} else {
+	    psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+	}
+
+	model->params->data.F32[PM_PAR_I0] *= FACTOR;
+	if (add) {
+	    pmModelAdd (source->weight, source->mask, model, false, false);
+	} else {
+	    pmModelSub (source->weight, source->mask, model, false, false);
+	}
+    }
+    if (add) {
+	psLogMsg ("psphot.noise", PS_LOG_INFO, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));
+    } else {
+	psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));
+    }
+    return true;
+}
Index: /trunk/psphot/src/psphotReadout.c
===================================================================
--- /trunk/psphot/src/psphotReadout.c	(revision 12625)
+++ /trunk/psphot/src/psphotReadout.c	(revision 12626)
@@ -126,6 +126,12 @@
     psphotImageMedian (config, view);
 
+    // add noise for subtracted objects
+    psphotAddNoise (readout, sources, recipe, true);
+
     // find the peaks in the image
     psArray *newPeaks = psphotFindPeaks (readout, recipe, 2);
+
+    // remove noise for subtracted objects
+    psphotAddNoise (readout, sources, recipe, false);
 
     // define new sources based on the new peaks
Index: /trunk/psphot/src/psphotSourcePlots.c
===================================================================
--- /trunk/psphot/src/psphotSourcePlots.c	(revision 12626)
+++ /trunk/psphot/src/psphotSourcePlots.c	(revision 12626)
@@ -0,0 +1,74 @@
+# include "psphot.h"
+
+bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe) {
+
+    bool status = false;
+    psTimerStart ("psphot");
+
+    // the source images are written to an image 10x the size of a PSF object
+    float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
+    PS_ASSERT (status);
+
+    int NX = 10*(2*OUTER+1);
+    int NY = 0;
+
+    // examine PSF sources in S/N order (brightest first)
+    sources = psArraySort (sources, psphotSortBySN);
+
+    bool newRow = true;
+
+    // first, examine the PSF and SAT stars:
+    // - generate radial plots (PS plots)
+    // - determine bounding boxes for summary image
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+	bool keep = false;
+        keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
+        keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
+	if (!keep) continue;
+
+	if (newRow) {
+	    dX = 0;
+	    dY = 0;
+	    newRow = false;
+	}
+
+	if (dX + source->pixels->numCols > NX) {
+	    // go to the next row
+	    // but, if we have only this one entry, we will truncate it.
+	}
+	dX += source->pixels->numCols;
+	dY = PS_MAX (dY, source->pixels->numRows);
+
+	// for a new row:
+	NY += dY;
+
+	RadialPlot (source);
+    }
+
+    // allocate output image
+    psImage *output = psImageAlloc (NX, NY, PS_TYPE_F32);
+
+    // first, examine the PSF and SAT stars:
+    // - generate radial plots (PS plots)
+    // - determine bounding boxes for summary image
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+	bool keep = false;
+        keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
+        keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
+	if (!keep) continue;
+
+	// add to image
+    }
+
+    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d failed: %f sec\n", sources->n, Nmoments, Nfail, psTimerMark ("psphot"));
+
+    return (sources);
+}
+
+// XXX EAM : filter out bad peaks (eg, no valid pixels available for sky)