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

Timestamp:

Apr 18, 2006, 8:31:55 PM (20 years ago)

Author:

magnier

Message:

fixed handling of masks, fixed negative-radius errors in models

File:

• branches/rel10_ifa/psModules/src/objects/models/pmModel_QGAUSS.c (modified) (3 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);
+    }
+
     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;
@@ -111 +131 @@
 psF64 pmModelFlux_QGAUSS(const psVector *params)
 {
-    float f, norm, z;
+    float z;
+    float norm;
 
     psF32 *PAR = params->data.F32;
@@ -122 +143 @@
 
     // the area needs to be multiplied by the integral of f(z)
+    // XXX this integral can be done faster and more accurately
     norm = 0.0;
-    for (z = 0.05; z < 50; z += 0.1) {
+
+    # define DZ 0.1
+
+    # if 0
+
+    float f;
+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));
         norm += f;
     }
-    norm *= 0.1;
+    norm *= DZ;
+    # else
+
+        float f0 = 1.0;
+    float f1, f2;
+    for (z = DZ; z < 50; z += DZ) {
+        f1 = 1.0 / (1 + PAR[7]*z + pow(z, 2.25));
+        z += DZ;
+        f2 = 1.0 / (1 + PAR[7]*z + pow(z, 2.25));
+        norm += f0 + 4*f1 + f2;
+        f0 = f2;
+    }
+    norm *= DZ / 3.0;
+    # endif
 
     psF64 Flux = PAR[1] * Area * norm;