Index: /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c
===================================================================
--- /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c	(revision 36020)
+++ /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c	(revision 36021)
@@ -63,7 +63,8 @@
 # define PM_MODEL_SET_LIMITS      pmModelSetLimits_DEV
 
-// f = exp(-z^0.125) 
+// f = exp(-kappa*r^(1/index)) 
+// f = exp(-kappa*z^(0.5/index)) 
+// index = 4, 0.5/index = 0.125
 # define ALPHA 0.125 
-// # define ALPHA 0.25 
 
 // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
@@ -86,5 +87,5 @@
 static float *paramsMinUse = paramsMinLax;
 static float *paramsMaxUse = paramsMaxLax;
-static float betaUse[] = { 2, 3e6, 5, 5, 3.0, 3.0, 0.5 };
+static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5 };
 
 static bool limitsApply = true;         // Apply limits?
@@ -113,19 +114,11 @@
     // r = sqrt(z)
     float q = kappa*pow(z,ALPHA);
-    psF32 f0 = exp(-q);
-
-    psF32 f1 = PAR[PM_PAR_I0]*f0;
-    psF32 f = PAR[PM_PAR_SKY] + f1;
+    float f0 = exp(-q);
 
     assert (isfinite(q));
     assert (isfinite(f0));
-    assert (isfinite(f1));
-    assert (isfinite(f));
-
-    // only worry about the central 4 pixels at most
-    // If I use DELTA = 0.2, I'm way off for the total flux
-    // If I use DELTA = 0.02, I'm totally good (but I am under on the total flux for R = 30 by 0.2 mags -- aperture failure)
-    // For DELTA = 0.02 & Rmin/Rmaj = 0.25, I'm over flux by 0.15 mags (due to the central pixel)
-    psF32 radius = hypot(X, Y);
+
+    // only worry about the central pixels at most
+    float radius = hypot(X, Y);
     if (radius <= 1.5) {
 	// Nsub ~ 10*index^2 + 1
@@ -135,7 +128,11 @@
 	Nsub = PS_MAX (Nsub, 11);
 	f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 4.0, Nsub);
-	f1 = PAR[PM_PAR_I0]*f0;
-	f = PAR[PM_PAR_SKY] + f1;
     }   
+
+    float f1 = PAR[PM_PAR_I0]*f0;
+    float f = PAR[PM_PAR_SKY] + f1;
+
+    assert (isfinite(f1));
+    assert (isfinite(f));
 
     if (deriv != NULL) {
@@ -304,5 +301,5 @@
     // f = Io exp(-kappa*z^n) -> z^n = ln(Io/f) / kappa
     psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
-    psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA);
+    psF64 radius = axes.major * pow(zn, 0.5 / ALPHA);
 
     psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
Index: /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_EXP.c
===================================================================
--- /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_EXP.c	(revision 36020)
+++ /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_EXP.c	(revision 36021)
@@ -79,5 +79,5 @@
 static float *paramsMinUse = paramsMinLax;
 static float *paramsMaxUse = paramsMaxLax;
-static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5};
+static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5};
 
 static bool limitsApply = true;         // Apply limits?
@@ -115,11 +115,5 @@
     psF32 f0 = exp(-q);
 
-    psF32 f1 = PAR[PM_PAR_I0]*f0;
-    psF32 f = PAR[PM_PAR_SKY] + f1;
-
     assert (isfinite(q));
-    assert (isfinite(f0));
-    assert (isfinite(f1));
-    assert (isfinite(f));
 
     // only worry about the central 4 pixels at most
@@ -127,7 +121,12 @@
     if (radius <= 1.5) {
 	f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 1.0, 51);
-	f1 = PAR[PM_PAR_I0]*f0;
-	f = PAR[PM_PAR_SKY] + f1;
-    }
+    }
+    assert (isfinite(f0));
+
+    psF32 f1 = PAR[PM_PAR_I0]*f0;
+    psF32 f = PAR[PM_PAR_SKY] + f1;
+
+    assert (isfinite(f1));
+    assert (isfinite(f));
 
     if (deriv != NULL) {
@@ -298,5 +297,5 @@
     // f = Io exp(-kappa*sqrt(z)) -> sqrt(z) = ln(Io/f) / kappa
     psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
-    psF64 radius = axes.major * sqrt (2.0) * zn;
+    psF64 radius = axes.major * zn;
 
     psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
Index: /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_SERSIC.c
===================================================================
--- /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_SERSIC.c	(revision 36020)
+++ /branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_SERSIC.c	(revision 36021)
@@ -20,5 +20,5 @@
    * note that a Sersic model is usually defined in terms of R_e, the half-light radius.  This
      construction does not include a factor of 2 in the X^2 term, etc, like for a Gaussian.
-     Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using setting:
+     Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using:
      shape.sx = SXX / sqrt(2), shape.sy = SYY / sqrt(2), shape.sxy = SXY, then calling
      psEllipseShapeToAxes, and multiplying the values of axes.major, axes.minor by sqrt(2)
@@ -55,4 +55,5 @@
 #include "pmPSFtry.h"
 #include "pmDetections.h"
+#include "pmModel_CentralPixel.h"
 
 #include "pmModel_SERSIC.h"
@@ -74,6 +75,6 @@
 
 // Lax parameter limits
-static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.05 };
-static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 4.0 };
+static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.1 };
+static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 1.0 };
 
 // Moderate parameter limits
@@ -88,9 +89,9 @@
 static float *paramsMinUse = paramsMinLax;
 static float *paramsMaxUse = paramsMaxLax;
-static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
+static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5, 1.0};
 
 static bool limitsApply = true;         // Apply limits?
 
-# include "pmModel_SERSIC.CP.h"
+// # include "pmModel_SERSIC.CP.h"
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -111,91 +112,33 @@
     psAssert (z >= 0, "do not allow negative z values in model");
 
-    float index = 0.5 / PAR[PM_PAR_7];
-    float par7 = PAR[PM_PAR_7];
-    float bn = 1.9992*index - 0.3271;
-    float Io = exp(bn);
-
-    psF32 f2 = bn*pow(z,par7);
-    psF32 f1 = Io*exp(-f2);
+    float Sindex = 0.5 / PAR[PM_PAR_7];
+    float kappa = pmSersicKappa (Sindex);
+
+    float q = kappa*pow(z,PAR[PM_PAR_7]);
+    psF32 f0 = exp(-q);
+
+    assert (isfinite(q));
 
     psF32 radius = hypot(X, Y);
-    if (radius < 1.0) {
-
-	// ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
-
-	// first, use Rmajor and index to find the central pixel flux (fraction of total flux)
-	psEllipseAxes axes;
-	pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
-
-	// get the central pixel flux from the lookup table
-	float xPix = (axes.major - centralPixelXo) / centralPixeldX;
-	xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
-	float yPix = (index - centralPixelYo) / centralPixeldY;
-	yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
-
-	// the integral of a Sersic has an analytical form as follows:
-	float logGamma = lgamma(2.0*index);
-	float bnFactor = pow(bn, 2.0*index);
-	float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
-
-	// XXX interpolate to get the value
-	// XXX for the moment, just integerize
-	// XXX I need to multiply by the integrated flux to get the flux in the central pixel
-	float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
-	
-	float px1 = 1.0 / PAR[PM_PAR_SXX];
-	float py1 = 1.0 / PAR[PM_PAR_SYY];
-	float z10 = PS_SQR(px1);
-	float z01 = PS_SQR(py1);
-
-	// which pixels do we need for this interpolation?
-	// (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
-	if ((X >= 0) && (Y >= 0)) {
-	    float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
-	    float V00 = Vcenter;
-	    float V10 = Io*exp(-bn*pow(z10,par7));
-	    float V01 = Io*exp(-bn*pow(z01,par7));
-	    float V11 = Io*exp(-bn*pow(z11,par7));
-	    f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
-	}
-	if ((X < 0) && (Y >= 0)) {
-	    float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
-	    float V00 = Io*exp(-bn*pow(z10,par7));
-	    float V10 = Vcenter;
-	    float V01 = Io*exp(-bn*pow(z11,par7));
-	    float V11 = Io*exp(-bn*pow(z01,par7));
-	    f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
-	}
-	if ((X >= 0) && (Y < 0)) {
-	    float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
-	    float V00 = Io*exp(-bn*pow(z01,par7));
-	    float V10 = Io*exp(-bn*pow(z11,par7));
-	    float V01 = Vcenter;
-	    float V11 = Io*exp(-bn*pow(z10,par7));
-	    f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
-	}
-	if ((X < 0) && (Y < 0)) {
-	    float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
-	    float V00 = Io*exp(-bn*pow(z11,par7));
-	    float V10 = Io*exp(-bn*pow(z10,par7));
-	    float V01 = Io*exp(-bn*pow(z01,par7));
-	    float V11 = Vcenter;
-	    f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
-	}
-    }   
-
-    psF32 z0 = PAR[PM_PAR_I0]*f1;
-    psF32 f0 = PAR[PM_PAR_SKY] + z0;
-
-    if (!isfinite(z0)) {
-        fprintf(stderr, "z0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, f1);
+    if (radius <= 1.5) {
+	// Nsub ~ 10*index^2 + 1
+	psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_SERSIC"));
+	int Nsub = 2 * ((int)(6.0*Sindex / axes.minor)) + 1;
+	Nsub = PS_MIN (Nsub, 121);
+	Nsub = PS_MAX (Nsub, 11);
+	f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], Sindex, Nsub);
+    }
+    if (!isfinite(f0)) {
+        fprintf(stderr, "f0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, q);
         fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4],
             PAR[5], PAR[6], PAR[7]);
     }
-
-    assert (isfinite(f2));
+    assert (isfinite(f0));
+
+    psF32 f1 = PAR[PM_PAR_I0]*f0;
+    psF32 f = PAR[PM_PAR_SKY] + f1;
+
     assert (isfinite(f1));
-    assert (isfinite(z0));
-    assert (isfinite(f0));
+    assert (isfinite(f));
 
     if (deriv != NULL) {
@@ -203,22 +146,27 @@
 
         dPAR[PM_PAR_SKY]  = +1.0;
-        dPAR[PM_PAR_I0]   = +f1;
-
-        // gradient is infinite for z = 0; saturate at z = 0.01
-        psF32 z1 = (z < 0.01) ? z0*bn*par7*pow(0.01,par7 - 1.0) : z0*bn*par7*pow(z,par7 - 1.0);
-
-        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,par7)*log(0.01) : -z0*f2*log(z);
-	dPAR[PM_PAR_7]   *= 3.0;
-
-        assert (isfinite(z1));
+        dPAR[PM_PAR_I0]   = +f0;
+
+	if (z > 0.01) {
+	  float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]-1.0);
+	  dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]);
+	  dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]);
+	  dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
+	  dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
+	  dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
+	  dPAR[PM_PAR_7]    = -1.0*f1*q*log(z);
+	} else {
+	  // gradient -> 0 for z -> 0, but has undef form
+	  float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]);
+	  dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]);
+	  dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]);
+	  dPAR[PM_PAR_SXX]  = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX];
+	  dPAR[PM_PAR_SYY]  = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY];
+	  dPAR[PM_PAR_SXY]  = -1.0*z1;
+	  dPAR[PM_PAR_7]    = -1.0*f1*q*log(z + 0.0001);
+	}
         assert (isfinite(dPAR[PM_PAR_7]));
-
-        dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
-        dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
-        dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag?
-        dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
-        dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
-    }
-    return (f0);
+    }
+    return (f);
 }
 
@@ -370,17 +318,11 @@
     psEllipseAxes axes;
     pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
-    float AspectRatio = axes.minor / axes.major;
-
-    float index = 0.5 / PAR[PM_PAR_7];
-    float bn = 1.9992*index - 0.3271;
-
-    // the integral of a Sersic has an analytical form as follows:
-    float logGamma = lgamma(2.0*index);
-    float bnFactor = pow(bn, 2.0*index);
-    float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
-    
-    psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
-
-    return(Flux);
+
+    float Sindex = 0.5 / PAR[PM_PAR_7];
+    float norm = pmSersicNorm (Sindex);
+
+    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
+
+    return(flux);
 }
 
@@ -401,7 +343,10 @@
     pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
+    float Sindex = 0.5 / PAR[PM_PAR_7];
+    float kappa = pmSersicKappa (Sindex);
+
     // f = Io exp(-z^n) -> z^n = ln(Io/f)
-    psF64 zn = log(PAR[PM_PAR_I0] / flux);
-    psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]);
+    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
+    psF64 radius = axes.major * pow(zn, Sindex);
 
     psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f", 
Index: /branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.c
===================================================================
--- /branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.c	(revision 36020)
+++ /branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.c	(revision 36021)
@@ -723,4 +723,51 @@
 }
 
+float pmSersicKappa (float Sindex) {
+    // this function is empirically derived from a fit to data for Sindex = 0.5 - 5.5
+    // constrain Sindex or kappa?
+    float kappa = -0.275552 + 1.972625*Sindex + 0.003487 * PS_SQR(Sindex);
+    return kappa;
+}
+
+float pmSersicNorm (float Sindex) {
+
+    float C0 = NAN;
+    float C1 = NAN;
+    float C2 = NAN;
+
+    // y = 0.201545 x^0 -0.950965 x^1 -0.315248 x^2 
+    // y = 0.402084 x^0 -1.357775 x^1 -0.105102 x^2 
+    // y = 0.619093 x^0 -1.591674 x^1 -0.041576 x^2 
+    // y = 0.770263 x^0 -1.696421 x^1 -0.023363 x^2 
+    // y = 0.885891 x^0 -1.755684 x^1 -0.015753 x^2 
+
+    if ((Sindex >= 0.0) && (Sindex < 1.0)) { 
+	C0 = 0.201545; C1 = -0.950965; C2 = -0.315248;
+	// y = 0.201545 x^0 -0.950965 x^1 -0.315248 x^2 
+    }
+    if ((Sindex >= 1.0) && (Sindex < 2.0)) { 
+	C0 = 0.402084; C1 = -1.357775; C2 = -0.105102;
+	// y = 0.402084 x^0 -1.357775 x^1 -0.105102 x^2 
+    }
+    if ((Sindex >= 2.0) && (Sindex < 3.0)) { 
+	C0 = 0.619093; C1 = -1.591674; C2 = -0.041576;
+	// y = 0.619093 x^0 -1.591674 x^1 -0.041576 x^2 
+    }
+    if ((Sindex >= 3.0) && (Sindex < 4.0)) { 
+	C0 = 0.770263; C1 = -1.696421; C2 = -0.023363;
+	// y = 0.770263 x^0 -1.696421 x^1 -0.023363 x^2 
+    }
+    if ((Sindex >= 4.0) && (Sindex < 5.5)) { 
+	C0 = 0.885891; C1 = -1.755684; C2 = -0.015753; 
+	// y = 0.885891 x^0 -1.755684 x^1 -0.015753 x^2 
+    }
+
+    if (isnan(C0)) return NAN;
+
+    float lnorm = C0 + C1*Sindex + C2*Sindex*Sindex;
+    float norm = exp(lnorm);
+    return norm;
+}
+
 # if (0)
 // create a vector containing only the unique entries in the input vector
Index: /branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.h
===================================================================
--- /branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.h	(revision 36020)
+++ /branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.h	(revision 36021)
@@ -64,3 +64,6 @@
 float pmModelCP_SersicSubpix (float dx, float dy, float Rxx, float Rxy, float Ryy, float Sindex, int Nsub);
 
+float pmSersicKappa (float Sindex);
+float pmSersicNorm (float Sindex);
+
 # endif
Index: /branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.c
===================================================================
--- /branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.c	(revision 36020)
+++ /branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.c	(revision 36021)
@@ -261,9 +261,9 @@
         // 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;
+	// constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
         if (params->n == 7) {
+	    nParams = params->n;
+	} else {
 	    nParams = params->n - 1;
-	} else {
-	    nParams = params->n - 2;
 	    constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
 	}
@@ -488,2 +488,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;
+}
+
Index: /branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.h
===================================================================
--- /branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.h	(revision 36020)
+++ /branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.h	(revision 36021)
@@ -98,4 +98,6 @@
 bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize);
 
+bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize);
+
 /// @}
 # endif /* PM_PCM_DATA_H */