Changeset 27531 for trunk/psModules/src/objects/models/pmModel_PS1_V1.c

Timestamp:

Mar 30, 2010, 1:29:50 PM (16 years ago)

Author:

eugene

Message:

add diff stats; add flux to CMF_PS1_DV1; output PSF residual mask; more accurate calculation of radius of source model

File:

• trunk/psModules/src/objects/models/pmModel_PS1_V1.c (modified) (5 diffs)

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

@@ -267 +267 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -273 +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;
@@ -305 +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;
     }
     psF64 radius = sigma * sqrt (2.0 * z);
@@ -350 +344 @@
     // 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;
@@ -475 +468 @@
 }
 
-
 # undef PM_MODEL_FUNC
 # undef PM_MODEL_FLUX