Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 26916)
+++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 27531)
@@ -236,5 +236,5 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -271,19 +271,22 @@
     z0 = 0.0;
 
-    // perform a type of bisection to find the value
-    float f0 = 1.0 / (1 + z0 + z0*z0/2.0 + z0*z0*z0/6.0);
-    float f1 = 1.0 / (1 + z1 + z1*z1/2.0 + z1*z1*z1/6.0);
-    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
-        z = 0.5*(z0 + z1);
-        f = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
-        if (f > limit) {
-            z0 = z;
-            f0 = f;
-        } else {
-            z1 = z;
-            f1 = f;
-        }
-        Nstep ++;
-    }
+    // starting guess:
+    z = 0.5*(z0 + z1);
+    float dz = 1.0;
+
+    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
+	// use Newton-Raphson to minimize f(z) - limit = 0
+	float dqdz = (1.0 + z + z*z/2.0);
+	float q = (dqdz + z*z*z/6.0);
+
+	float f = 1.0 / q;
+	float dfdz = -dqdz * f / q;
+
+	dz = (f - limit) / dfdz;
+
+	// fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
+	z -= dz;
+    }
+
     psF64 radius = sigma * sqrt (2.0 * z);
 
Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_PS1_V1.c	(revision 26916)
+++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c	(revision 27531)
@@ -267,5 +267,5 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -273,16 +273,8 @@
     psF32 *PAR = params->data.F32;
 
-    if (flux <= 0) {
-        return 1.0;
-    }
-    if (PAR[PM_PAR_I0] <= 0) {
-        return 1.0;
-    }
-    if (flux >= PAR[PM_PAR_I0]) {
-        return 1.0;
-    }
-    if (PAR[PM_PAR_7] == 0.0) {
-        return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
-    }
+    if (flux <= 0) return 1.0;
+    if (PAR[PM_PAR_I0] <= 0) return 1.0;
+    if (flux >= PAR[PM_PAR_I0]) return 1.0;
+    if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
 
     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
@@ -305,18 +297,20 @@
     if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
 
-    // perform a type of bisection to find the value
-    float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, ALPHA));
-    float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, ALPHA));
-    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
-        z = 0.5*(z0 + z1);
-        f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
-        if (f > limit) {
-            z0 = z;
-            f0 = f;
-        } else {
-            z1 = z;
-            f1 = f;
-        }
-        Nstep ++;
+    // starting guess:
+    z = 0.5*(z0 + z1);
+    float dz = 1.0;
+
+    // use Newton-Raphson to minimize f(z) - limit = 0
+    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
+	float q = (1.0 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+	float dqdz = (PAR[PM_PAR_7] + ALPHA*pow(z, ALPHA - 1.0));
+
+	float f = 1.0 / q;
+	float dfdz = -dqdz * f / q;
+
+	dz = (f - limit) / dfdz;
+
+	// fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
+	z -= dz;
     }
     psF64 radius = sigma * sqrt (2.0 * z);
@@ -350,5 +344,4 @@
     // convert to shape terms (SXX,SYY,SXY)
     if (!pmPSF_FitToModel (out, 0.1)) {
-        // psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -475,5 +468,4 @@
 }
 
-
 # undef PM_MODEL_FUNC
 # undef PM_MODEL_FLUX
Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 26916)
+++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 27531)
@@ -39,4 +39,7 @@
 # define PM_MODEL_SET_LIMITS      pmModelSetLimits_QGAUSS
 
+# define ALPHA   2.250
+# define ALPHA_M 1.250
+
 // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
 // 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
@@ -44,19 +47,21 @@
 
 // Lax parameter limits
-static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.1 };
+static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
 
 // Moderate parameter limits
-static float *paramsMinModerate = paramsMinLax;
-static float *paramsMaxModerate = paramsMaxLax;
+// Tolerate a small divot (k < 0)
+static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 };
+static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
 
 // Strict parameter limits
-static float *paramsMinStrict = paramsMinLax;
-static float *paramsMaxStrict = paramsMaxLax;
+// k = PAR_7 < 0 is very undesirable (big divot in the middle)
+static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
+static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
 
 // Parameter limits to use
 static float *paramsMinUse = paramsMinLax;
 static float *paramsMaxUse = paramsMaxLax;
-static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 };
+static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
 
 static bool limitsApply = true;         // Apply limits?
@@ -81,5 +86,5 @@
     assert (z >= 0);
 
-    psF32 zp = pow(z,1.25);
+    psF32 zp = pow(z,ALPHA_M);
     psF32 r  = 1.0 / (1 + PAR[PM_PAR_7]*z + z*zp);
 
@@ -92,5 +97,5 @@
         // note difference from a pure gaussian: q = params->data.F32[PM_PAR_I0]*r
         psF32 t = r1*r;
-        psF32 q = t*(PAR[PM_PAR_7] + 2.25*zp);
+        psF32 q = t*(PAR[PM_PAR_7] + ALPHA*zp);
 
         dPAR[PM_PAR_SKY]  = +1.0;
@@ -246,7 +251,7 @@
     float f1, f2;
     for (z = DZ; z < 50; z += DZ) {
-        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
+        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, 2.25));
+        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
         norm += f0 + 4*f1 + f2;
         f0 = f2;
@@ -263,5 +268,5 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -269,10 +274,8 @@
     psF32 *PAR = params->data.F32;
 
-    if (flux <= 0)
-        return (1.0);
-    if (PAR[PM_PAR_I0] <= 0)
-        return (1.0);
-    if (flux >= PAR[PM_PAR_I0])
-        return (1.0);
+    if (flux <= 0) return 1.0;
+    if (PAR[PM_PAR_I0] <= 0) return 1.0;
+    if (flux >= PAR[PM_PAR_I0]) return 1.0;
+    if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
 
     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
@@ -290,25 +293,25 @@
 
     // choose a z value guaranteed to be beyond our limit
-    float z0 = pow((1.0 / limit), (1.0 / 2.25));
-    psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_7]);
-    float z1 = (1.0 / limit) / PAR[PM_PAR_7];
+    float z0 = 0.0;
+    float z1 = pow((1.0 / limit), (1.0 / ALPHA));
     psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
-    z1 = PS_MAX (z0, z1);
-    z0 = 0.0;
-
-    // perform a type of bisection to find the value
-    float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, 2.25));
-    float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, 2.25));
-    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
-        z = 0.5*(z0 + z1);
-        f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
-        if (f > limit) {
-            z0 = z;
-            f0 = f;
-        } else {
-            z1 = z;
-            f1 = f;
-        }
-        Nstep ++;
+    if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
+
+    // starting guess:
+    z = 0.5*(z0 + z1);
+    float dz = 1.0;
+
+    // use Newton-Raphson to minimize f(z) - limit = 0
+    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
+	float q = (1.0 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+	float dqdz = (PAR[PM_PAR_7] + ALPHA*pow(z, ALPHA - 1.0));
+
+	float f = 1.0 / q;
+	float dfdz = -dqdz * f / q;
+
+	dz = (f - limit) / dfdz;
+
+	// fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
+	z -= dz;
     }
     psF64 radius = sigma * sqrt (2.0 * z);
@@ -475,2 +478,4 @@
 # undef PM_MODEL_FIT_STATUS
 # undef PM_MODEL_SET_LIMITS
+# undef ALPHA
+# undef ALPHA_M
Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_RGAUSS.c	(revision 26916)
+++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c	(revision 27531)
@@ -257,5 +257,5 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -291,19 +291,22 @@
     z0 = 0.0;
 
-    // perform a type of bisection to find the value
-    float f0 = 1.0 / (1 + z0 + pow(z0, PAR[PM_PAR_7]));
-    float f1 = 1.0 / (1 + z1 + pow(z1, PAR[PM_PAR_7]));
-    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
-        z = 0.5*(z0 + z1);
-        f = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
-        if (f > limit) {
-            z0 = z;
-            f0 = f;
-        } else {
-            z1 = z;
-            f1 = f;
-        }
-        Nstep ++;
-    }
+    // starting guess:
+    z = 0.5*(z0 + z1);
+    float dz = 1.0;
+
+    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
+	// use Newton-Raphson to minimize f(z) - limit = 0
+	float q = (1 + z + pow(z,PAR[PM_PAR_7]));
+	float dqdz = (1.0 + PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0));
+
+	float f = 1.0 / q;
+	float dfdz = -dqdz * f / q;
+
+	dz = (f - limit) / dfdz;
+
+	// fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
+	z -= dz;
+    }
+
     psF64 radius = sigma * sqrt (2.0 * z);