Changeset 27838 for branches/tap_branches/psModules/src/objects/models/pmModel_QGAUSS.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/objects/models/pmModel_QGAUSS.c (modified) (10 diffs)

branches/tap_branches/psModules/src/objects/models/pmModel_QGAUSS.c

@@ -39 +39 @@
 # 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 +47 @@
 
 // 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
+// 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?
@@ -77 +86 @@
     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);
 
@@ -88 +97 @@
         // 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;
@@ -242 +251 @@
     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;
@@ -259 +268 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -265 +274 @@
     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;
@@ -286 +293 @@
 
     // 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;
+        z = PS_MAX(z, 0.0);
     }
     psF64 radius = sigma * sqrt (2.0 * z);
@@ -446 +454 @@
         limitsApply = true;
         break;
+      case PM_MODEL_LIMITS_MODERATE:
+        paramsMinUse = paramsMinModerate;
+        paramsMaxUse = paramsMaxModerate;
+        limitsApply = true;
+        break;
       case PM_MODEL_LIMITS_STRICT:
         paramsMinUse = paramsMinStrict;
@@ -466 +479 @@
 # undef PM_MODEL_FIT_STATUS
 # undef PM_MODEL_SET_LIMITS
+# undef ALPHA
+# undef ALPHA_M