Changeset 27838 for branches/tap_branches/psModules/src/objects/models/pmModel_PS1_V1.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_PS1_V1.c (modified) (7 diffs)

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

@@ -49 +49 @@
 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
@@ -262 +267 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -268 +273 @@
     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;
@@ -300 +297 @@
     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;
+        z = PS_MAX(z, 0.0);
     }
     psF64 radius = sigma * sqrt (2.0 * z);
@@ -345 +345 @@
     // 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;
@@ -454 +453 @@
         limitsApply = true;
         break;
+      case PM_MODEL_LIMITS_MODERATE:
+        paramsMinUse = paramsMinModerate;
+        paramsMaxUse = paramsMaxModerate;
+        limitsApply = true;
+        break;
       case PM_MODEL_LIMITS_STRICT:
         paramsMinUse = paramsMinStrict;
@@ -464 +468 @@
     return;
 }
-
 
 # undef PM_MODEL_FUNC