Index: trunk/psModules/src/objects/models/pmModel_DEV.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_DEV.c	(revision 31670)
+++ trunk/psModules/src/objects/models/pmModel_DEV.c	(revision 32347)
@@ -89,4 +89,6 @@
 static bool limitsApply = true;         // Apply limits?
 
+# include "pmModel_SERSIC.CP.h"
+
 psF32 PM_MODEL_FUNC (psVector *deriv,
                      const psVector *params,
@@ -94,8 +96,4 @@
 {
     psF32 *PAR = params->data.F32;
-
-    float index = 0.5 / ALPHA;
-    float bn = 1.9992*index - 0.3271;
-    float Io = exp(bn);
 
     psF32 X  = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS];
@@ -105,8 +103,89 @@
     psF32 z  = (PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y);
 
-    assert (z >= 0);
+    // If the elliptical contour is defined in a valid way, we should never trigger this
+    // assert.  Other models (like PGAUSS) don't use fractional powers, and thus do not have
+    // NaN values for negative values of z
+    psAssert (z >= 0, "do not allow negative z values in model");
+
+    float index = 0.5 / ALPHA;
+    float par7 = ALPHA;
+    float bn = 1.9992*index - 0.3271;
+    float Io = exp(bn);
 
     psF32 f2 = bn*pow(z,ALPHA);
     psF32 f1 = Io*exp(-f2);
+
+    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)
+	psEllipseShape shape;
+
+	shape.sx  = PAR[PM_PAR_SXX];
+	shape.sy  = PAR[PM_PAR_SYY];
+	shape.sxy = PAR[PM_PAR_SXY];
+
+	// for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
+	psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+
+	// 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;
@@ -120,9 +199,9 @@
         psF32 *dPAR = deriv->data.F32;
 
+        dPAR[PM_PAR_SKY]  = +1.0;
+        dPAR[PM_PAR_I0]   = +2.0*f1; // XXX extra damping..
+
         // gradient is infinite for z = 0; saturate at z = 0.01
         psF32 z1 = (z < 0.01) ? z0*bn*ALPHA*pow(0.01,ALPHA - 1.0) : z0*bn*ALPHA*pow(z,ALPHA - 1.0);
-
-        dPAR[PM_PAR_SKY]  = +1.0;
-        dPAR[PM_PAR_I0]   = +2.0*f1;
 
         assert (isfinite(z1));
@@ -223,4 +302,5 @@
 
     // set the shape parameters
+    // XXX adjust this?
     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
       return false;
@@ -246,36 +326,28 @@
 }
 
+// A DeVaucouleur model is equivalent to a Sersic with index = 4.0
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    float z, norm;
     psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
+    // for a non-circular DeVaucouleur, the flux of the Rmajor equivalent is scaled by the AspectRatio
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
-    norm = 0.0;
-
-    # define DZ 0.25
-
-    float f0 = 1.0;
-    float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-	f1 = exp(-pow(z,ALPHA));
-        z += DZ;
-	f2 = exp(-pow(z,ALPHA));
-        norm += f0 + 4*f1 + f2;
-        f0 = f2;
-    }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    float AspectRatio = axes.minor / axes.major;
+
+    float index = 4.0;
+    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);
@@ -297,6 +369,6 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 
Index: trunk/psModules/src/objects/models/pmModel_EXP.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_EXP.c	(revision 31670)
+++ trunk/psModules/src/objects/models/pmModel_EXP.c	(revision 32347)
@@ -81,4 +81,6 @@
 static bool limitsApply = true;         // Apply limits?
 
+# include "pmModel_SERSIC.CP.h"
+
 psF32 PM_MODEL_FUNC (psVector *deriv,
                      const psVector *params,
@@ -93,13 +95,89 @@
     psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
 
-    // XXX if the elliptical contour is defined in valid way, this step should not be required.
-    // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values
-    // for negative values of z
-    // XXX use an assert here to force the elliptical parameters to be correctly determined
-    // if (z < 0) z = 0;
-    assert (z >= 0);
-
-    psF32 f2 = sqrt(z);
-    psF32 f1 = exp(-f2);
+    // If the elliptical contour is defined in a valid way, we should never trigger this
+    // assert.  Other models (like PGAUSS) don't use fractional powers, and thus do not have
+    // NaN values for negative values of z
+    psAssert (z >= 0, "do not allow negative z values in model");
+
+    float index = 1.0;
+    float par7 = 0.5;
+    float bn = 1.9992*index - 0.3271;
+    float Io = exp(bn);
+
+    psF32 f2 = bn*sqrt(z);
+    psF32 f1 = Io*exp(-f2);
+
+    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)
+	psEllipseShape shape;
+
+	shape.sx  = PAR[PM_PAR_SXX];
+	shape.sy  = PAR[PM_PAR_SYY];
+	shape.sxy = PAR[PM_PAR_SXY];
+
+	// for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
+	psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+
+	// 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;
@@ -118,5 +196,5 @@
         // gradient is infinite for z = 0; saturate at z = 0.01
 	// z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i]
-        psF32 z1 = (z < 0.01) ? 0.5*z0/sqrt(0.01) : 0.5*z0/sqrt(z);
+        psF32 z1 = (z < 0.01) ? 0.5*bn*z0/sqrt(0.01) : 0.5*bn*z0/sqrt(z);
 
 	// XXX dampen SXX and SYY as in GAUSS?
@@ -216,4 +294,5 @@
 
     // set the shape parameters
+    // XXX adjust this?
     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
       return false;
@@ -233,36 +312,28 @@
 }
 
+// An exponential model is equivalent to a Sersic with index = 1.0
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    float z, norm;
     psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
+    // for a non-circular Exponential, the flux of the Rmajor equivalent is scaled by the AspectRatio
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z) = exp(-sqrt(z)) [0 to infinity]
-    norm = 0.0;
-
-    # define DZ 0.25
-
-    float f0 = 1.0;
-    float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        f1 = exp(-sqrt(z));
-        z += DZ;
-        f2 = exp(-sqrt(z));
-        norm += f0 + 4*f1 + f2;
-        f0 = f2;
-    }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    float AspectRatio = axes.minor / axes.major;
+
+    float index = 1.0;
+    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);
@@ -284,6 +355,6 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 
Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 31670)
+++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 32347)
@@ -217,7 +217,8 @@
 }
 
-psF64 PM_MODEL_FLUX(const psVector *params)
-{
-    float norm, z;
+// integrate the model to get the full flux
+psF64 PM_MODEL_FLUX (const psVector *params)
+{
+    float z, norm;
     psEllipseShape shape;
 
@@ -228,25 +229,27 @@
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
+    float AspectRatio = axes.minor / axes.major;
+
+    // flux = 2 \pi \int f(r) r dr
     norm = 0.0;
 
-    # define DZ 0.25
-
-    float f0 = 1.0;
+    # define DR 0.25
+
+    // f = f(r) * r
+    float f0 = 0.0;
     float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        f1 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
-        z += DZ;
-        f2 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
+    for (float r = DR; r < 150; r += DR) {
+	z = 0.5 * PS_SQR(r / axes.major);
+        f1 = r / (1 + z + z*z/2.0 + z*z*z/6.0);
+        r += DR;
+	z = 0.5 * PS_SQR(r / axes.major);
+        f2 = r / (1 + z + z*z/2.0 + z*z*z/6.0);
         norm += f0 + 4*f1 + f2;
         f0 = f2;
     }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    norm *= DR / 3.0;
+
+    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
 
     return(Flux);
Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_PS1_V1.c	(revision 31670)
+++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c	(revision 32347)
@@ -14,6 +14,6 @@
    * PM_PAR_XPOS 2  - X center of object
    * PM_PAR_YPOS 3  - Y center of object
-   * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
-   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
+   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
+   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
    * PM_PAR_SXY 6   - X*Y term of elliptical contour
    * PM_PAR_7   7   - amplitude of the linear component (k)
@@ -239,4 +239,5 @@
 }
 
+// integrate the model to get the full flux
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
@@ -250,25 +251,27 @@
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
+    float AspectRatio = axes.minor / axes.major;
+
+    // flux = 2 \pi \int f(r) r dr
     norm = 0.0;
 
-    # define DZ 0.25
-
-    float f0 = 1.0;
+    # define DR 0.25
+
+    // f = f(r) * r
+    float f0 = 0.0;
     float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
-        z += DZ;
-        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+    for (float r = DR; r < 150; r += DR) {
+	z = 0.5 * PS_SQR(r / axes.major);
+        f1 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+        r += DR;
+	z = 0.5 * PS_SQR(r / axes.major);
+        f2 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
         norm += f0 + 4*f1 + f2;
         f0 = f2;
     }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    norm *= DR / 3.0;
+
+    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
 
     return(Flux);
Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 31670)
+++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 32347)
@@ -240,4 +240,5 @@
 }
 
+// integrate the model to get the full flux
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
@@ -251,25 +252,27 @@
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
+    float AspectRatio = axes.minor / axes.major;
+
+    // flux = 2 \pi \int f(r) r dr
     norm = 0.0;
 
-    # define DZ 0.25
-
-    float f0 = 1.0;
+    # define DR 0.25
+
+    // f = f(r) * r
+    float f0 = 0.0;
     float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
-        z += DZ;
-        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+    for (float r = DR; r < 150; r += DR) {
+	z = 0.5 * PS_SQR(r / axes.major);
+        f1 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+        r += DR;
+	z = 0.5 * PS_SQR(r / axes.major);
+        f2 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
         norm += f0 + 4*f1 + f2;
         f0 = f2;
     }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    norm *= DR / 3.0;
+
+    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
 
     return(Flux);
Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_RGAUSS.c	(revision 31670)
+++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c	(revision 32347)
@@ -229,7 +229,8 @@
 }
 
+// integrate the model to get the full flux
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    float norm, z;
+    float z, norm;
     psEllipseShape shape;
 
@@ -240,25 +241,27 @@
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
+    float AspectRatio = axes.minor / axes.major;
+
+    // flux = 2 \pi \int f(r) r dr
     norm = 0.0;
 
-    # define DZ 0.25
-
-    float f0 = 1.0;
+    # define DR 0.25
+
+    // f = f(r) * r
+    float f0 = 0.0;
     float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        f1 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
-        z += DZ;
-        f2 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
+    for (float r = DR; r < 150; r += DR) {
+	z = 0.5 * PS_SQR(r / axes.major);
+        f1 = r / (1 + z + pow(z, PAR[PM_PAR_7]));
+        r += DR;
+	z = 0.5 * PS_SQR(r / axes.major);
+        f2 = r / (1 + z + pow(z, PAR[PM_PAR_7]));
         norm += f0 + 4*f1 + f2;
         f0 = f2;
     }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    norm *= DR / 3.0;
+
+    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
 
     return(Flux);
Index: trunk/psModules/src/objects/models/pmModel_SERSIC.CP.h
===================================================================
--- trunk/psModules/src/objects/models/pmModel_SERSIC.CP.h	(revision 32347)
+++ trunk/psModules/src/objects/models/pmModel_SERSIC.CP.h	(revision 32347)
@@ -0,0 +1,30 @@
+static int centralPixelNX = 29;
+static int centralPixelNY = 12;
+static float centralPixelXo = 2.0;
+static float centralPixelYo = 0.5;
+static float centralPixeldX = 0.5;
+static float centralPixeldY = 0.5;
+
+static float centralPixel[12][29] = {
+{  0.0521,  0.0336,  0.0235,  0.0173,  0.0133,  0.0105,  0.0085,  0.0070,  0.0059,  0.0051,  0.0044,  0.0038,  0.0033,  0.0030,  0.0026,  0.0024,  0.0021,  0.0019,  0.0018,  0.0016,  0.0015,  0.0014,  0.0013,  0.0012,  0.0011,  0.0010,  0.0010,  0.0009,  0.0008, }, 
+{  0.0817,  0.0556,  0.0402,  0.0304,  0.0238,  0.0191,  0.0157,  0.0131,  0.0111,  0.0096,  0.0083,  0.0073,  0.0064,  0.0057,  0.0051,  0.0046,  0.0042,  0.0038,  0.0035,  0.0032,  0.0029,  0.0027,  0.0025,  0.0023,  0.0022,  0.0020,  0.0019,  0.0018,  0.0017, }, 
+{  0.1096,  0.0785,  0.0591,  0.0462,  0.0371,  0.0305,  0.0255,  0.0217,  0.0187,  0.0162,  0.0143,  0.0126,  0.0113,  0.0101,  0.0091,  0.0083,  0.0076,  0.0069,  0.0064,  0.0059,  0.0054,  0.0050,  0.0047,  0.0044,  0.0041,  0.0038,  0.0036,  0.0034,  0.0032, }, 
+{  0.1339,  0.0998,  0.0777,  0.0624,  0.0513,  0.0431,  0.0367,  0.0317,  0.0277,  0.0244,  0.0217,  0.0194,  0.0175,  0.0158,  0.0144,  0.0132,  0.0121,  0.0112,  0.0103,  0.0096,  0.0089,  0.0083,  0.0078,  0.0073,  0.0069,  0.0065,  0.0061,  0.0058,  0.0055, }, 
+{  0.1547,  0.1190,  0.0951,  0.0781,  0.0655,  0.0559,  0.0484,  0.0423,  0.0374,  0.0334,  0.0300,  0.0271,  0.0246,  0.0225,  0.0206,  0.0190,  0.0176,  0.0163,  0.0152,  0.0142,  0.0133,  0.0125,  0.0117,  0.0110,  0.0104,  0.0099,  0.0093,  0.0089,  0.0084, }, 
+{  0.1726,  0.1361,  0.1110,  0.0928,  0.0791,  0.0685,  0.0600,  0.0531,  0.0475,  0.0427,  0.0387,  0.0353,  0.0323,  0.0297,  0.0275,  0.0255,  0.0237,  0.0221,  0.0207,  0.0195,  0.0183,  0.0172,  0.0163,  0.0154,  0.0146,  0.0139,  0.0132,  0.0126,  0.0120, }, 
+{  0.1881,  0.1513,  0.1255,  0.1065,  0.0919,  0.0805,  0.0713,  0.0637,  0.0574,  0.0521,  0.0476,  0.0437,  0.0403,  0.0373,  0.0347,  0.0323,  0.0303,  0.0284,  0.0267,  0.0252,  0.0238,  0.0225,  0.0214,  0.0203,  0.0193,  0.0184,  0.0176,  0.0168,  0.0161, }, 
+{  0.2017,  0.1649,  0.1387,  0.1191,  0.1040,  0.0919,  0.0821,  0.0740,  0.0672,  0.0614,  0.0564,  0.0521,  0.0483,  0.0450,  0.0420,  0.0394,  0.0370,  0.0349,  0.0329,  0.0312,  0.0296,  0.0281,  0.0268,  0.0255,  0.0244,  0.0233,  0.0223,  0.0214,  0.0205, }, 
+{  0.2138,  0.1771,  0.1507,  0.1308,  0.1152,  0.1027,  0.0924,  0.0839,  0.0767,  0.0705,  0.0651,  0.0604,  0.0563,  0.0527,  0.0494,  0.0465,  0.0439,  0.0415,  0.0393,  0.0374,  0.0356,  0.0339,  0.0324,  0.0310,  0.0296,  0.0284,  0.0273,  0.0262,  0.0252, }, 
+{  0.2246,  0.1883,  0.1618,  0.1416,  0.1257,  0.1128,  0.1022,  0.0933,  0.0857,  0.0792,  0.0735,  0.0686,  0.0642,  0.0603,  0.0568,  0.0536,  0.0508,  0.0482,  0.0458,  0.0436,  0.0416,  0.0398,  0.0381,  0.0365,  0.0351,  0.0337,  0.0324,  0.0313,  0.0301, }, 
+{  0.2347,  0.1986,  0.1721,  0.1517,  0.1355,  0.1224,  0.1115,  0.1023,  0.0945,  0.0877,  0.0817,  0.0765,  0.0719,  0.0677,  0.0640,  0.0606,  0.0576,  0.0548,  0.0522,  0.0499,  0.0477,  0.0458,  0.0439,  0.0422,  0.0406,  0.0391,  0.0377,  0.0364,  0.0352, }, 
+{  0.2448,  0.2086,  0.1819,  0.1614,  0.1450,  0.1316,  0.1204,  0.1110,  0.1029,  0.0958,  0.0897,  0.0842,  0.0793,  0.0750,  0.0711,  0.0675,  0.0643,  0.0613,  0.0586,  0.0561,  0.0538,  0.0517,  0.0497,  0.0479,  0.0462,  0.0445,  0.0430,  0.0416,  0.0403, }, 
+}; 
+
+static float interpolatePixels (float V00, float V10, float V01, float V11, float dx, float dy) {
+    // bilinear interpolation
+    float rx = 1.0 - dx;
+    float ry = 1.0 - dy;
+    float value = V00*rx*ry + V10*dx*ry + V01*rx*dy + V11*dx*dy;
+    return value;
+}
+
Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_SERSIC.c	(revision 31670)
+++ trunk/psModules/src/objects/models/pmModel_SERSIC.c	(revision 32347)
@@ -13,8 +13,14 @@
    * PM_PAR_XPOS 2  - X center of object
    * PM_PAR_YPOS 3  - Y center of object
-   * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
-   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
+   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
+   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
    * PM_PAR_SXY 6   - X*Y term of elliptical contour
    * PM_PAR_7   7   - normalized sersic parameter
+
+   * 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:
+     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)
 
    note that a standard sersic model uses exp(-K*(z^(1/2n) - 1).  the additional elements (K,
@@ -85,4 +91,6 @@
 static bool limitsApply = true;         // Apply limits?
 
+# include "pmModel_SERSIC.CP.h"
+
 psF32 PM_MODEL_FUNC (psVector *deriv,
                      const psVector *params,
@@ -97,17 +105,89 @@
     psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
 
-    // XXX if the elliptical contour is defined in valid way, this step should not be required.
-    // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values
-    // for negative values of z
-    // XXX use an assert here to force the elliptical parameters to be correctly determined
-    // if (z < 0) z = 0;
-    assert (z >= 0);
+    // If the elliptical contour is defined in a valid way, we should never trigger this
+    // assert.  Other models (like PGAUSS) don't use fractional powers, and thus do not have
+    // NaN values for negative values of z
+    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,PAR[PM_PAR_7]);
+    psF32 f2 = bn*pow(z,par7);
     psF32 f1 = Io*exp(-f2);
+
+    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)
+	psEllipseShape shape;
+
+	shape.sx  = PAR[PM_PAR_SXX];
+	shape.sy  = PAR[PM_PAR_SYY];
+	shape.sxy = PAR[PM_PAR_SXY];
+
+	// for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
+	psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+
+	// 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;
@@ -121,10 +201,11 @@
         psF32 *dPAR = deriv->data.F32;
 
-        // gradient is infinite for z = 0; saturate at z = 0.01
-        psF32 z1 = (z < 0.01) ? z0*bn*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*bn*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0);
-
         dPAR[PM_PAR_SKY]  = +1.0;
         dPAR[PM_PAR_I0]   = +f1;
-        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -z0*f2*log(z);
+
+        // 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;
 
@@ -269,8 +350,4 @@
     float Io = exp(0.5*bn);
 
-    // XXX do we need this factor of sqrt(2)?
-    // float Sxx = PS_MAX(0.5, M_SQRT2*shape.sx);
-    // float Syy = PS_MAX(0.5, M_SQRT2*shape.sy);
-
     float Sxx = PS_MAX(0.5, shape.sx);
     float Syy = PS_MAX(0.5, shape.sy);
@@ -294,37 +371,28 @@
 }
 
+// A sersic model has an integral flux which can be analytically represented
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    float z, norm;
     psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
+    // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
-    norm = 0.0;
-
-    # define DZ 0.25
-
-    float f0 = 1.0;
-    float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
-	f1 = exp(-pow(z,PAR[PM_PAR_7]));
-        z += DZ;
-	f2 = exp(-pow(z,PAR[PM_PAR_7]));
-        norm += f0 + 4*f1 + f2;
-        f0 = f2;
-    }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    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);
@@ -346,6 +414,6 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];