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

Timestamp:

Dec 9, 2008, 3:37:57 PM (17 years ago)

Author:

Paul Price

Message:

First cut at relphot. Still needs work, but basics look OK.

File:

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

trunk/psModules/src/objects/pmSourceMatch.c

@@ -176 +176 @@
         if (!boundsMaster) {
             // First run through --- can just copy
-            boundsMaster = boundsImage;
-            xMaster = xImage;
-            yMaster = yImage;
+            boundsMaster = psMemIncrRefCounter(boundsImage);
+            xMaster = psMemIncrRefCounter(xImage);
+            yMaster = psMemIncrRefCounter(yImage);
             matches = psArrayAlloc(numSources);
             for (int j = 0; j < numSources; j++) {
@@ -266 +266 @@
             if (i != numGood) {
                 psFree(matches->data[numGood]);
-                matches->data[numGood] = match;
-                psFree(match);
+                matches->data[numGood] = psMemIncrRefCounter(match);
+                //                psFree(match);
             }
             numGood++;
@@ -273 +273 @@
     }
     matches->n = numGood;
+    for (int i = numGood; i < numMaster; i++) {
+        psFree(matches->data[i]);
+        matches->data[i] = NULL;
+    }
 
     return matches;
@@ -281 +285 @@
 
 
-psVector *pmSourceRelphot(const psArray *matches, // Array of matches
-                          const psVector *zp, // Zero points for each image (including airmass term)
-                          float cloudClip // Clipping for clouds
-    )
+psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches
+                               const psVector *zp, // Zero points for each image (including airmass term)
+                               float cloudClip // Clipping for clouds
+                               )
 {
     PS_ASSERT_ARRAY_NON_NULL(matches, NULL);
@@ -291 +295 @@
     PS_ASSERT_FLOAT_LARGER_THAN(cloudClip, 0.0, NULL);
 
-    return NULL;
-}
+    int numImages = zp->n;              // Number of images
+    int numStars = matches->n;          // Number of stars
+    psVector *trans = psVectorAlloc(numImages, PS_TYPE_F32); // Transparencies for each image, magnitudes
+    psVector *transErr = psVectorAlloc(numImages, PS_TYPE_F32); // Errors in the transparency
+    psVector *stars = psVectorAlloc(numStars, PS_TYPE_F32); // Magnitudes for each star
+    psVector *starsErr = psVectorAlloc(numStars, PS_TYPE_F32); // Errors in the star magnitudes
+
+
+    // XXX Need to apply provided ZPs
+
+    psVectorInit(trans, 0.0);
+    psVectorInit(transErr, 0.0);
+
+    // Solve the star magnitudes
+    psVectorInit(stars, 0.0);
+    psVectorInit(starsErr, 0.0);
+    for (int i = 0; i < numStars; i++) {
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        for (int j = 0; j < match->num; j++) {
+            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 invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error
+            stars->data.F32[i] += (mag - trans->data.F32[index]) * invErr2;
+            starsErr->data.F32[i] += invErr2;
+        }
+    }
+    for (int i = 0; i < numStars; i++) {
+        float inverse = 1.0 / starsErr->data.F32[i]; // Inverse error
+        stars->data.F32[i] *= inverse;
+        starsErr->data.F32[i] = sqrtf(inverse);
+    }
+
+    // Solve for the transparencies
+    psVectorInit(trans, 0.0);
+    psVectorInit(transErr, 0.0);
+    for (int i = 0; i < numStars; i++) {
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        for (int j = 0; j < match->num; j++) {
+            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 invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error
+            trans->data.F32[index] += (mag - stars->data.F32[i]) * invErr2;
+            transErr->data.F32[index] += invErr2;
+        }
+    }
+    for (int i = 0; i < numImages; i++) {
+        float inverse = 1.0 / transErr->data.F32[i]; // Inverse error
+        trans->data.F32[i] *= inverse;
+        transErr->data.F32[i] = sqrtf(inverse);
+
+        psTrace("psModules.objects", 1, "Transparency for image %d: %f +/- %f\n",
+                i, trans->data.F32[i], transErr->data.F32[i]);
+    }
+
+    psFree(transErr);
+    psFree(stars);
+    psFree(starsErr);
+
+    return trans;
+}