Changeset 30811

Timestamp:

Mar 6, 2011, 12:39:13 PM (15 years ago)

Author:

eugene

Message:

update output headers to include more info on the input images (airmasses, exptimes, zpts)

Location:

branches/eam_branches/ipp-20110213/ppStack/src

Files:

• ppStackConvolve.c (modified) (1 diff)
• ppStackOptions.h (modified) (1 diff)
• ppStackSources.c (modified) (9 diffs)
• ppStackUpdateHeader.c (modified) (2 diffs)

branches/eam_branches/ipp-20110213/ppStack/src/ppStackConvolve.c

@@ -323 +323 @@
 
         // UPDATE CELL.SATURATION here
-        UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",  options->sumExposure);
-        UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", options->sumExposure);
-        UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN);
+        UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",    options->sumExposure);
+        UPDATE_CONCEPT(outCell, "CELL.EXPOSURE",   options->sumExposure);
+        UPDATE_CONCEPT(outCell, "CELL.DARKTIME",   NAN);
         UPDATE_CONCEPT(outCell, "CELL.SATURATION", satBest);
-        UPDATE_CONCEPT(outFPA,  "FPA.ZP",        options->zp);
-
-        UPDATE_CONCEPT(unconvFPA,  "FPA.EXPOSURE",  options->sumExposure);
-        UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE", options->sumExposure);
-        UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME", NAN);
+        UPDATE_CONCEPT(outFPA,  "FPA.ZP",          options->zp);
+        UPDATE_CONCEPT(outFPA,  "FPA.AIRMASS",     options->airmass);
+
+        UPDATE_CONCEPT(unconvFPA,  "FPA.EXPOSURE",    options->sumExposure);
+        UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE",   options->sumExposure);
+        UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME",   NAN);
         UPDATE_CONCEPT(unconvCell, "CELL.SATURATION", satBest);
-        UPDATE_CONCEPT(unconvFPA,  "FPA.ZP",        options->zp);
+        UPDATE_CONCEPT(unconvFPA,  "FPA.ZP",          options->zp);
+        UPDATE_CONCEPT(unconvFPA,  "FPA.AIRMASS",     options->airmass);
 
         psLogMsg("ppStack", PS_LOG_INFO, "stack adjust metadata values : zp %f, exptime %f, sat %f", options->zp, options->sumExposure, satBest);

branches/eam_branches/ipp-20110213/ppStack/src/ppStackOptions.h

@@ -22 +22 @@
     float sumExposure;                  // Sum of exposure times
     float zp;                           // Zero point for output
+    float zpErr;                        // Zero point error for output
+    float airmass;                      // airmass for output image
+    float airmassSlope;                 // airmass slope used for zp analysis
     psVector *inputMask;                // Mask for inputs
     psArray *sourceLists;               // Individual lists of sources for matching
     psVector *norm;                     // Normalisation for each image
+    psVector *zpInput;                  // reported zero point of input exposure
+    psVector *expTimeInput;             // exposure time of input exposure
+    psVector *airmassInput;             // airmass of input exposure
     psArray *sources;                   // Matched sources
     // Convolve

branches/eam_branches/ipp-20110213/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)) {
@@ -277 +293 @@
 
         // 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
+        // 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
@@ -288 +306 @@
         // 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
+        // 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
@@ -305 +323 @@
         // 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
+        // m_inst_o - m_inst_i = zp[i] - zpTarget - c1 * airmassTarget - 2.5log(sumExpTime) - trans_i
 
         if (options->matchZPs) {
@@ -314 +332 @@
                 }
                 psArray *sources = sourceLists->data[i]; // Sources of interest
-                float magCorr = zp->data.F32[i] - trans->data.F32[i] - 2.5*log10(options->sumExposure);
+                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
@@ -341 +359 @@
             // Producing image with target zero point
             options->zp = zpTarget;
+            options->airmass = airmassTarget;
+            options->airmassSlope = airmassTerm;
         } else {
             options->zp = NAN;

branches/eam_branches/ipp-20110213/ppStack/src/ppStackUpdateHeader.c

@@ -70 +70 @@
     psMetadataAddStr(hdu->header, PS_LIST_TAIL, "TESS_ID", PS_META_REPLACE, "type of stack", tessID);
 
-    // add data about the inputs
+    psMetadataAddF32(hdu->header, PS_LIST_TAIL, "AIRM_SLP", PS_META_REPLACE, "airmass slope", options->airmassSlope);
+    
+
+    // write the following fields to the output headers:
+    // INP_NNNN : input file names
+    // SCL_NNNN : scale factor for each input file
+    // ZPT_NNNN : zero point of input file 
+    // EXT_NNNN : exptime of input file
+    // AIR_NNNN : airmass of input file
+
+    // use separate loops so they are blocked together in the headers (more readable)
+
+    // save the input filenames
     for (int i = 0; i < options->num; i++) {
         char field[64];
@@ -78 +90 @@
         psMetadataAddStr(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image name", basename);
         psFree(basename);
-
-        float inScale = pow(10.0, -0.4*options->norm->data.F32[i]);
-        snprintf (field, 64, "SCL_%04d", i);
-        psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image scale factor", inScale);
     }   
 
+    // save the input normalizations
+    for (int i = 0; i < options->num; i++) {
+        char field[64];
+
+        float value = options->norm ? pow(10.0, -0.4*options->norm->data.F32[i]) : NAN;
+        snprintf (field, 64, "SCL_%04d", i);
+        psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image scale factor", value);
+    }   
+    // save the input exptimes
+    for (int i = 0; i < options->num; i++) {
+        char field[64];
+
+        float value = options->zpInput ? options->zpInput->data.F32[i] : NAN;
+        snprintf (field, 64, "ZPT_%04d", i);
+        psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image zero point", value);
+    }   
+    // save the input normalizations
+    for (int i = 0; i < options->num; i++) {
+        char field[64];
+
+        float value = options->expTimeInput ? options->expTimeInput->data.F32[i] : NAN;
+        snprintf (field, 64, "EXP_%04d", i);
+        psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image exptime", value);
+    }   
+    // save the input normalizations
+    for (int i = 0; i < options->num; i++) {
+        char field[64];
+
+        float value = options->airmassInput ? options->airmassInput->data.F32[i] : NAN;
+        snprintf (field, 64, "AIR_%04d", i);
+        psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image airmass", value);
+    }   
     return true;
 }