Index: trunk/psastro/src/psastroConvert.c =================================================================== --- trunk/psastro/src/psastroConvert.c (revision 31661) +++ trunk/psastro/src/psastroConvert.c (revision 36837) @@ -16,5 +16,5 @@ static psArray *chooseStars(psArray *inStars, char *listName, psArray *sources, psVector *index, int nMax, float iMagMin, float iMagMax, pmSourceMode skip); -bool psastroConvertFPA (pmFPA *fpa, psMetadata *recipe) { +bool psastroConvertFPA (pmConfig *config, pmFPA *fpa, psMetadata *recipe) { pmChip *chip; @@ -35,5 +35,8 @@ if (! readout->data_exists) { continue; } - psastroConvertReadout (readout, recipe); + if (!psastroConvertReadout (config, view, readout, recipe)) { + psError(PSASTRO_ERR_CONFIG, true, "problem converting readout\n"); + return false; + } } } @@ -44,5 +47,5 @@ // pass/apply the WCS information? -bool psastroConvertReadout (pmReadout *readout, psMetadata *recipe) { +bool psastroConvertReadout (pmConfig *config, pmFPAview *view, pmReadout *readout, psMetadata *recipe) { bool status; @@ -57,4 +60,10 @@ // convert the pmSource objects into pmAstromObj objects (drop !STAR and SATSTAR?) psArray *inStars = pmSourceToAstromObj (sources); + + // apply Koppenhoefer correction if needed + if (!psastroCorrectKH (config, view, readout, recipe, inStars)){ + psError(PSASTRO_ERR_CONFIG, true, "failed to correct Koppenhoefer Effect\n"); + return false; + } // sort in ascending magnitude order @@ -183,4 +192,13 @@ psF32 *dPAR = model->dparams->data.F32; + // Note from EAM: I measured the KH correction from the nightly science database + // The KH effect was active until 2011/05/11 when the camera voltages were modified + // Between 2011/09/06 and 2013/04/30, there was a bug in pmPSF_ModelToAxes with the + // result that the reported values were too large by a factor of sqrt(2). Fortunately, + // this did not affect the data in the nightly science DVO used to measure the + // effect (all affected data was processed BEFORE the bug was introduced), and it + // does not affect any of the PV2 or PV3 data processed AFTER the bug was fixed. + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type); + pmAstromObj *obj = pmAstromObjAlloc (); @@ -192,4 +210,5 @@ obj->Mag = source->psfMag; obj->dMag = source->psfMagErr; + obj->SBinst = source->psfMag + 5.0*log10(axes.major); // XXX do we have the information giving the readout and cell offset? @@ -221,2 +240,225 @@ } +static float pltScale[] = { ++0.0, +-0.255347177386, +-0.255585625172, +-0.255733260512, +-0.255747352242, +-0.255611001253, +-0.255311933756, + +0.0, + +0.0, + +0.0, +-0.255296201766, +-0.255688636720, +-0.256161080003, +-0.256464074850, +-0.256324860930, +-0.256147363186, +-0.255734743476, +-0.255422729909, + +0.0, + +0.0, +-0.255508323163, +-0.256142899513, +-0.255493895173, +-0.256997376740, +-0.256935591102, +-0.256644588679, +-0.256123321533, +-0.255515242517, + +0.0, + +0.0, +-0.255633033872, +-0.256356916189, +-0.256929253668, +-0.257250442505, +-0.257260403216, +-0.256962024927, +-0.256387320280, +-0.255681164920, + +0.0, + +0.0, +-0.255671992302, +-0.256371312678, +-0.256936249495, +-0.257237892330, +-0.257248004675, +-0.256911942035, +-0.256365648448, +-0.255632205397, + +0.0, + +0.0, +-0.255455072403, +-0.256080205411, +-0.256570428520, +-0.256908982217, +-0.256899240971, +-0.256607972383, +-0.256107793808, +-0.255482232273, + +0.0, + +0.0, +-0.255237901688, +-0.255669494152, +-0.256101093233, +-0.256372770309, +-0.256373448133, +-0.256118197262, +-0.255689391643, +-0.255279636681, + +0.0, + +0.0, + +0.0, +-0.255286350489, +-0.255517200172, +-0.255689459741, +-0.255668706447, +-0.255539349556, +-0.255363308907, +}; + +bool psastroCorrectKH (pmConfig *config, pmFPAview *view, pmReadout *readout, psMetadata *recipe, psArray *inStars) { + + bool status; + + psAssert (readout, "missing readout"); + psAssert (readout->parent, "missing cell"); + psAssert (readout->parent->parent, "missing chip"); + + // should we stay or should we go? + bool apply = psMetadataLookupBool (&status, recipe, "KH.CORRECT.APPLY.EXP"); + if (!apply) return true; + + // get the chip ID from the chip name: + pmChip *chip = readout->parent->parent; + char *chipName = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + + psAssert (strlen(chipName) == 4, "error in chip name"); + psAssert (chipName[0] == 'X', "error in chip name"); + psAssert (chipName[1] == 'Y', "error in chip name"); + + int chipID = atoi (&chipName[2]); + + // XXX hardwired list of chips to correct + // some notes: + // XY24 surpisingly does not need correction + // XY27 has poor astrometry with a weird correction for all mags + // XY36 surpisingly does not need correction + // XY67 has poor astrometry with a weird correction for all mags + + // correct this chip? (raise an error for unexpected chipID values) + switch (chipID) { + case 1: // do not correct + case 2: // do not correct + case 3: // do not correct + case 10: // do not correct + case 11: // do not correct + case 12: // do not correct + case 13: // do not correct + case 20: // do not correct + case 21: // do not correct + case 22: // do not correct + case 23: // do not correct + case 24: // do not correct + case 27: // do not correct + case 30: // do not correct + case 31: // do not correct + case 32: // do not correct + case 33: // do not correct + case 36: // do not correct + case 44: // do not correct + case 45: // do not correct + case 46: // do not correct + case 47: // do not correct + case 54: // do not correct + case 55: // do not correct + case 56: // do not correct + case 57: // do not correct + case 64: // do not correct + case 65: // do not correct + case 66: // do not correct + case 67: // do not correct + case 74: // do not correct + case 75: // do not correct + case 76: // do not correct + return true; + case 4: // correct + case 5: // correct + case 6: // correct + case 14: // correct + case 15: // correct + case 16: // correct + case 17: // correct + case 25: // correct + case 26: // correct + case 34: // correct + case 35: // correct + case 37: // correct + case 40: // correct + case 41: // correct + case 42: // correct + case 43: // correct + case 50: // correct + case 51: // correct + case 52: // correct + case 53: // correct + case 60: // correct + case 61: // correct + case 62: // correct + case 63: // correct + case 71: // correct + case 72: // correct + case 73: // correct + break; + default: + psAbort ("chipID is invalid"); + } + + // grab the KH correction file + pmFPAfile *KHfile = psMetadataLookupPtr (NULL, config->files, "PSASTRO.KH.CORRECT"); + if (!KHfile) { + psError(PM_ERR_CONFIG, false, "KH correction file not found"); + return false; + } + + // grab the corresponding chip + pmChip *KHchip = pmFPAviewThisChip (view, KHfile->fpa); + psAssert (KHchip, "found KH file, but not chip?"); + + // grab the correction spline + KHcorrectData *spline = psMetadataLookupPtr (&status, KHchip->analysis, "KH.CORRECT"); + if (!spline) { + psError(PM_ERR_CONFIG, false, "KH correction not found"); + return false; + } + + // the hard-wired array of plate-scales per chip above come from + // a single smf. in there, a negative plate scale means parity is flipped + // on the sky. + float myPltScale = fabs(pltScale[chipID]); + + // apply the correction to the detections + for (int i = 0; i < inStars->n; i++) { + pmAstromObj *obj = inStars->data[i]; + + float Xraw = obj->pix->x; + float SBinst = obj->SBinst; + + float dX = KHcorrectApply (spline, SBinst); + + float Xfix = Xraw + dX / myPltScale; + + // note that we carry around pix, cell, chip but the real analysis only operates on chip + // if in the future we add transformations between chip->cell->pix, then we will need to + // make these consistent as well. + obj->pix->x = Xfix; + obj->cell->x = Xfix; + obj->chip->x = Xfix; + } + + return true; + +} +