Changeset 31158 for trunk/ppStack/src/ppStackSources.c

Timestamp:

Apr 4, 2011, 1:16:41 PM (15 years ago)

Author:

eugene

Message:

improvements to the output headers; some code organization; plug leaks; user-specified limits to input seeing ranges

Location:

trunk/ppStack

Files:

• . (modified) (2 props)
• src/ppStackSources.c (modified) (9 diffs)

trunk/ppStack

@@ -17 +17 @@
 autom4te.cache
 Doxyfile
+a.out.dSYM

@@ -17 +17 @@
 autom4te.cache
 Doxyfile
+a.out.dSYM

trunk/ppStack/src/ppStackSources.c

@@ -115 +115 @@
     float fracMatch = psMetadataLookupF32(NULL, recipe, "ZP.MATCH"); // Fraction of images to match for star
 
-    psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms
+    bool mdok = false;
+    float airmassTarget = psMetadataLookupF32(&mdok, recipe, "ZP.AIRMASS.TARGET"); // output airmass value 
+    if (!mdok) {
+        airmassTarget = 1.0;
+    }
+
+    psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms (slopes) by filter
     if (!airmassZP) {
         psError(PPSTACK_ERR_CONFIG, false, "Unable to find ZP.AIRMASS in recipe.");
@@ -128 +134 @@
     int num = options->num;             // Number of inputs
     psAssert(num == sourceLists->n, "Wrong number of source lists: %ld\n", sourceLists->n);
+
+    // vectors to store these inputs so they may be recorded in the output headers
+    options->zpInput      = psVectorAlloc(num, PS_TYPE_F32);
+    options->expTimeInput = psVectorAlloc(num, PS_TYPE_F32);
+    options->airmassInput = psVectorAlloc(num, PS_TYPE_F32);
 
     psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Relative zero points for each image
@@ -167 +178 @@
         const char *expFilter = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.FILTER"); // Filter name
         zpExp->data.F32[i] = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.ZP"); // Exposure zero point
+        // XXX need to get the zero point error values and propagate to get the FPA.ZP.ERR value
+
+        options->zpInput->data.F32[i] = zpExp->data.F32[i]; // NOTE zpExp may be re-assigned below using relative photometry
+        options->expTimeInput->data.F32[i] = exptime;
+        options->airmassInput->data.F32[i] = airmass;
+
         psLogMsg("ppStack", PS_LOG_INFO,
                  "Image %d: %.2f sec exposure in %s at airmass %.2f with zero point %.2f",
@@ -185 +202 @@
         if (!filter) {
             filter = expFilter;
-            bool mdok;
             airmassTerm = psMetadataLookupF32(&mdok, airmassZP, filter);
             if (!mdok || !isfinite(airmassTerm)) {
@@ -206 +222 @@
         }
 
+        // XXX this is wrong, or at least inconsistent with the above: this needs to include 
+        // a value for the nominal system zero point to be consistent with zpExp
         zp->data.F32[i] = airmassTerm * airmass + 2.5 * log10(exptime);
     }
@@ -236 +254 @@
 
     if (zpExpNum == numGoodImages) {
+        psLogMsg ("ppStack", PS_LOG_INFO, "all zero points are finite; using reported zero points listed above");
         for (int i = 0; i < num; i++) {
             zp->data.F32[i] = zpExp->data.F32[i];
+        }
+    } else {
+        psLogMsg ("ppStack", PS_LOG_INFO, "missing some zero points; using guess values:");
+        for (int i = 0; i < num; i++) {
+            psLogMsg("ppStack", PS_LOG_INFO, "Image %d: %.2f sec exposure with zero point %.2f", i, options->exposures->data.F32[i], zp->data.F32[i]);
         }
     }
@@ -265 +289 @@
             }
         }
-        // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
-        // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
-        // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
-        // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
-        // We don't need to know the magnitude zero point for the filter, since it cancels out
+
+        // EAM : the discussion here was not quite right (or at least sloppy).  Here is a replacement explanation:
+
+        // For any star, the observed instrumental magnitude on an image and the apparent magnitude are related by:
+        // M_app = m_inst + zp + c1 * airmass + 2.5log(t) - transparency
+        // NOTE the sign of 'transparency'  this must agree with the definition in pmSourceMatch.c. see, eg, line 457 where 
+        // transparency = m_inst + zp + c1 * airmass + 2.5log(t) - M_app 
+
+        // we want to adjust the input images to be in a consistent flux system so that the
+        // final stack can be generated with a specific target zero point.  Any adjustment to
+        // the flux scale of the image must be made in coordination with the resulting
+        // zeropoint, exposure time, and airmass such that the above relationship yields the
+        // same apparent magnitude for a given star:
+
+        // m_inst_i : instrumental mags on input image (in)
+        // m_inst_o : instrumental mags on re-normalized image (out)
+
+        // m_inst_o + zp_o + c1 * airmass_o + 2.5log(t_o) - trans_o = m_inst_i + zp_i + c1 * airmass_i + 2.5log(t_i) - trans_i
+
+        // m_inst_o = m_inst_i + (zp_i - zp_o) + c1 * (airmass_i - airmass_o) + 2.5log(t_i) - 2.5log(t_o) - trans_i + trans_o
+
+        // zp_i, airmass_i, t_i, trans_i : reported or measured for input image
+
+        // zp_o      = zpTarget      (from recipe)
+        // airmass_o = airmassTarget (from recipe)
+        // t_o       = sumExpTime    [sum of input exposure times: once images are scale to this time, they can be avereaged]
+        // trans_o   = 0.0           [obviously!]
+
+        // we have 2 cases: (a) all reported ZPs are good or (b) some are bad:
+        // (a) FPA.ZP = zp_i + c1 * airmass_i
+        //  --> zp[i] = zp_i + c1 * airmass_i + 2.5log(exptime_i)
+        // (b)  zp[i] = c1 * airmass_i + 2.5log(exptime_i)
+        // NOTE: in case (b), the current code is equating the TARGET zp with the NOMINAL zp, which is wrong.
+
+        // m_inst_o - m_inst_i = zp[i] - zpTarget - c1 * airmassTarget - 2.5log(sumExpTime) - trans_i
+
         if (options->matchZPs) {
             options->norm = psVectorAlloc(num, PS_TYPE_F32);
@@ -277 +332 @@
                 }
                 psArray *sources = sourceLists->data[i]; // Sources of interest
-                float magCorr = zp->data.F32[i] - trans->data.F32[i] - 2.5*log10(options->sumExposure);
-                if (zpExpNum == numGoodImages) {
+                float magCorr = zp->data.F32[i] - trans->data.F32[i] - 2.5*log10(options->sumExposure) - airmassTerm * airmassTarget;
+                if (zpExpNum == numGoodImages) { // case (a)
                     // Using measured zero points, so attempt to set target zero point
+                    // XXX see NOTE above regarding case (b) : this is wrong.  the code should load a nominal zero point and supply it above
+                    // 
                     magCorr -= zpTarget;
                 }
@@ -302 +359 @@
             // Producing image with target zero point
             options->zp = zpTarget;
+            options->airmass = airmassTarget;
+            options->airmassSlope = airmassTerm;
         } else {
             options->zp = NAN;

@@ -17 +17 @@
 autom4te.cache
 Doxyfile
+a.out.dSYM