Index: /branches/eam_branches/ipp-20120627/psModules/src/objects/models/pmModel_TRAIL.c
===================================================================
--- /branches/eam_branches/ipp-20120627/psModules/src/objects/models/pmModel_TRAIL.c	(revision 34140)
+++ /branches/eam_branches/ipp-20120627/psModules/src/objects/models/pmModel_TRAIL.c	(revision 34141)
@@ -59,6 +59,6 @@
 
 // Lax parameter limits 
-static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0 };
-static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 };
+static float paramsMinLax[] = { -1.0e3, 1.0e-2, -1.0e2, -1.0e2,  0.5, -0.1, -1.0 };
+static float paramsMaxLax[] = {  1.0e5, 1.0e+8, +1.0e4, +1.0e4, 50.0, +3.3 , 1.0 };
 
 // Moderate parameter limits
@@ -73,5 +73,5 @@
 static float *paramsMinUse = paramsMinLax;
 static float *paramsMaxUse = paramsMaxLax;
-static float betaUse[] = { 1000, 3e6, 5, 5, 2.0, 2.0, 0.5 };
+static float betaUse[] = { 1000, 3e6, 5, 5, 2.0, 0.1, 0.1 };
 
 static bool limitsApply = true;         // Apply limits?
@@ -166,5 +166,5 @@
         dPAR[PM_PAR_LENGTH] = PAR[PM_PAR_I0] * dPdL;
         dPAR[PM_PAR_THETA]  = PAR[PM_PAR_I0] * dPdT;
-        dPAR[PM_PAR_SIGMA]  = NAN;	// we don't actually allow this to vary, so we do not need to calculate it
+        dPAR[PM_PAR_SIGMA]  = 0;	// we don't actually allow this to vary, so we do not need to calculate it
     }
     return(f);
@@ -184,23 +184,8 @@
     psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
 
-    // we need to calculate the limits for SXY specially
-    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 q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
-        q1 = (q1 < 0.0) ? 0.0 : q1;
-        // 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);
-    }
-
     switch (mode) {
       case PS_MINIMIZE_BETA_LIMIT: {
           psAssert(beta, "Require beta to limit beta");
           float limit = betaUse[nParam];
-          if (nParam == PM_PAR_SXY) {
-              limit *= q2;
-          }
           if (fabs(beta[nParam]) > fabs(limit)) {
               beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
@@ -215,7 +200,4 @@
           psAssert(paramsMinUse, "Require parameter limits to limit parameters");
           float limit = paramsMinUse[nParam];
-          if (nParam == PM_PAR_SXY) {
-              limit *= q2;
-          }
           if (params[nParam] < limit) {
               params[nParam] = limit;
@@ -230,7 +212,4 @@
           psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
           float limit = paramsMaxUse[nParam];
-          if (nParam == PM_PAR_SXY) {
-              limit *= q2;
-          }
           if (params[nParam] > limit) {
               params[nParam] = limit;
@@ -257,8 +236,26 @@
     PAR[PM_PAR_SKY]  = 0.0;
 
-    // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
-      return false;
-    }
+    // 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);
+
+    psEllipseMoments emoments;
+    emoments.x2 = source->moments->Mxx;
+    emoments.xy = source->moments->Mxy;
+    emoments.y2 = source->moments->Myy;
+
+    // force the axis ratio to be < 20.0
+    psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
+
+    if (!isfinite(axes.major)) return false;
+    if (!isfinite(axes.minor)) return false;
+    if (!isfinite(axes.theta)) return false;
+
+    PAR[PM_PAR_LENGTH] = axes.major; // XXX major or 2*major?
+    PAR[PM_PAR_THETA] = axes.theta; // XXX theta or theta * RAD_DEG?
+    PAR[PM_PAR_SIGMA] = psfAxes.major; // XXX major or 0.5*major?
 
     // set the model normalization