Changeset 6947

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

Location:

trunk/psModules/src/objects/models

Files:

• pmModel_GAUSS.c (modified) (1 diff)
• pmModel_PGAUSS.c (modified) (3 diffs)
• pmModel_QGAUSS.c (modified) (7 diffs)

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

@@ -84 +84 @@
     params[2] = moments->x;
     params[3] = moments->y;
-    params[4] = 1.2 / moments->Sx;
-    params[5] = 1.2 / moments->Sy;
+    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);

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);
 }
 

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

@@ -95 +95 @@
 
 // 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)
 {
@@ -107 +106 @@
     params[2] = peak->x;
     params[3] = peak->y;
-    params[4] = 1.2 / moments->Sx;
-    params[5] = 1.2 / moments->Sy;
+    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;
     params[7] = 1.0;
@@ -129 +128 @@
 
     // 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;
 
-    # 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 *= DZ;
-    # else
-
-        float f0 = 1.0;
+    # define DZ 0.25
+
+    float f0 = 1.0;
     float f1, f2;
     for (z = DZ; z < 50; z += DZ) {
@@ -155 +142 @@
     }
     norm *= DZ / 3.0;
-    # endif
 
     psF64 Flux = PAR[1] * Area * norm;
@@ -168 +154 @@
     psF64 z, f;
     psF32 *PAR = params->data.F32;
+    int Nstep = 0;
 
     if (flux <= 0)
@@ -177 +164 @@
 
     // 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 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));
@@ -189 +179 @@
             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 */
+        z = 0;
+    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);