Index: trunk/ppSim/src/Makefile.am
===================================================================
--- trunk/ppSim/src/Makefile.am	(revision 28006)
+++ trunk/ppSim/src/Makefile.am	(revision 28125)
@@ -43,4 +43,5 @@
 	ppSimBounds.c             \
 	ppSimStars.c              \
+	ppSimSmoothReadout.c      \
 	ppSimSetPSF.c             \
 	ppSimUtils.c              \
Index: trunk/ppSim/src/ppSim.h
===================================================================
--- trunk/ppSim/src/ppSim.h	(revision 28006)
+++ trunk/ppSim/src/ppSim.h	(revision 28125)
@@ -106,4 +106,6 @@
 bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config);
 
+bool ppSimSmoothReadout(pmReadout *input, psMetadata *recipe);
+
 bool ppSimInitHeader(pmConfig *config,
                      pmFPA *fpa,
Index: trunk/ppSim/src/ppSimInsertStars.c
===================================================================
--- trunk/ppSim/src/ppSimInsertStars.c	(revision 28006)
+++ trunk/ppSim/src/ppSimInsertStars.c	(revision 28125)
@@ -23,4 +23,6 @@
     // XXX this is missing the gain.
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
+
+    bool psfConvolve = psMetadataLookupBool(NULL, recipe, "PSF.CONVOLVE"); // smooth the image with the PSF?
 
     float expTime   = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time
@@ -103,43 +105,50 @@
         float flux = star->flux * expCorr->data.F32[(int)yCell][(int)xCell];
 
-        // instantiate a model for the PSF at this location, set desired flux
-        pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
-        pmModelSetFlux (model, flux);
+	// if psfConvolve is TRUE, we will (elsewhere) convolve the image we a PSF
+	// in this case, simply place delta functions in the image
+	if (psfConvolve) {
+	    readout->image->data.F32[(int)(yCell)][(int)(xCell)] = flux;
+	    continue;
+	}
 
-        float radius = model->modelRadius (model->params, nSigmaLim * roughNoise);
-        radius = PS_MAX (radius, 1.0);
+	// instantiate a model for the PSF at this location, set desired flux
+	pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
+	pmModelSetFlux (model, flux);
 
-        // construct a source, with model flux pixels set, based on the model
-        pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
+	float radius = model->modelRadius (model->params, nSigmaLim * roughNoise);
+	radius = PS_MAX (radius, 1.0);
 
-        // XXX set the mag & err values (should this be done in pmSourceFromModel?)
-        // XXX i should be applying the gain and the correct effective area
-        psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
-        psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
+	// construct a source, with model flux pixels set, based on the model
+	pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
 
-        // this value is the pure (input) flux, and is saved in the output source cmf files
-        source->psfMag = -2.5*log10(star->flux);
-        source->errMag = sqrt(Area*PS_SQR(roughNoise) + flux) / flux;
+	// XXX set the mag & err values (should this be done in pmSourceFromModel?)
+	// XXX i should be applying the gain and the correct effective area
+	psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
+	psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
 
-        fprintf (outfile, "%8.3f %8.3f  %10.2f  %7.3f %5.3f\n", star->x, star->y, star->flux, source->psfMag, source->errMag);
+	// this value is the pure (input) flux, and is saved in the output source cmf files
+	source->psfMag = -2.5*log10(star->flux);
+	source->errMag = sqrt(Area*PS_SQR(roughNoise) + flux) / flux;
 
-        // insert the source flux in the image
-        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
+	fprintf (outfile, "%8.3f %8.3f  %10.2f  %7.3f %5.3f\n", star->x, star->y, star->flux, source->psfMag, source->errMag);
 
-        // insert the source flux in the noise image
-        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, dX, dY);
+	// insert the source flux in the image
+	pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
 
-        // Blow away the image parts of the source, which makes the memory explode
-        RESET(source->pixels);
-        RESET(source->variance);
-        RESET(source->maskObj);
-        RESET(source->maskView);
-        RESET(source->modelFlux);
-        RESET(source->psfImage);
-        RESET(source->blends);
+	// insert the source flux in the noise image
+	pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, dX, dY);
 
-        // add the sources to the source array
-        psArrayAdd (sources, 100,source);
-        psFree(source);                 // Drop reference
+	// Blow away the image parts of the source, which makes the memory explode
+	RESET(source->pixels);
+	RESET(source->variance);
+	RESET(source->maskObj);
+	RESET(source->maskView);
+	RESET(source->modelFlux);
+	RESET(source->psfImage);
+	RESET(source->blends);
+
+	// add the sources to the source array
+	psArrayAdd (sources, 100,source);
+	psFree(source);                 // Drop reference
     }
     fclose (outfile);
Index: trunk/ppSim/src/ppSimLoop.c
===================================================================
--- trunk/ppSim/src/ppSimLoop.c	(revision 28006)
+++ trunk/ppSim/src/ppSimLoop.c	(revision 28125)
@@ -30,4 +30,6 @@
     ppSimType type = ppSimTypeFromString (typeStr); // Type of image to simulate
     int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y
+
+    bool psfConvolve = psMetadataLookupBool(NULL, recipe, "PSF.CONVOLVE"); // smooth the image with the PSF?
 
     ppSimUpdateConceptsFPA (fpa, config);
@@ -155,4 +157,11 @@
                 psFree(expCorr);
 
+		// we have two options for generating images which have a PSF:
+		// 1) lay down stars with the PSF model applied : in this case, galaxies do NOT have the PSF
+		// 2) lay down delta functions for stars and smooth the image with a PSF : in this case, the stars land at integer pixel locations
+		if (psfConvolve) {
+		    if (!ppSimSmoothReadout(readout, recipe)) ESCAPE (PS_ERR_UNKNOWN, "problem smoothing image");
+		}
+
             done:
                 if (!ppSimAddNoise(readout->image, readout->variance, cell, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding noise");
Index: trunk/ppSim/src/ppSimSmoothReadout.c
===================================================================
--- trunk/ppSim/src/ppSimSmoothReadout.c	(revision 28125)
+++ trunk/ppSim/src/ppSimSmoothReadout.c	(revision 28125)
@@ -0,0 +1,29 @@
+#include "ppSim.h"
+
+bool ppSimSmoothReadout(pmReadout *input, psMetadata *recipe)
+{
+    bool status;
+
+    // XXX use these defaults?
+    // float minGauss = 0.1;
+    float nSigma = 3.0;
+    float sigma = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels)
+
+    // bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading in psImageConvolve
+
+    psTimerStart ("ppSmooth");
+
+    // smooth the image in place, applying the mask as we go
+    psImageSmooth(input->image, sigma, nSigma);
+    psLogMsg("ppSmooth", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark("ppSmooth"));
+
+    // psImageConvolveSetThreads(oldThreads);
+    return true;
+}
+
+// The variance is smoothed by the PSF^2, renomalized to maintain the input level of the
+// variance.  We achieve this by smoothing with a Gaussian with sigma = SIGMA_SMTH/sqrt(2) with
+// unity normalization.  Note that this process yields a smoothed image with correlated errors.
+// The pixel-to-pixel variations in smooth_im will be decreased by a factor of
+// 4*pi*SIGMA_SMTH^2, but for measurements based on apertures comparable to or larger than the
+// smoothing kernel, the effective per-pixel variance is maintained.