Changeset 27531 for trunk/psModules/src/objects/models/pmModel_PGAUSS.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_PGAUSS.c (modified) (2 diffs)

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

@@ -236 +236 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -271 +271 @@
     z0 = 0.0;
 
-    // perform a type of bisection to find the value
-    float f0 = 1.0 / (1 + z0 + z0*z0/2.0 + z0*z0*z0/6.0);
-    float f1 = 1.0 / (1 + z1 + z1*z1/2.0 + z1*z1*z1/6.0);
-    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
-        z = 0.5*(z0 + z1);
-        f = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
-        if (f > limit) {
-            z0 = z;
-            f0 = f;
-        } else {
-            z1 = z;
-            f1 = f;
-        }
-        Nstep ++;
-    }
+    // starting guess:
+    z = 0.5*(z0 + z1);
+    float dz = 1.0;
+
+    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
+        // use Newton-Raphson to minimize f(z) - limit = 0
+        float dqdz = (1.0 + z + z*z/2.0);
+        float q = (dqdz + z*z*z/6.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);