Changeset 6947 for trunk/psModules/src/objects/models/pmModel_PGAUSS.c

Timestamp:

Apr 21, 2006, 11:33:00 AM (20 years ago)

Author:

magnier

Message:

faster calculation of radius, simpsons rule for flux integration, rational range for Sx,Sy guess

File:

• trunk/psModules/src/objects/models/pmModel_PGAUSS.c (modified) (3 diffs)

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

@@ -76 +76 @@
 }
 
+// make an initial guess for parameters
+bool pmModelGuess_PGAUSS (pmModel *model, pmSource *source)
+{
+
+    pmMoments *moments = source->moments;
+    psF32     *params  = model->params->data.F32;
+
+    params[0] = moments->Sky;
+    params[1] = moments->Peak - moments->Sky;
+    params[2] = moments->x;
+    params[3] = moments->y;
+    params[4] = (moments->Sx < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sx);
+    params[5] = (moments->Sy < (1.2 / 2.0)) ? 2.0 : (1.2 / moments->Sy);
+    params[6] = 0.0;
+
+    return(true);
+}
+
 psF64 pmModelFlux_PGAUSS(const psVector *params)
 {
-    float f, norm, z;
+    float norm, z;
 
     psF32 *PAR = params->data.F32;
@@ -90 +108 @@
     // the area needs to be multiplied by the integral of f(z)
     norm = 0.0;
-    for (z = 0.005; z < 50; z += 0.01) {
-        f = 1.0 / (1 + z + z*z/2 + z*z*z/6);
-        norm += f;
+
+    # define DZ 0.25
+
+    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 *= 0.01;
+    norm *= DZ / 3.0;
 
-    psF64 Flux = params->data.F32[1] * Area * norm;
+    psF64 Flux = PAR[1] * Area * norm;
 
     return(Flux);
@@ -118 +144 @@
     }
     return (radius);
-}
-
-bool pmModelGuess_PGAUSS (pmModel *model, pmSource *source)
-{
-
-    pmMoments *moments = source->moments;
-    psF32     *params  = model->params->data.F32;
-
-    params[0] = moments->Sky;
-    params[1] = moments->Peak - moments->Sky;
-    params[2] = moments->x;
-    params[3] = moments->y;
-    params[4] = 1.2 / moments->Sx;
-    params[5] = 1.2 / moments->Sy;
-    params[6] = 0.0;
-
-    return(true);
 }