Index: trunk/psModules/src/objects/models/pmModel_DEV.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_DEV.c	(revision 35577)
+++ trunk/psModules/src/objects/models/pmModel_DEV.c	(revision 35768)
@@ -123,12 +123,6 @@
 
 	// 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);
+	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
@@ -238,5 +232,6 @@
         // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
         // angle and let f2,f1 fight it out
-        q2 = 0.5*sqrtf(q1);
+	// NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        q2 = 2.0*0.5*sqrtf(q1);
     }
 
@@ -303,6 +298,5 @@
 
     // set the shape parameters
-    // XXX adjust this?
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, true)) {
       return false;
     }
@@ -330,14 +324,8 @@
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX];
-    shape.sy  = PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    // for a non-circular DeVaucouleur, the flux of the Rmajor equivalent is scaled by the AspectRatio
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     float AspectRatio = axes.minor / axes.major;
 
@@ -359,6 +347,4 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
@@ -370,9 +356,6 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX];
-    shape.sy  = PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    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)
@@ -382,5 +365,4 @@
     psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
 	      PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
-
     return (radius);
 }
@@ -407,5 +389,6 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -460,5 +443,6 @@
     // convert to shape terms (SXX,SYY,SXY)
     // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;
Index: trunk/psModules/src/objects/models/pmModel_EXP.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_EXP.c	(revision 35577)
+++ trunk/psModules/src/objects/models/pmModel_EXP.c	(revision 35768)
@@ -115,12 +115,6 @@
 
 	// 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);
+	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
@@ -230,5 +224,6 @@
         // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
         // angle and let f2,f1 fight it out
-        q2 = 0.5*sqrtf(q1);
+	// NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        q2 = 2.0*0.5*sqrtf(q1);
     }
 
@@ -295,6 +290,5 @@
 
     // set the shape parameters
-    // XXX adjust this?
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, true)) {
       return false;
     }
@@ -316,14 +310,8 @@
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX];
-    shape.sy  = PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    // for a non-circular Exponential, the flux of the Rmajor equivalent is scaled by the AspectRatio
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     float AspectRatio = axes.minor / axes.major;
 
@@ -345,6 +333,4 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
@@ -356,9 +342,6 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX];
-    shape.sy  = PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    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)
@@ -392,5 +375,6 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -445,5 +429,6 @@
     // convert to shape terms (SXX,SYY,SXY)
     // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;
Index: trunk/psModules/src/objects/models/pmModel_GAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_GAUSS.c	(revision 35577)
+++ trunk/psModules/src/objects/models/pmModel_GAUSS.c	(revision 35768)
@@ -129,6 +129,7 @@
     float q2 = NAN;
     if (nParam == PM_PAR_SXY) {
-        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
+	// NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
         float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
         q1 = (q1 < 0.0) ? 0.0 : q1;
@@ -200,5 +201,5 @@
 
     // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
       return false;
     }
@@ -219,16 +220,10 @@
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
 
     // Area is equivalent to 2 pi sigma^2
-    // axes ratio < 20
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
 
@@ -242,6 +237,4 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
@@ -253,9 +246,7 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
+
     psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
     psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
@@ -285,12 +276,8 @@
     }
 
-    // the OLD 2D model for SXY actually fitted SXY / (SXX^-2 + SYY^-2); correct here
-    // out[PM_PAR_SXY] = pmPSF_SXYtoModel (out);
-
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    // XXX user-defined value for limit?
-    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]);
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -343,6 +330,6 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;
Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 35577)
+++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 35768)
@@ -129,6 +129,7 @@
     float q2 = NAN;
     if (nParam == PM_PAR_SXY) {
-        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
+	// NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
         float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
         q1 = (q1 < 0.0) ? 0.0 : q1;
@@ -201,5 +202,5 @@
 
     // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
       return false;
     }
@@ -222,13 +223,10 @@
 {
     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.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
+
     float AspectRatio = axes.minor / axes.major;
 
@@ -262,5 +260,4 @@
 {
     psF64 z;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
@@ -273,10 +270,7 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    // this estimates the radius assuming f(z) is roughly exp(-z)
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
+
     psF64 sigma = axes.major;
 
@@ -347,10 +341,8 @@
     }
 
-    // the OLD 2D model for SXY actually fitted SXY / (SXX^-2 + SYY^-2); correct here
-    // out[PM_PAR_SXY] = pmPSF_SXYtoModel (out);
-
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -403,6 +395,6 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;
Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_PS1_V1.c	(revision 35577)
+++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c	(revision 35768)
@@ -1,9 +1,10 @@
 /******************************************************************************
- * this file defines the PS1_V1 source shape model.  Note that these model functions are loaded
- * by pmModelClass.c using 'include', and thus need no 'include' statements of their own.  The
- * models use a psVector to represent the set of parameters, with the sequence used to specify
- * the meaning of the parameter.  The meaning of the parameters may thus vary depending on the
- * specifics of the model.  All models which are used as a PSF representations share a few
- * parameters, for which # define names are listed in pmModel.h:
+ * this file defines the PS1_V1 source shape model.  Note that these model functions are
+ * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
+ * their own.  The models use a psVector to represent the set of parameters, with the
+ * sequence used to specify the meaning of the parameter.  The meaning of the parameters
+ * may thus vary depending on the specifics of the model.  All models which are used as a
+ * PSF representations share a few parameters, for which # define names are listed in
+ * pmModel.h:
 
    power-law with fitted linear term
@@ -148,6 +149,6 @@
     float q2 = NAN;
     if (nParam == PM_PAR_SXY) {
-        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
         float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
         q1 = (q1 < 0.0) ? 0.0 : q1;
@@ -220,5 +221,5 @@
 
     // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
       return false;
     }
@@ -244,13 +245,9 @@
 {
     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.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     float AspectRatio = axes.minor / axes.major;
 
@@ -284,5 +281,4 @@
 {
     psF64 z;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
@@ -292,9 +288,6 @@
     if (flux >= PAR[PM_PAR_I0]) return 1.0;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     psF64 sigma = axes.major;
 
@@ -363,5 +356,6 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -416,5 +410,6 @@
     // convert to shape terms (SXX,SYY,SXY)
     // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;
Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 35577)
+++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 35768)
@@ -1,9 +1,10 @@
 /******************************************************************************
- * this file defines the QGAUSS source shape model (XXX need a better name!).  Note that these
- * model functions are loaded by pmModelClass.c using 'include', and thus need no 'include'
- * statements of their own.  The models use a psVector to represent the set of parameters, with
- * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
- * may thus vary depending on the specifics of the model.  All models which are used a PSF
- * representations share a few parameters, for which # define names are listed in pmModel.h:
+ * this file defines the QGAUSS source shape model.  Note that these model functions are
+ * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
+ * their own.  The models use a psVector to represent the set of parameters, with the
+ * sequence used to specify the meaning of the parameter.  The meaning of the parameters
+ * may thus vary depending on the specifics of the model.  All models which are used as a
+ * PSF representations share a few parameters, for which # define names are listed in
+ * pmModel.h:
 
    power-law with fitted linear term
@@ -14,6 +15,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)
@@ -138,5 +139,4 @@
 # define AR_MAX 20.0
 # define AR_RATIO 0.99
-
 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
 {
@@ -149,6 +149,7 @@
     float q2 = NAN;
     if (nParam == PM_PAR_SXY) {
-        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
+	// NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
         float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
         q1 = (q1 < 0.0) ? 0.0 : q1;
@@ -203,5 +204,5 @@
           return true;
       }
-    default:
+      default:
         psAbort("invalid choice for limits");
     }
@@ -221,5 +222,5 @@
 
     // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
       return false;
     }
@@ -245,13 +246,9 @@
 {
     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.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     float AspectRatio = axes.minor / axes.major;
 
@@ -285,5 +282,4 @@
 {
     psF64 z;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
@@ -293,11 +289,6 @@
     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;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     psF64 sigma = axes.major;
 
@@ -307,5 +298,4 @@
         return ( sigma * sqrt (2.0 * z) );
     }
-
     psF64 limit = flux / PAR[PM_PAR_I0];
 
@@ -367,5 +357,6 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -424,6 +415,6 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;
Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_RGAUSS.c	(revision 35577)
+++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c	(revision 35768)
@@ -138,6 +138,7 @@
     float q2 = NAN;
     if (nParam == PM_PAR_SXY) {
-        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
+	// NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
         float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
         q1 = (q1 < 0.0) ? 0.0 : q1;
@@ -210,5 +211,5 @@
 
     // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
       return false;
     }
@@ -234,13 +235,9 @@
 {
     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.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     float AspectRatio = axes.minor / axes.major;
 
@@ -274,5 +271,4 @@
 {
     psF64 z;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
@@ -285,9 +281,6 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     psF64 sigma = axes.major;
 
@@ -357,5 +350,6 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -409,6 +403,6 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;
Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_SERSIC.c	(revision 35577)
+++ trunk/psModules/src/objects/models/pmModel_SERSIC.c	(revision 35768)
@@ -125,12 +125,6 @@
 
 	// 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);
+	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
@@ -250,5 +244,6 @@
         // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
         // angle and let f2,f1 fight it out
-        q2 = 0.5*sqrtf(q1);
+	// NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        q2 = 2.0*0.5*sqrtf(q1);
     }
 
@@ -347,9 +342,6 @@
     axes.major = Rmajor;
     axes.minor = Rminor;
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
-    if (!isfinite(shape.sx))  return false;
-    if (!isfinite(shape.sy))  return false;
-    if (!isfinite(shape.sxy)) return false;
+
+    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], axes, true);
 
     float bn = 1.9992*index - 0.3271;
@@ -357,11 +349,4 @@
     float Io = exp(0.5*bn);
 
-    float Sxx = PS_MAX(0.5, shape.sx);
-    float Syy = PS_MAX(0.5, shape.sy);
-
-    PAR[PM_PAR_SXX]  = Sxx;
-    PAR[PM_PAR_SYY]  = Syy;
-    PAR[PM_PAR_SXY]  = shape.sxy;
-
     // set the model normalization (adjust for Sersic best guess)
     if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
@@ -381,14 +366,8 @@
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
-    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);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     float AspectRatio = axes.minor / axes.major;
 
@@ -410,6 +389,4 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
@@ -421,9 +398,6 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX];
-    shape.sy  = PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    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)
@@ -457,5 +431,6 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -510,5 +485,6 @@
     // convert to shape terms (SXX,SYY,SXY)
     // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;
Index: trunk/psModules/src/objects/models/pmModel_TRAIL.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_TRAIL.c	(revision 35577)
+++ trunk/psModules/src/objects/models/pmModel_TRAIL.c	(revision 35768)
@@ -350,10 +350,9 @@
     PAR[PM_PAR_SKY]  = 0.0;
 
-    // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
-    psEllipseShape psfShape;
-    psfShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
-    psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
-    psfShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
-    psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0);
+    psF32 *psfPAR  = source->modelPSF->params->data.F32;
+    bool useReff = pmModelUseReff (source->modelPSF->type);
+
+    psEllipseAxes psfAxes;
+    pmModelParamsToAxes (&psfAxes, psfPAR[PM_PAR_SXX], psfPAR[PM_PAR_SXY], psfPAR[PM_PAR_SYY], useReff);
 
     psEllipseMoments emoments;
@@ -369,7 +368,10 @@
     if (!isfinite(axes.theta)) return false;
 
-    float size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf);
-    //    if (size > psfAxes.major) { size -= psfAxes.major; }
-    //else { size = psfAxes.major; }
+    float size = NAN;
+    if (!isfinite(source->moments->Mrf)) {
+      size = axes.major;
+    } else {
+      size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf);
+    }
 
     float theta, peak;