Changeset 6908

Timestamp:

Apr 19, 2006, 10:33:06 AM (20 years ago)

Author:

magnier

Message:

added simpsons rule for flux integral

File:

• trunk/psModules/src/objects/models/pmModel_QGAUSS.c (modified) (7 diffs)

trunk/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)
+        z = 0;
+
     psF32 zp = pow(z,1.25);
-
     psF32 r  = 1.0 / (1 + PAR[7]*z + z*zp);
-    // test: psF32 r  = 1.0 / (1 + PAR[7]*z + PS_SQR(z));
+
     psF32 r1 = PAR[1]*r;
     psF32 f  = r1 + PAR[0];
@@ -39 +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;
@@ -93 +94 @@
 }
 
+// make an initial guess for parameters
+// XXX we could probably do better with params[6] and params[7]
 bool pmModelGuess_QGAUSS (pmModel *model, pmSource *source)
 {
@@ -114 +117 @@
 psF64 pmModelFlux_QGAUSS(const psVector *params)
 {
-    float f, norm, z;
+    float z;
+    float norm;
 
     psF32 *PAR = params->data.F32;
@@ -125 +129 @@
 
     // 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 1
+
+    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;
@@ -153 +177 @@
 
     // if Sx == Sy, sigma = Sx == Sy
+    // XXX we should return the major axis length...??
     psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0);
     psF64 dz = 1.0 / (2.0 * sigma*sigma);
@@ -164 +189 @@
             break;
     }
+
     psF64 radius = sigma * sqrt (2.0 * z);
     if (isnan(radius)) {