Index: trunk/ppStack/src/ppStackSources.c =================================================================== --- trunk/ppStack/src/ppStackSources.c (revision 30685) +++ trunk/ppStack/src/ppStackSources.c (revision 31158) @@ -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)) { @@ -206,4 +222,6 @@ } + // 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,6 +254,12 @@ 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,9 +289,40 @@ } } - // 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,7 +332,9 @@ } 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,4 +359,6 @@ // Producing image with target zero point options->zp = zpTarget; + options->airmass = airmassTarget; + options->airmassSlope = airmassTerm; } else { options->zp = NAN;