Changeset 6916 for branches/rel10_ifa/psModules/src/objects/models/pmModel_QGAUSS.c

Timestamp:

Apr 19, 2006, 4:38:34 PM (20 years ago)

Author:

magnier

Message:

fixed error in z limits (neg goes to 0); better radius search

File:

• branches/rel10_ifa/psModules/src/objects/models/pmModel_QGAUSS.c (modified) (6 diffs)

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

@@ -26 +26 @@
     psF32 py = PAR[5]*Y;
     psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y;
-    if (z <= 0) {
-        if (deriv) {
-            psF32 *dPAR = deriv->data.F32;
-            dPAR[0] = 0.0;
-            dPAR[1] = 0.0;
-            dPAR[2] = 0.0;
-            dPAR[3] = 0.0;
-            dPAR[4] = 0.0;
-            dPAR[5] = 0.0;
-            dPAR[6] = 0.0;
-            dPAR[7] = 0.0;
-        }
-        return(0.0);
-    }
+    if (z < 0)
+        z = 0;
 
     psF32 zp = pow(z,1.25);
-    if (isnan(zp))
-        psAbort ("psMinLMM", "nan in zp");
-
     psF32 r  = 1.0 / (1 + PAR[7]*z + z*zp);
-    if (isnan(r))
-        psAbort ("psMinLMM", "nan in r");
-
-    // test: psF32 r  = 1.0 / (1 + PAR[7]*z + PS_SQR(z));
+
     psF32 r1 = PAR[1]*r;
     psF32 f  = r1 + PAR[0];
@@ -59 +41 @@
         psF32 t = r1*r;
         psF32 q = t*(PAR[7] + 2.25*zp);
-        // test: psF32 q = t*(PAR[7] + 2*z);
 
         dPAR[0] = +1.0;
@@ -146 +127 @@
     norm = 0.0;
 
-    # define DZ 0.1
+    # define DZ 0.25
 
     # if 0
 
     float f;
+float zp;
 for (z = 0.5*DZ; z < 50; z += DZ) {
-        f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25));
-        // test: f = 1.0 / (1 + PAR[7]*z + PS_SQR(z));
+        zp = pow(z,2.25);
+        f = 1.0 / (1 + PAR[7]*z + zp);
         norm += f;
     }
@@ -182 +164 @@
     psF64 z, f;
     psF32 *PAR = params->data.F32;
+    int Nstep = 0;
 
     if (flux <= 0)
@@ -192 +175 @@
     // if Sx == Sy, sigma = Sx == Sy
     psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0);
+    psF64 limit = flux / PAR[1];
+
+    # if 0
+    /* test example will just use both, printing both */
     psF64 dz = 1.0 / (2.0 * sigma*sigma);
-    psF64 limit = flux / PAR[1];
 
     // we can do this much better with intelligent choices here
     for (z = 0.0; z < 30.0; z += dz) {
+        Nstep ++;
         f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25));
         // test: f = 1.0 / (1 + PAR[7]*z + PS_SQR(z));
@@ -202 +189 @@
             break;
     }
+    // fprintf (stderr, "rad 1: %f, want: %f, got: %f (%d steps)\n", z, limit, f, Nstep);
+
+    # else
+
+        /* use the fact that f is monotonically decreasing */
+        Nstep = 0;
+
+    // choose a z value guaranteed to be beyond our limit
+    float z0 = pow((1.0 / limit), (1.0 / 2.25));
+    float z1 = (1.0 / limit) / PAR[7];
+    z1 = PS_MAX (z0, z1);
+    z0 = 0.0;
+
+    float f0 = 1.0 / (1 + PAR[7]*z0 + pow(z0, 2.25));
+    float f1 = 1.0 / (1 + 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[7]*z + pow(z, 2.25));
+        // fprintf (stderr, "%f  %f  %f   :   %f  %f  %f\n", f0, f, f1, z0, z, z1);
+        if (f > limit) {
+            z0 = z;
+            f0 = f;
+        } else {
+            z1 = z;
+            f1 = f;
+        }
+        Nstep ++;
+    }
+    // fprintf (stderr, "rad 2: %f, want: %f, got: %f (%d steps)\n", z, limit, f, Nstep);
+    # endif
+
     psF64 radius = sigma * sqrt (2.0 * z);
     if (isnan(radius)) {