Index: /branches/eam_branches/ipp-20110213/ppStack/src/ppStackConvolve.c =================================================================== --- /branches/eam_branches/ipp-20110213/ppStack/src/ppStackConvolve.c (revision 30810) +++ /branches/eam_branches/ipp-20110213/ppStack/src/ppStackConvolve.c (revision 30811) @@ -323,15 +323,17 @@ // 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); Index: /branches/eam_branches/ipp-20110213/ppStack/src/ppStackOptions.h =================================================================== --- /branches/eam_branches/ipp-20110213/ppStack/src/ppStackOptions.h (revision 30810) +++ /branches/eam_branches/ipp-20110213/ppStack/src/ppStackOptions.h (revision 30811) @@ -22,7 +22,13 @@ 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 Index: /branches/eam_branches/ipp-20110213/ppStack/src/ppStackSources.c =================================================================== --- /branches/eam_branches/ipp-20110213/ppStack/src/ppStackSources.c (revision 30810) +++ /branches/eam_branches/ipp-20110213/ppStack/src/ppStackSources.c (revision 30811) @@ -115,5 +115,11 @@ 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,4 +134,9 @@ 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,4 +178,10 @@ 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,5 +202,4 @@ if (!filter) { filter = expFilter; - bool mdok; airmassTerm = psMetadataLookupF32(&mdok, airmassZP, filter); if (!mdok || !isfinite(airmassTerm)) { @@ -277,5 +293,7 @@ // 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,7 +306,7 @@ // 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,5 +323,5 @@ // 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,5 +332,5 @@ } 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,4 +359,6 @@ // Producing image with target zero point options->zp = zpTarget; + options->airmass = airmassTarget; + options->airmassSlope = airmassTerm; } else { options->zp = NAN; Index: /branches/eam_branches/ipp-20110213/ppStack/src/ppStackUpdateHeader.c =================================================================== --- /branches/eam_branches/ipp-20110213/ppStack/src/ppStackUpdateHeader.c (revision 30810) +++ /branches/eam_branches/ipp-20110213/ppStack/src/ppStackUpdateHeader.c (revision 30811) @@ -70,5 +70,17 @@ 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,10 +90,38 @@ 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; }