Index: trunk/psModules/src/objects/pmPCMdata.c
===================================================================
--- trunk/psModules/src/objects/pmPCMdata.c	(revision 35768)
+++ trunk/psModules/src/objects/pmPCMdata.c	(revision 36085)
@@ -173,4 +173,78 @@
 }
 
+int pmPCMsetParams (psMinConstraint *constraint, pmSourceFitMode mode) {
+
+    // set parameter mask based on fitting mode
+    int nParams = 0;
+
+    switch (mode) {
+      case PM_SOURCE_FIT_NORM:
+        // fits only source normalization (Io)
+        nParams = 1;
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+        break;
+
+      case PM_SOURCE_FIT_PSF:
+        // fits only x,y,Io
+        nParams = 3;
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
+        break;
+
+      case PM_SOURCE_FIT_EXT:
+        // fits all params except sky
+        nParams = params->n - 1;
+        psVectorInit (constraint->paramMask, 0);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
+        break;
+
+      case PM_SOURCE_FIT_EXT_AND_SKY:
+        // fits all params including sky
+        nParams = params->n;
+        psVectorInit (constraint->paramMask, 0);
+        break;
+
+      case PM_SOURCE_FIT_SHAPE:
+	// fits shape (Sxx, Sxy, Syy) and Io
+	nParams = 5;
+	psVectorInit (constraint->paramMask, 1);
+	constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 0;
+	constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+	constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SXX] = 0;
+	constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SXY] = 0;
+	constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SYY] = 0;
+	break;
+
+      case PM_SOURCE_FIT_INDEX:
+        // fits only Io, index (PAR7) -- only Io for models with < 8 params
+	psVectorInit (constraint->paramMask, 1);
+	constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+        if (params->n == 7) {
+	    nParams = 1;
+	} else {
+	    nParams = 2;
+	    constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
+	}
+	break;
+
+      case PM_SOURCE_FIT_NO_INDEX:
+        // fits all but index (PAR7) including sky
+	psVectorInit (constraint->paramMask, 0);
+        if (params->n == 7) {
+	    nParams = params->n;
+	} else {
+	    nParams = params->n - 1;
+	    constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
+	}
+	break;
+      default:
+	psAbort("invalid fitting mode");
+    }
+    return nParams;
+}
+
 pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) {
 
@@ -219,52 +293,5 @@
     constraint->checkLimits = model->modelLimits;
 
-    // set parameter mask based on fitting mode
-    int nParams = 0;
-    switch (fitOptions->mode) {
-      case PM_SOURCE_FIT_NORM:
-        // NORM-only model fits only source normalization (Io)
-        nParams = 1;
-        psVectorInit (constraint->paramMask, 1);
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-        break;
-      case PM_SOURCE_FIT_PSF:
-        // PSF model only fits x,y,Io
-        nParams = 3;
-        psVectorInit (constraint->paramMask, 1);
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
-        break;
-      case PM_SOURCE_FIT_EXT:
-        // EXT model fits all params (except sky)
-        nParams = params->n - 1;
-        psVectorInit (constraint->paramMask, 0);
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
-        break;
-      case PM_SOURCE_FIT_INDEX:
-        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
-	psVectorInit (constraint->paramMask, 1);
-	constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-        if (params->n == 7) {
-	    nParams = 1;
-	} else {
-	    nParams = 2;
-	    constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
-	}
-	break;
-      case PM_SOURCE_FIT_NO_INDEX:
-        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
-	psVectorInit (constraint->paramMask, 0);
-	constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
-        if (params->n == 7) {
-	    nParams = params->n - 1;
-	} else {
-	    nParams = params->n - 2;
-	    constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
-	}
-	break;
-      default:
-	psAbort("invalid fitting mode");
-    }
+    int nParams = pmPCMsetParams (constraint, fitOptions->mode);
 
     if (nPix <  nParams + 1) {
@@ -341,53 +368,5 @@
     }
 
-    // if we changed the fit mode, we need to update nDOF
-    int nParams = 0;
-    // set parameter mask based on fitting mode
-    switch (fitOptions->mode) {
-      case PM_SOURCE_FIT_NORM:
-	// NORM-only model fits only source normalization (Io)
-	nParams = 1;
-	psVectorInit (pcm->constraint->paramMask, 1);
-	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-	break;
-      case PM_SOURCE_FIT_PSF:
-	// PSF model only fits x,y,Io
-	nParams = 3;
-	psVectorInit (pcm->constraint->paramMask, 1);
-	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
-	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
-	break;
-      case PM_SOURCE_FIT_EXT:
-	// EXT model fits all params (except sky)
-	nParams = model->params->n - 1;
-	psVectorInit (pcm->constraint->paramMask, 0);
-	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
-	break;
-      case PM_SOURCE_FIT_INDEX:
-	// PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
-	psVectorInit (pcm->constraint->paramMask, 1);
-	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-	if (model->params->n == 7) {
-	    nParams = 1;
-	} else {
-	    nParams = 2;
-	    pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
-	}
-	break;
-      case PM_SOURCE_FIT_NO_INDEX:
-	// PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
-	psVectorInit (pcm->constraint->paramMask, 0);
-	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
-	if (model->params->n == 7) {
-	    nParams = model->params->n - 1;
-	} else {
-	    nParams = model->params->n - 2;
-	    pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
-	}
-	break;
-      default:
-	psAbort("invalid fitting mode");
-    }
+    int nParams = pmPCMsetParams (pcm->constraint, fitOptions->mode);
 
     if (pcm->nPix <  nParams + 1) {
@@ -478,2 +457,54 @@
 }
 
+// construct a realization of the source model
+bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize) {
+
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
+    // 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");
+    
+	psEllipseAxes axes;
+	bool useReff = pmModelUseReff (modelPSF->type);
+	psF32 *PAR = modelPSF->params->data.F32;
+	pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
+    
+	float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[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;
+}
+