Index: branches/eam_branches/ipp-20130711/psModules/src/objects/Makefile.am
===================================================================
--- branches/eam_branches/ipp-20130711/psModules/src/objects/Makefile.am	(revision 35948)
+++ branches/eam_branches/ipp-20130711/psModules/src/objects/Makefile.am	(revision 35961)
@@ -20,4 +20,5 @@
 	pmModelClass.c \
 	pmModelUtils.c \
+	pmModel_CentralPixel.c \
 	pmSource.c \
 	pmPhotObj.c \
@@ -97,4 +98,5 @@
 	pmModelClass.h \
 	pmModelUtils.h \
+	pmModel_CentralPixel.h \
 	pmSource.h \
 	pmPhotObj.h \
@@ -110,5 +112,5 @@
 	pmSourceOutputs.h \
 	pmSourceIO.h \
-	pmSourceSatstar.h \ 
+	pmSourceSatstar.h \
 	pmSourcePlots.h \
 	pmSourceVisual.h \
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 35948)
+++ branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c	(revision 35961)
@@ -16,5 +16,4 @@
    * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
    * PM_PAR_SXY 6   - X*Y term of elliptical contour
-   * PM_PAR_7   7   - normalized dev parameter
 
    note that a standard dev model uses exp(-K*(z^(1/2n) - 1).  the additional elements (K,
@@ -86,9 +85,9 @@
 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, 3.0, 3.0, 0.5 };
 
 static bool limitsApply = true;         // Apply limits?
 
-# include "pmModel_SERSIC.CP.h"
+// # include "pmModel_SERSIC.CP.h"
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -109,85 +108,52 @@
     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);
-
+    // for DEV, we can hard-wire kappa(4):
+    // float index = 4.0;
+    float kappa = 7.670628;
+
+    // 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;
+
+    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);
     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);
+      // subdivide the central 2,3,4 pixels by Nx,Ny 
+      float Npix = 0.0;
+      float Fpix = 0.0;
+      float Xpix = floor(pixcoord->data.F32[0]) - PAR[PM_PAR_XPOS];
+      float Ypix = floor(pixcoord->data.F32[1]) - PAR[PM_PAR_YPOS];
+      # define DELTA 0.02
+      for (float ix = 0.1; ix <= 0.9; ix += DELTA) {
+	for (float iy = 0.1; iy <= 0.9; iy += DELTA) {
+	  psF32 X  = Xpix + ix;
+	  psF32 Y  = Ypix + iy;
+	  psF32 px = X / PAR[PM_PAR_SXX];
+	  psF32 py = Y / PAR[PM_PAR_SYY];
+	  psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
+	  
+	  // sqrt(z) is r
+	  float q = kappa*pow(z,ALPHA);
+	  psF32 f0 = exp(-q);
+	  
+	  psF32 f1 = PAR[PM_PAR_I0]*f0;
+	  psF32 fx = PAR[PM_PAR_SKY] + f1;
+	  Fpix += fx;
+	  Npix += 1.0;
 	}
-	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));
-	}
+      }
+      f = Fpix / Npix;
     }   
-
-    psF32 z0 = PAR[PM_PAR_I0]*f1;
-    psF32 f0 = PAR[PM_PAR_SKY] + z0;
-
-    assert (isfinite(f2));
-    assert (isfinite(f1));
-    assert (isfinite(z0));
-    assert (isfinite(f0));
 
     if (deriv != NULL) {
@@ -195,18 +161,24 @@
 
         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);
-
-        assert (isfinite(z1));
-
-        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];
-        dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
-        dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
-    }
-    return (f0);
+        dPAR[PM_PAR_I0]   = +f0;
+
+	if (z > 0.01) {
+	  float z1 = f1*kappa*ALPHA*pow(z,ALPHA-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;
+	} else {
+	  // gradient -> 0 for z -> 0, but has undef form
+	  float z1 = f1*kappa*ALPHA*pow(z,ALPHA);
+	  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;
+	}
+    }
+    return (f);
 }
 
@@ -302,14 +274,8 @@
     }
 
-    // the normalization is modified by the slope
-    float index = 0.5 / ALPHA;
-    float bn = 1.9992*index - 0.3271;
-    float Io = exp(0.5*bn);
-
     // set the model normalization
     if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
       return false;
     }
-    PAR[PM_PAR_I0] /= Io;
 
     // set the model position
@@ -328,17 +294,9 @@
     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 = 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);
+
+    float norm = 0.00168012;
+    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
+
+    return(flux);
 }
 
@@ -359,6 +317,9 @@
     pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
-    // f = Io exp(-z^n) -> z^n = ln(Io/f)
-    psF64 zn = log(PAR[PM_PAR_I0] / flux);
+    // static value for DEV:
+    float kappa = 7.670628;
+
+    // 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);
 
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 35948)
+++ branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_EXP.c	(revision 35961)
@@ -82,5 +82,10 @@
 static bool limitsApply = true;         // Apply limits?
 
-# include "pmModel_SERSIC.CP.h"
+// # include "pmModel_SERSIC.CP.h"
+
+// the problems I'm having with the SERSIC-like functions are:
+// 1) making sure I have the right functional form so that PAR[SXX,etc] represent R_eff (half-light radius)
+// 2) getting the central pixel right
+// 3) getting the derivaties right.
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -101,85 +106,48 @@
     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);
-
+    // for EXP, we can hard-wire kappa(1):
+    // float index = 1.0;
+    float kappa = 1.70056;
+
+    // sqrt(z) is r
+    float q = kappa*sqrt(z);
+    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
     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);
+      // subdivide the central 2,3,4 pixels by Nx,Ny 
+      float Npix = 0.0;
+      float Fpix = 0.0;
+      float Xpix = floor(pixcoord->data.F32[0]) - PAR[PM_PAR_XPOS];
+      float Ypix = floor(pixcoord->data.F32[1]) - PAR[PM_PAR_YPOS];
+      for (float ix = 0.1; ix < 1.0; ix += 0.2) {
+	for (float iy = 0.1; iy < 1.0; iy += 0.2) {
+	  psF32 X  = Xpix + ix;
+	  psF32 Y  = Ypix + iy;
+	  psF32 px = X / PAR[PM_PAR_SXX];
+	  psF32 py = Y / PAR[PM_PAR_SYY];
+	  psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
+	  
+	  // sqrt(z) is r
+	  float q = kappa*sqrt(z);
+	  psF32 f0 = exp(-q);
+	  
+	  psF32 f1 = PAR[PM_PAR_I0]*f0;
+	  psF32 fx = PAR[PM_PAR_SKY] + f1;
+	  Fpix += fx;
+	  Npix += 1.0;
 	}
-	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;
-
-    assert (isfinite(f2));
-    assert (isfinite(f1));
-    assert (isfinite(z0));
-    assert (isfinite(f0));
+      }
+      f = Fpix / Npix;
+    }
 
     if (deriv != NULL) {
@@ -187,18 +155,24 @@
 
         dPAR[PM_PAR_SKY]  = +1.0;
-        dPAR[PM_PAR_I0]   = +f1;
-
-        // 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*bn*z0/sqrt(0.01) : 0.5*bn*z0/sqrt(z);
-
-	// XXX dampen SXX and SYY as in GAUSS?
-        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];
-        dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
-        dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
-    }
-    return (f0);
+        dPAR[PM_PAR_I0]   = +f0;
+
+	if (z > 0.01) {
+	  float z1 = 0.5*f1*kappa/sqrt(z);
+	  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;
+	} else {
+	  // gradient -> 0 for z -> 0, but has undef form
+	  float z1 = 0.5*f1*kappa;
+	  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;
+	}
+    }
+    return (f);
 }
 
@@ -314,17 +288,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 = 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);
+
+    // static value for EXP:
+    float norm = 0.34578; // \int exp(-kappa*sqrt(z)) r dr
+
+    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
+
+    return(flux);
 }
 
@@ -345,6 +313,9 @@
     pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
-    // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f)
-    psF64 zn = log(PAR[PM_PAR_I0] / flux);
+    // static value for EXP:
+    float kappa = 1.70056;
+
+    // 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;
 
@@ -501,2 +472,66 @@
     return;
 }
+
+# if (0)
+void bilin_inter_function () {
+	// 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));
+	}
+}
+# endif
Index: branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.c
===================================================================
--- branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.c	(revision 35948)
+++ branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.c	(revision 35961)
@@ -209,9 +209,11 @@
 
 // XXX for test purposes only:
-# define TEST_IMAGE 1
+# define TEST_IMAGE 0
 # if (TEST_IMAGE)
 static psImage *map = NULL;
 # endif
 
+float pmModelCP_GetFlux_RotSquare (pmModelCP *cp, float dx, float dy, float theta);
+
 float pmModelCP_GetFlux (pmModelCP *cp, float dx, float dy, float theta) {
 
@@ -221,6 +223,11 @@
 # endif
 
-    float flux = pmModelCP_GetFlux_Bresen (cp, dx, dy, theta);
-    
+    // float flux = pmModelCP_GetFlux_Bresen (cp, dx, dy, theta);
+    // float flux = pmModelCP_GetFlux_Old (cp, dx, dy, theta);
+    float flux = pmModelCP_GetFlux_RotSquare (cp, dx, dy, theta);
+    
+    // RotSquare for theta = 0.0 & Bresen give the same answer 
+    // if I count from x[0] <= ix < x[1]
+
 # if (TEST_IMAGE) 
     psFits *fits = psFitsOpen ("map.fits", "w");
@@ -309,5 +316,28 @@
 }
 
-float pmModelCP_GetFlux_RotSqaure (pmModelCP *cp, float dx, float dy, float theta) {
+// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
+# define COMPARE_INDEX(A,B) (y[A] < y[B])
+# define SWAP_INDEX(TYPE,A,B) {				\
+	int tmp;					\
+	if (A != B) {					\
+	    tmp = x[A];					\
+	    x[A] = x[B];				\
+	    x[B] = tmp;					\
+	    tmp = y[A];					\
+	    y[A] = y[B];				\
+	    y[B] = tmp;					\
+	}						\
+    }
+
+bool pmModelCP_SortCorners (int *x, int *y, int Npar) {
+
+    if (Npar < 2) return true;
+
+    // sort the vector set by the radius
+    PSSORT (Npar, COMPARE_INDEX, SWAP_INDEX, NONE);
+    return true;
+}
+
+float pmModelCP_GetFlux_RotSquare (pmModelCP *cp, float dx, float dy, float theta) {
 
     // the cp data is defined for the central 3x3 pixels.  we allow dx,dy to have values of
@@ -324,67 +354,146 @@
     float cs = cos(theta*PS_RAD_DEG);
     float sn = sin(theta*PS_RAD_DEG);
-
     float Nsub = 11.0;
-    int Xsub00 = ((dx - 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub;
-    int Ysub00 = ((dx - 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub;
-    int Xsub01 = ((dx - 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub;
-    int Ysub01 = ((dx - 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub;
-    int Xsub10 = ((dx + 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub;
-    int Ysub10 = ((dx + 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub;
-    int Xsub11 = ((dx + 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub;
-    int Ysub11 = ((dx + 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub;
-
-    /* generic rotated square:
-       
-     */
-
-    int Xmin, Xmax, Ymin, Ymax;
-
-    Xmin = PS_MIN(Xsub00,Xsub01);
-    Xmin = PS_MIN(Xsub10,Xmin);
-    Xmin = PS_MIN(Xsub11,Xmin);
-    Xmin = PS_MIN(Xmin, cp->flux->numCols - 1);
-    Xmin = PS_MAX(Xmin, 0);
-    Xmax = PS_MAX(Xsub00,Xsub01);
-    Xmax = PS_MAX(Xsub10,Xmax);
-    Xmax = PS_MAX(Xsub11,Xmax);
-    Xmax = PS_MIN(Xmax, cp->flux->numCols - 1);
-    Xmax = PS_MAX(Xmax, 0);
-    Ymin = PS_MIN(Ysub00,Ysub01);
-    Ymin = PS_MIN(Ysub10,Ymin);
-    Ymin = PS_MIN(Ysub11,Ymin);
-    Ymin = PS_MIN(Ymin, cp->flux->numRows - 1);
-    Ymin = PS_MAX(Ymin, 0);
-    Ymax = PS_MAX(Ysub00,Ysub01);
-    Ymax = PS_MAX(Ysub10,Ymax);
-    Ymax = PS_MAX(Ysub11,Ymax);
-    Ymax = PS_MIN(Ymax, cp->flux->numRows - 1);
-    Ymax = PS_MAX(Ymax, 0);
-
-    // integrate pixels from Xmin,Ymin to Xmax,Ymax, only include pixels contained in the
-    // target pixel
+
+    int Xsub[4], Ysub[4];
+
+    Xsub[0] = ((dx - 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub;
+    Ysub[0] = ((dx - 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub;
+    Xsub[1] = ((dx - 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub;
+    Ysub[1] = ((dx - 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub;
+    Xsub[2] = ((dx + 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub;
+    Ysub[2] = ((dx + 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub;
+    Xsub[3] = ((dx + 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub;
+    Ysub[3] = ((dx + 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub;
+
+    // first, sort the corners in order of the Y coordinate:
+    pmModelCP_SortCorners (Xsub, Ysub, 4);
 
     float flux = 0.0;
-    int   npix = 0;
-    for (int i = Xmin; i < Xmax; i++) {
-	float dX = i / Nsub - 1.5;
-	for (int j = Ymin; j < Ymax; j++) {
-	    float dY = j / Nsub - 1.5;
-
-	    float Xim =  dX*cs + dY*sn;
-	    if (Xim < (dx - 0.5)) continue;
-	    if (Xim > (dx + 0.5)) continue;
-
-	    float Yim = -dX*sn + dY*cs;
-	    if (Yim < (dy - 0.5)) continue;
-	    if (Yim > (dy + 0.5)) continue;
-
-	    flux += cp->flux->data.F32[j][i];
-	    npix ++;
-	}
-    }
-	   
-    float normFlux = flux / npix;
-    return normFlux;
+    float npix = 0.0;
+
+    // if (Ysub[0] == Ysub[1]), we have a simple square
+    if (Ysub[0] == Ysub[1]) {
+	psAssert (Ysub[2] == Ysub[3], "not square?");
+	int Xmin = PS_MIN(Xsub[0], Xsub[1]);
+	int Xmax = PS_MAX(Xsub[0], Xsub[1]);
+	for (int iy = Ysub[0]; iy < Ysub[3]; iy++) {
+	    for (int ix = Xmin; ix < Xmax; ix++) {
+		flux += cp->flux->data.F32[iy][ix];
+		npix += 1.0;
+# if (TEST_IMAGE) 
+		fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
+		map->data.S32[iy][ix] ++;
+# endif
+	    }
+	}
+	float normFlux = flux / npix;
+	return normFlux;
+    }
+    
+    // second case: Xsub[1] > Xsub[2]:
+    if (Xsub[1] > Xsub[2]) {
+	float dYdXp, dYdXm;
+	// first segment, Ysub[0] to Ysub[1]:
+	dYdXp = (Ysub[1] - Ysub[0]) / (float) (Xsub[1] - Xsub[0]);
+	dYdXm = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]);
+	for (int iy = Ysub[0]; iy < Ysub[1]; iy++) {
+	    int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0];
+	    int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0];
+	    for (int ix = Xs; ix < Xe; ix ++) {
+		flux += cp->flux->data.F32[iy][ix];
+		npix += 1.0;
+# if (TEST_IMAGE) 
+		fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
+		map->data.S32[iy][ix] ++;
+# endif
+	    }
+	}
+	// 2nd segment, Ysub[1] to Ysub[2]:
+	dYdXp = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]);
+	dYdXm = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]);
+	for (int iy = Ysub[1]; iy < Ysub[2]; iy++) {
+	    int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0];
+	    int Xe = (iy - Ysub[1]) / dYdXp + Xsub[1];
+	    for (int ix = Xs; ix < Xe; ix ++) {
+		flux += cp->flux->data.F32[iy][ix];
+		npix += 1.0;
+# if (TEST_IMAGE) 
+		fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
+		map->data.S32[iy][ix] ++;
+# endif
+	    }
+	}
+	// first segment, Ysub[0] to Ysub[1]:
+	dYdXp = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]);
+	dYdXm = (Ysub[3] - Ysub[2]) / (float) (Xsub[3] - Xsub[2]);
+	for (int iy = Ysub[2]; iy < Ysub[3]; iy++) {
+	    int Xs = (iy - Ysub[2]) / dYdXm + Xsub[2];
+	    int Xe = (iy - Ysub[1]) / dYdXp + Xsub[1];
+	    for (int ix = Xs; ix < Xe; ix ++) {
+		flux += cp->flux->data.F32[iy][ix];
+		npix += 1.0;
+# if (TEST_IMAGE) 
+		fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
+		map->data.S32[iy][ix] ++;
+# endif
+	    }
+	}
+	float normFlux = flux / npix;
+	return normFlux;
+    }
+
+    // third case: Xsub[1] < Xsub[2]:
+    if (Xsub[2] > Xsub[1]) {
+	// first segment, Ysub[0] to Ysub[1]:
+	float dYdXp, dYdXm;
+	dYdXp = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]);
+	dYdXm = (Ysub[1] - Ysub[0]) / (float) (Xsub[1] - Xsub[0]);
+	for (int iy = Ysub[0]; iy < Ysub[1]; iy++) {
+	    int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0];
+	    int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0];
+	    for (int ix = Xs; ix < Xe; ix ++) {
+		flux += cp->flux->data.F32[iy][ix];
+		npix += 1.0;
+# if (TEST_IMAGE) 
+		fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
+		map->data.S32[iy][ix] ++;
+# endif
+	    }
+	}
+	// 2nd segment, Ysub[1] to Ysub[2]:
+	dYdXp = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]);
+	dYdXm = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]);
+	for (int iy = Ysub[1]; iy < Ysub[2]; iy++) {
+	    int Xs = (iy - Ysub[1]) / dYdXm + Xsub[1];
+	    int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0];
+	    for (int ix = Xs; ix < Xe; ix ++) {
+		flux += cp->flux->data.F32[iy][ix];
+		npix += 1.0;
+# if (TEST_IMAGE) 
+		fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
+		map->data.S32[iy][ix] ++;
+# endif
+	    }
+	}
+	// first segment, Ysub[0] to Ysub[1]:
+	dYdXp = (Ysub[3] - Ysub[2]) / (float) (Xsub[3] - Xsub[2]);
+	dYdXm = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]);
+	for (int iy = Ysub[2]; iy < Ysub[3]; iy++) {
+	    int Xs = (iy - Ysub[1]) / dYdXm + Xsub[1];
+	    int Xe = (iy - Ysub[2]) / dYdXp + Xsub[2];
+	    for (int ix = Xs; ix < Xe; ix ++) {
+		flux += cp->flux->data.F32[iy][ix];
+		npix += 1.0;
+# if (TEST_IMAGE) 
+		fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
+		map->data.S32[iy][ix] ++;
+# endif
+	    }
+	}
+	float normFlux = flux / npix;
+	return normFlux;
+    }
+    myAbort ("impossible case?");
 }
 
@@ -478,5 +587,5 @@
     }
     float normFlux = flux / npix;
-    fprintf (stderr, "bres: %f %f %f\n", flux, (float) npix, normFlux);
+    // fprintf (stderr, "bres: %f %f %f\n", flux, (float) npix, normFlux);
     return normFlux;
 }
@@ -574,5 +683,5 @@
     }
     float normFlux = flux / npix;
-    fprintf (stderr, "full : %f %f %f\n", flux, (float) npix, normFlux);
+    // fprintf (stderr, "full : %f %f %f\n", flux, (float) npix, normFlux);
     return normFlux;
 }
Index: branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCM_MinimizeChisq.c
===================================================================
--- branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCM_MinimizeChisq.c	(revision 35948)
+++ branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCM_MinimizeChisq.c	(revision 35961)
@@ -42,4 +42,9 @@
 #include "pmPCMdata.h"
 
+# define SAVE_IMAGES 0
+# if (SAVE_IMAGES) 
+int psphotSaveImage (psMetadata *header, psImage *image, char *filename);
+# endif
+
 # define FACILITY "psModules.objects"
 
@@ -130,6 +135,30 @@
 	}
 
+	char key[10]; // used for interactive responses
+	bool testValue = false;
+
         // set a new guess for Alpha, Beta, Params
         if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) {
+	    if (min->isInteractive) {
+		fprintf (stdout, "guess failed (singular matrix or NaN values), continue? [Y,n] ");
+		if (!fgets(key, 8, stdin)) {
+		    psWarning("Unable to read option");
+		}
+		switch (key[0]) {
+		  case 'n':
+		  case 'N':
+		    done = true;
+		    break;
+		  case 'y':
+		  case 'Y':
+		  case '\n':
+		    lambda *= 10.0;
+		    continue;
+		  default:
+		    lambda *= 10.0;
+		    continue;
+		}
+		if (done) break;
+	    }
             min->iter ++;
 	    if (min->iter >=  min->maxIter) break;
@@ -138,4 +167,37 @@
         }
 
+	if (min->isInteractive) {
+            p_psVectorPrint(psTraceGetDestination(), Params, "current parameters: ");
+	    fprintf (stdout, "last chisq : %f\n", min->value);
+	    bool getOptions = true;
+	    while (getOptions) {
+		fprintf (stdout, "options: (m)odify, (g)o, (q)uit: ");
+		if (!fgets(key, 8, stdin)) {
+		    psWarning("Unable to read option");
+		}
+		switch (key[0]) {
+		  case 'm':
+		  case 'M':
+		    testValue = TRUE;
+		    fprintf (stdout, "enter (Npar) (value): ");
+		    int Npar = 0;
+		    float value= 0;
+		    fscanf (stdin, "%d %f", &Npar, &value);
+		    Params->data.F32[Npar] = value;
+		    break;
+		  case 'g':
+		  case 'G':
+		  case '\n':
+		    getOptions = false;
+		    break;
+		  default:
+		    done = true;
+		    break;
+		}
+		fprintf (stderr, "foo\n");
+	    }
+	    if (done) break;
+	}
+	    
         // dump some useful info if trace is defined
         if (psTraceGetLevel(FACILITY) >= 6) {
@@ -202,5 +264,5 @@
 	// XXX : Madsen gives suggestion for better use of rho
         // rho is positive if the new chisq is smaller
-        if (rho >= -1e-6) {
+        if (testValue || (rho >= -1e-6)) {
             min->value = Chisq;
             alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
@@ -474,10 +536,18 @@
     // XXX TEST : SAVE IMAGES
 # if (SAVE_IMAGES)
-    psphotSaveImage (NULL, pcm->psf->image, "psf.fits");
-    psphotSaveImage (NULL, pcm->modelFlux, "model.fits");
-    psphotSaveImage (NULL, pcm->modelConvFlux, "modelConv.fits");
-    psphotSaveImage (NULL, source->pixels, "obj.fits");
-    psphotSaveImage (NULL, source->maskObj, "mask.fits");
-    psphotSaveImage (NULL, source->variance, "variance.fits");
+    static int Npass = 0;
+    char name[128]; 
+    snprintf (name, 128, "psf.%03d.fits", Npass); psphotSaveImage (NULL, pcm->psf->image, name);
+    snprintf (name, 128, "mod.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelFlux, name);
+    snprintf (name, 128, "cnv.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelConvFlux, name);
+    snprintf (name, 128, "obj.%03d.fits", Npass); psphotSaveImage (NULL, source->pixels, name);
+    snprintf (name, 128, "msk.%03d.fits", Npass); psphotSaveImage (NULL, source->maskObj, name);
+    snprintf (name, 128, "var.%03d.fits", Npass); psphotSaveImage (NULL, source->variance, name);
+    for (int n = 0; n < pcm->dmodelsFlux->n; n++) {
+        psImage *dmodelConv = pcm->dmodelsConvFlux->data[n];
+	if (!dmodelConv) continue;
+	snprintf (name, 128, "dpar.%01d.%03d.fits", n, Npass); psphotSaveImage (NULL, dmodelConv, name);
+    }
+    Npass ++;
 # endif
 
Index: branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.c
===================================================================
--- branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.c	(revision 35948)
+++ branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.c	(revision 35961)
@@ -242,4 +242,9 @@
         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
         break;
+      case PM_SOURCE_FIT_EXT_AND_SKY:
+        // EXT model fits all params (including sky)
+        nParams = params->n;
+        psVectorInit (constraint->paramMask, 0);
+        break;
       case PM_SOURCE_FIT_INDEX:
         // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
@@ -365,4 +370,9 @@
 	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
 	break;
+      case PM_SOURCE_FIT_EXT_AND_SKY:
+        // EXT model fits all params (including sky)
+        nParams = model->params->n;
+        psVectorInit (pcm->constraint->paramMask, 0);
+        break;
       case PM_SOURCE_FIT_INDEX:
 	// PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
Index: branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitModel.c
===================================================================
--- branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitModel.c	(revision 35948)
+++ branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitModel.c	(revision 35961)
@@ -66,4 +66,5 @@
     opt->gainFactorMode = 0;
     opt->chisqConvergence = true;
+    opt->isInteractive = false;
 
     return opt;
@@ -247,4 +248,5 @@
     myMin->gainFactorMode = options->gainFactorMode;
     myMin->chisqConvergence = options->chisqConvergence;
+    myMin->isInteractive = options->isInteractive;
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
Index: branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitModel.h
===================================================================
--- branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitModel.h	(revision 35948)
+++ branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitModel.h	(revision 35961)
@@ -37,4 +37,5 @@
     int gainFactorMode;
     bool chisqConvergence; 
+    bool isInteractive;
 } pmSourceFitOptions;
 
Index: branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitPCM.c
===================================================================
--- branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitPCM.c	(revision 35948)
+++ branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitPCM.c	(revision 35961)
@@ -68,4 +68,5 @@
     myMin->chisqConvergence = fitOptions->chisqConvergence;
     myMin->gainFactorMode = fitOptions->gainFactorMode;
+    myMin->isInteractive = fitOptions->isInteractive;
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
Index: branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitSet.c
===================================================================
--- branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitSet.c	(revision 35948)
+++ branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitSet.c	(revision 35961)
@@ -570,4 +570,5 @@
     myMin->gainFactorMode = options->gainFactorMode;
     myMin->chisqConvergence = options->chisqConvergence;
+    myMin->isInteractive = options->isInteractive;
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);