Index: trunk/psModules/src/objects/pmPCMdata.c
===================================================================
--- trunk/psModules/src/objects/pmPCMdata.c	(revision 36089)
+++ trunk/psModules/src/objects/pmPCMdata.c	(revision 36375)
@@ -43,5 +43,4 @@
 
 # define USE_DELTA_PSF 0
-# define USE_1D_GAUSS 1
 
 static void pmPCMdataFree (pmPCMdata *pcm) {
@@ -58,5 +57,7 @@
     psFree (pcm->psfFFT);
     psFree (pcm->constraint);
+
     psFree (pcm->smdata); // pre-allocated data for psImageSmooth_PreAlloc
+    psFree (pcm->smdata2d); // pre-allocated data for psImageSmooth_PreAlloc
     return;
 }
@@ -88,4 +89,7 @@
     }
 
+    pcm->smdata = NULL;
+    pcm->smdata2d = NULL;
+
     pcm->modelConv = NULL;
     pcm->psf = NULL;
@@ -94,8 +98,10 @@
     pcm->nDOF = 0;
 
+    pcm->poissonErrors = true;
+
     // 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->nsigma = NAN; // this is set to something defined by the user
     pcm->sigma = 1.0; // this should be set to something sensible when the psf is known
 
@@ -248,21 +254,115 @@
 }
 
+static int modelType_GAUSS = -1;
+static int modelType_PS1_V1 = -1;
+
+// generate a Gaussian smoothing kernel for supplied sigma.  sigma here does not need to match
+// that used to allocate the structure, but it is recommended
+bool psImageSmoothCacheKernel_PS1_V1 (psImageSmoothCacheData *smdata, float sigma, float kappa) {
+    // check for NULL structure elements?
+
+    int size = smdata->Nrange;
+
+    psFree (smdata->kernel);
+    smdata->kernel = psVectorAlloc(2 * smdata->Nrange + 1, PS_TYPE_F32);
+
+    double sum = 0.0;			// Sum of Gaussian, for normalization
+    double factor = 1.0 / (sigma * M_SQRT2);	// Multiplier for i -> z
+
+    // PS1_V1 is a power-law with fitted linear term:
+    // 1 / (1 + kappa z + z^1.666)  where z = (r/sigma)^2
+
+    // generate the kernel (not normalized)
+    for (int i = -size, j = 0; i <= size; i++, j++) {
+	float z = PS_SQR(i * factor);
+        sum += smdata->kernel->data.F32[j] = 1.0 / (1 + kappa * z + pow(z,1.666));
+    }
+
+    // renormalize kernel to integral of 1.0
+    for (int i = 0; i < 2 * size + 1; i++) {
+        smdata->kernel->data.F32[i] /= sum;
+    }
+
+    return true;
+}
+
+psImageSmoothCacheData *psImageSmoothCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *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);
+    
+    *sigma = NAN;
+    *kappa = NAN;
+
+    // XXX need to do this more carefully
+    if (modelPSF->type == modelType_GAUSS) {
+	float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
+	float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+	*sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+	*sigma = 0.5 * (axes.major + axes.minor);
+	*kappa = PAR[PM_PAR_7];
+    }
+    psAssert (isfinite(*sigma), "invalid model type");
+
+    // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
+    psImageSmoothCacheData *smdata = psImageSmoothCacheAlloc (flux, *sigma, nsigma);
+
+    if (modelPSF->type == modelType_GAUSS) {
+	psImageSmoothCacheKernel_Gauss (smdata, *sigma);
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+	psImageSmoothCacheKernel_PS1_V1 (smdata, *sigma, *kappa);
+    }
+
+    return smdata;
+}
+
+psImageSmooth2dCacheData *psImageSmooth2dCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *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);
+    
+    *sigma = NAN;
+    *kappa = NAN;
+
+    // XXX need to do this more carefully
+    if (modelPSF->type == modelType_GAUSS) {
+	float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
+	float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+	*sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+	*sigma = 0.5 * (axes.major + axes.minor);
+	*kappa = PAR[PM_PAR_7];
+    }
+    psAssert (isfinite(*sigma), "invalid model type");
+
+    // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
+    psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (nsigma);
+
+    if (modelPSF->type == modelType_GAUSS) {
+	psImageSmooth2dCacheKernel_Gauss (smdata, *sigma);
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+	psImageSmooth2dCacheKernel_PS1_V1 (smdata, *sigma, *kappa);
+    }
+
+    return smdata;
+}
+
 pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) {
 
-    // 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;
-    }
-
-# if (USE_DELTA_PSF)
-    psImageInit (psf->image, 0.0);
-    psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0;
-# endif
+    modelType_GAUSS = pmModelClassGetType ("PS_MODEL_GAUSS");
+    modelType_PS1_V1 = pmModelClassGetType ("PS_MODEL_PS1_V1");
 
     // count the number of unmasked pixels:
@@ -298,5 +398,4 @@
     if (nPix <  nParams + 1) {
         psTrace ("psModules.objects", 4, "insufficient valid pixels\n");
-	psFree (psf);
 	psFree (constraint);
         model->flags |= PM_MODEL_STATUS_BADARGS;
@@ -306,8 +405,9 @@
     // generate PCM data storage structure
     pmPCMdata *pcm = pmPCMdataAlloc (params, constraint->paramMask, source);
-
-    pcm->psf = psf;
     pcm->modelConv = psMemIncrRefCounter(model);
     pcm->constraint = constraint;
+
+    pcm->poissonErrors = fitOptions->poissonErrors;
+    pcm->nsigma = fitOptions->nsigma;
 
     pcm->nPix = nPix;
@@ -316,21 +416,29 @@
 
 # if (USE_1D_GAUSS)
-    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);
 
     pcm->use1Dgauss = true;
-    pcm->sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
-    pcm->nsigma = 2.0;
-
-    pcm->smdata = psImageSmooth_PreAlloc_DataAlloc (source->pixels, pcm->sigma, pcm->nsigma);
+    if (USE_1D_CACHE) {
+	pcm->smdata = psImageSmoothCacheSetKernel (&pcm->sigma, &pcm->kappa, pcm->nsigma, source->pixels, source->modelPSF);
+    } else {
+	pcm->smdata2d = psImageSmooth2dCacheSetKernel (&pcm->sigma, &pcm->kappa, pcm->nsigma, source->pixels, source->modelPSF);
+    }
+
 # 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;
+    }
+
+# if (USE_DELTA_PSF)
+    psImageInit (psf->image, 0.0);
+    psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0;
+# endif
+    pcm->psf = psf;
     pcm->smdata = NULL;
 # endif
@@ -395,6 +503,29 @@
 	    pcm->dmodelsConvFlux->data[n] = psImageCopy (pcm->dmodelsConvFlux->data[n], source->pixels, PS_TYPE_F32);
 	}
-	psFree(pcm->smdata);
-	pcm->smdata = psImageSmooth_PreAlloc_DataAlloc (source->pixels, pcm->sigma, pcm->nsigma);
+
+	// If we have changed the window, we need to redefine the smoothing target vectors (but pcm->sigma,kappa,nsigma remain)
+	if (USE_1D_CACHE) {
+	    psFree(pcm->smdata);
+	    pcm->smdata = psImageSmoothCacheAlloc (source->pixels, pcm->sigma, pcm->nsigma);
+
+	    pmModel *modelPSF = source->modelPSF;
+	    if (modelPSF->type == modelType_GAUSS) {
+		psImageSmoothCacheKernel_Gauss (pcm->smdata, pcm->sigma);
+	    }
+	    if (modelPSF->type == modelType_PS1_V1) {
+		psImageSmoothCacheKernel_PS1_V1 (pcm->smdata, pcm->sigma, pcm->kappa);
+	    }
+	} else {
+	    psFree(pcm->smdata2d);
+	    pcm->smdata2d = psImageSmooth2dCacheAlloc (pcm->nsigma);
+
+	    pmModel *modelPSF = source->modelPSF;
+	    if (modelPSF->type == modelType_GAUSS) {
+		// psImageSmooth2dCacheKernel_Gauss (pcm->smdata2d, pcm->sigma);
+	    }
+	    if (modelPSF->type == modelType_PS1_V1) {
+		psImageSmooth2dCacheKernel_PS1_V1 (pcm->smdata2d, pcm->sigma, pcm->kappa);
+	    }
+	}
     }
 
@@ -403,5 +534,5 @@
 
 // construct a realization of the source model
-bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize) {
+bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize, float nsigma) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -420,24 +551,18 @@
     // 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);
+	float sigma = NAN;
+	float kappa = NAN;
+
+	if (USE_1D_CACHE) {
+	    psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, nsigma, source->modelFlux, source->modelPSF);
+	    psImageSmoothCache_F32 (source->modelFlux, smdata);
+	    psFree (smdata);
+	} else {
+	    psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, nsigma, source->modelFlux, source->modelPSF);
+	    psImageSmooth2dCache_F32 (source->modelFlux, smdata);
+	    psFree (smdata);
+	}
+	// old call: psImageSmooth (source->modelFlux, sigma, nsigma);
     } else {
 	// make sure we save a cached copy of the psf flux
@@ -459,5 +584,5 @@
 
 // construct a realization of the source model
-bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize) {
+bool pmPCMMakeModel (pmSource *source, pmModel *model, float Nsigma, psImageMaskType maskVal, int psfSize) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -468,28 +593,23 @@
 
     // modelFlux always has unity normalization (I0 = 1.0)
-    pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+    // pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+    pmModelAdd (source->modelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_SKY | 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);
+
+	float sigma = NAN;
+	float kappa = NAN;
+
+	if (USE_1D_CACHE) {
+	    psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, Nsigma, source->modelFlux, source->modelPSF);
+	    psImageSmoothCache_F32 (source->modelFlux, smdata);
+	    psFree (smdata);
+	} else {
+	    psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, Nsigma, source->modelFlux, source->modelPSF);
+	    psImageSmooth2dCache_F32 (source->modelFlux, smdata);
+	    psFree (smdata);
+	}
+	// old call: psImageSmooth (source->modelFlux, sigma, nsigma);
     } else {
 	// make sure we save a cached copy of the psf flux
@@ -509,3 +629,2 @@
     return true;
 }
-