Index: trunk/psModules/src/objects/pmPCMdata.c
===================================================================
--- trunk/psModules/src/objects/pmPCMdata.c	(revision 31451)
+++ trunk/psModules/src/objects/pmPCMdata.c	(revision 32347)
@@ -41,4 +41,7 @@
 #include "pmPCMdata.h"
 
+# define USE_DELTA_PSF 0
+# define USE_1D_GAUSS 1
+
 static void pmPCMdataFree (pmPCMdata *pcm) {
 
@@ -88,4 +91,10 @@
     pcm->constraint = NULL;
     pcm->nDOF = 0;
+
+    // full convolution with the PSF is expensive.  if we have to save time, we can do a 1D
+    // convolution with a Gaussian approximation to the kernel
+    pcm->use1Dgauss = false;
+    pcm->nsigma = 3.0; 
+    pcm->sigma = 1.0; // this should be set to something sensible when the psf is known
 
     return pcm;
@@ -257,4 +266,24 @@
     pcm->nDOF = nPix - nParams;
 
+# if (USE_1D_GAUSS)
+    pmModel *modelPSF = source->modelPSF;
+    psAssert (modelPSF, "psf model must be defined");
+    
+    psEllipseShape shape;
+    psEllipseAxes axes;
+
+    shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
+    shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
+    shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
+    axes = psEllipseShapeToAxes (shape, 20.0);
+    
+    float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[PM_PAR_I0]);
+    float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+
+    pcm->use1Dgauss = true;
+    pcm->sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
+    pcm->nsigma = 2.0;
+# endif
+
     return pcm;
 }
@@ -368,2 +397,62 @@
     return true;
 }
+
+// construct a realization of the source model
+bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize) {
+
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
+    // select appropriate model
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+
+    // if we already have a cached image, re-use that memory
+    source->modelFlux = psImageCopy (source->modelFlux, source->pixels, PS_TYPE_F32);
+    psImageInit (source->modelFlux, 0.0);
+
+    // modelFlux always has unity normalization (I0 = 1.0)
+    pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+
+    // convolve the model image with the PSF
+    if (USE_1D_GAUSS) {
+	// do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
+	// * the model flux is not masked
+	// * threading takes place above this level
+	
+	// define the Gauss parameters from the psf
+	pmModel *modelPSF = source->modelPSF;
+	psAssert (modelPSF, "psf model must be defined");
+    
+	psEllipseShape shape;
+	psEllipseAxes axes;
+
+	shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
+	shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
+	shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
+	axes = psEllipseShapeToAxes (shape, 20.0);
+    
+	float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[PM_PAR_I0]);
+	float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+
+	float sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
+	float nsigma = 2.0;
+
+	psImageSmooth (source->modelFlux, sigma, nsigma);
+    } else {
+	// make sure we save a cached copy of the psf flux
+	pmSourceCachePSF (source, maskVal);
+
+	// convert the cached cached psf model for this source to a psKernel
+	psKernel *psf = pmPCMkernelFromPSF (source, psfSize);
+	if (!psf) {
+	    // NOTE: this only happens if the source is too close to an edge
+	    model->flags |= PM_MODEL_STATUS_BADARGS;
+	    return NULL;
+	}
+
+	// XXX not sure if I can place the output on top of the input
+	psImageConvolveFFT (source->modelFlux, source->modelFlux, NULL, 0, psf);
+    }
+    return true;
+}
+