Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 6907)
+++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 6908)
@@ -26,8 +26,10 @@
     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,5 +41,4 @@
         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,4 +94,6 @@
 }
 
+// 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,5 +117,6 @@
 psF64 pmModelFlux_QGAUSS(const psVector *params)
 {
-    float f, norm, z;
+    float z;
+    float norm;
 
     psF32 *PAR = params->data.F32;
@@ -125,11 +129,31 @@
 
     // 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,4 +177,5 @@
 
     // 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,4 +189,5 @@
             break;
     }
+
     psF64 radius = sigma * sqrt (2.0 * z);
     if (isnan(radius)) {