Changeset 20963 for trunk/psModules/src/objects/pmSourceMatch.c

Timestamp:

Dec 11, 2008, 3:54:17 PM (17 years ago)

Author:

Paul Price

Message:

Soften errors during iteration. Without this, it gets the wrong solution.

File:

• trunk/psModules/src/objects/pmSourceMatch.c (modified) (10 diffs)

trunk/psModules/src/objects/pmSourceMatch.c

@@ -290 +290 @@
                                 const psArray *matches, // Array of matches
                                 const psVector *zp, // Zero points for each image (including airmass term)
-                                const psVector *photo // Photometric image?
+                                const psVector *photo, // Photometric image?
+                                float sysErr2 // Systematic error, squared
                                 )
 {
@@ -321 +322 @@
             int index = match->image->data.U32[j]; // Image index
             float mag = match->mag->data.F32[j]; // Measured magnitude
-            double magErr = match->magErr->data.F32[j]; // Error in measured magnitude
-            double invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error
+            double magErr2 = PS_SQR(match->magErr->data.F32[j]) + sysErr2; // Error in measured magnitude
+            double invErr2 = 1.0 / magErr2; // Inverse square error
             float cal = zp->data.F32[index]; // Calibration to apply to image
             if (!photo || !photo->data.U8[index]) {
@@ -351 +352 @@
             int index = match->image->data.U32[j]; // Image index
             float mag = match->mag->data.F32[j]; // Measured magnitude
-            double magErr = match->magErr->data.F32[j]; // Error in measured magnitude
-            double invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error
+            double magErr2 = PS_SQR(match->magErr->data.F32[j]) + sysErr2; // Error in measured magnitude
+            double invErr2 = 1.0 / magErr2; // Inverse square error
             float cal = zp->data.F32[index]; // Calibration to apply to image
             accum->data.F64[index] += (mag + cal - stars->data.F32[i]) * invErr2;
@@ -377 +378 @@
             int index = match->image->data.U32[j]; // Image index
             float mag = match->mag->data.F32[j]; // Measured magnitude
-            float magErr = match->magErr->data.F32[j]; // Error in measured magnitude
+            float magErr2 = PS_SQR(match->magErr->data.F32[j]) + sysErr2; // Error in measured magnitude
             float cal = zp->data.F32[index]; // Calibration to apply to image
             if (!photo || !photo->data.U8[index]) {
                 cal -= trans->data.F32[index];
             }
-            float dev = (mag + cal - stars->data.F32[index]) / magErr; // Deviation from fit
-            if (isfinite(dev)) {
-                chi2 += PS_SQR(dev);
+            float dev2 = mag + cal - stars->data.F32[index]; // Deviation from fit
+            if (isfinite(dev2)) {
+                chi2 += PS_SQR(dev2) / magErr2;
                 dof++;
             }
@@ -449 +450 @@
                                       const psVector *photo, // Photometric image?
                                       float starClip, // Clipping for stars
-                                      float sysErr // Systematic error
+                                      float sysErr2 // Systematic error squared
                                 )
 {
@@ -467 +468 @@
 
     starClip = PS_SQR(starClip);
-    sysErr = PS_SQR(sysErr);
 
     int numRejected = 0;                // Number rejected
@@ -486 +486 @@
             }
             float dev = mag + cal - stars->data.F32[i]; // Deviation
-            if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr)) {
+            if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr2)) {
                 numRejected++;
                 match->mask->data.PS_TYPE_MASK_DATA[j] = 0xFF;
@@ -514 +514 @@
     PS_ASSERT_FLOAT_LARGER_THAN(transClip, 0.0, NULL);
 
+    sysErr *= sysErr;
+
     int numImages = zp->n;              // Number of images
     int numStars = matches->n;          // Number of stars
@@ -522 +524 @@
     psVector *stars = psVectorAlloc(numStars, PS_TYPE_F32); // Magnitudes for each star
 
-    float chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo); // chi^2 for solution
+    float chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo, sysErr); // chi^2 for solution
     psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2);
     float lastChi2 = INFINITY;          // chi^2 on last iteration
@@ -544 +546 @@
         psTrace("psModules.objects", 3, "%f%% of measurements rejected", fracRej * 100);
 
-        chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo); // chi^2 for solution
+        chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo, sysErr); // chi^2 for solution
         psTrace("psModules.objects", 1, "iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej);
     }