Changeset 14220 for trunk/psModules/src/objects/models/pmModel_PGAUSS.c

Timestamp:

Jul 15, 2007, 1:30:42 PM (19 years ago)

Author:

magnier

Message:

further attempt to set valid limits for sigma_xy

File:

• trunk/psModules/src/objects/models/pmModel_PGAUSS.c (modified) (6 diffs)

trunk/psModules/src/objects/models/pmModel_PGAUSS.c

@@ -39 +39 @@
     psF32 py = Y / PAR[PM_PAR_SYY];
     psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
+    assert (z >= 0.0);
+
     psF32 t  = 1 + z + z*z/2.0;
     psF32 r  = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */
@@ -63 +65 @@
 # define AR_MAX 20.0
 # define AR_RATIO 0.99
+
+// we constraint limits by the min valid minor axis:
+# define MIN_MINOR_AXIS 0.5
+
+// f3 = (s_b^-2 - s_a^-2); F3_SQ_MAX is MIN_MINOR_AXIS^-4
+# define F3_SQ_MAX 16.0 
+
+static float saveParams[8];
+
 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
 {
@@ -73 +84 @@
         f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
         q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
-        assert (q1 > 0);
+        q1 = PS_MAX (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*sqrt (q1);
     }
@@ -93 +106 @@
             break;
         case PM_PAR_SXX:
-            beta_lim = 0.5;
+            beta_lim = 2.0;
             break;
         case PM_PAR_SYY:
-            beta_lim = 0.5;
+            beta_lim = 2.0;
             break;
         case PM_PAR_SXY:
-            beta_lim =  q2;
+            beta_lim =  0.5*q2;
             break;
         default:
@@ -137 +150 @@
             psAbort("invalid parameter %d for param min test", nParam);
         }
-        if (params[nParam] < params_min) {
+        saveParams[nParam] = params[nParam];
+        if (params[nParam] < params_min) {
             params[nParam] = params_min;
             psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
@@ -170 +184 @@
             psAbort("invalid parameter %d for param max test", nParam);
         }
+        saveParams[nParam] = params[nParam];
         if (params[nParam] > params_max) {
             params[nParam] = params_max;