Index: /trunk/psastro/src/psastro.h =================================================================== --- /trunk/psastro/src/psastro.h (revision 36836) +++ /trunk/psastro/src/psastro.h (revision 36837) @@ -66,6 +66,8 @@ bool psastroAnalysis (pmConfig *config, psMetadata *stats); -bool psastroConvertFPA (pmFPA *fpa, psMetadata *recipe); -bool psastroConvertReadout (pmReadout *readout, psMetadata *recipe); +bool psastroConvertFPA (pmConfig *config, pmFPA *fpa, psMetadata *recipe); +bool psastroConvertReadout (pmConfig *config, pmFPAview *view, pmReadout *readout, psMetadata *recipe); +bool psastroCorrectKH (pmConfig *config, pmFPAview *view, pmReadout *readout, psMetadata *recipe, psArray *inStars); + psArray *pmSourceToAstromObj (psArray *sources); bool psastroAstromGuess (int *nStars, pmConfig *config); @@ -90,4 +92,6 @@ psArray *psastroRemoveClumpsIterate (psArray *input, int scale, int nIter); bool psastroRemoveClumpsRawstars (pmConfig *config); + +bool psastroChipFailureHeader (psMetadata *updates); Index: /trunk/psastro/src/psastroArguments.c =================================================================== --- /trunk/psastro/src/psastroArguments.c (revision 36836) +++ /trunk/psastro/src/psastroArguments.c (revision 36837) @@ -60,6 +60,9 @@ psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.USE.MODEL", PS_META_REPLACE, "", true); } - // define the reference astrometry file + // define the reference astrometry file (add a container ASTROM.MODEL to config->arguments if -astrommodel (file) is found) status = pmConfigFileSetsMD (config->arguments, &argc, argv, "ASTROM.MODEL", "-astrommodel", "-astrommodellist"); + + // define the koppenhoefer correction file (add a container KH.CORRECT to config->arguments if -kh-correct (file) is found) + status = pmConfigFileSetsMD (config->arguments, &argc, argv, "KH.CORRECT", "-kh-correct", NULL); // define the reference astrometry file @@ -108,4 +111,7 @@ } psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.SKIP.ASTRO", PS_META_REPLACE, "", true); + + // SKIP.ASTRO and USE.MODEL are fundamentally incompatible + psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.USE.MODEL", PS_META_REPLACE, "", false); } Index: /trunk/psastro/src/psastroCleanup.c =================================================================== --- /trunk/psastro/src/psastroCleanup.c (revision 36836) +++ /trunk/psastro/src/psastroCleanup.c (revision 36837) @@ -17,4 +17,5 @@ psFree (config); pmVisualClose (); + pmVisualCleanup (); psTimerStop (); Index: /trunk/psastro/src/psastroConvert.c =================================================================== --- /trunk/psastro/src/psastroConvert.c (revision 36836) +++ /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; + +} + Index: /trunk/psastro/src/psastroDataLoad.c =================================================================== --- /trunk/psastro/src/psastroDataLoad.c (revision 36836) +++ /trunk/psastro/src/psastroDataLoad.c (revision 36837) @@ -53,4 +53,10 @@ pmFPAfileActivate (config->files, true, "PSASTRO.INPUT"); + // if we request KH Correction, activate that file as well + bool applyKH = psMetadataLookupBool (NULL, recipe, "KH.CORRECT.APPLY.EXP"); + if (applyKH) { + pmFPAfileActivate (config->files, true, "PSASTRO.KH.CORRECT"); + } + pmFPAview *view = pmFPAviewAlloc (0); @@ -73,5 +79,5 @@ if (!readout->data_exists) { continue; } - if (!psastroConvertReadout (readout, recipe)) ESCAPE; + if (!psastroConvertReadout (config, view, readout, recipe)) ESCAPE; if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; Index: /trunk/psastro/src/psastroDefineFiles.c =================================================================== --- /trunk/psastro/src/psastroDefineFiles.c (revision 36836) +++ /trunk/psastro/src/psastroDefineFiles.c (revision 36837) @@ -43,7 +43,44 @@ if (!psastroDefineFile (config, input->fpa, "PSASTRO.MODEL", "ASTROM.MODEL", PM_FPA_FILE_ASTROM_MODEL, PM_DETREND_TYPE_ASTROM)) { psError (PS_ERR_IO, false, "Can't find an astrometry model file"); - return NULL; - } - } + return false; + } + } + + bool KHapply = psMetadataLookupBool (&status, recipe, "KH.CORRECT.APPLY"); + if (KHapply) { + // Get the time from FPA.TIME + psTime *time = psMetadataLookupPtr(NULL, input->fpa->concepts, "FPA.TIME"); + if (time->sec == 0 && time->nsec == 0) { + psError (PM_ERR_CONFIG, false, "KH.CORRECT.APPLY requested, but date/time of exposure not found"); + return false; + } + // what is the appropriate date range for the correction + char *KHendDate = psMetadataLookupStr (&status, recipe, "KH.CORRECT.ENDDATE"); + if (!KHendDate) { + psError (PM_ERR_CONFIG, false, "KH.CORRECT.APPLY requested, but KH.CORRECT.ENDDATE not in recipe"); + return false; + } + psTime *endtime = psTimeFromString (KHendDate, PS_TIME_TAI); + if (!endtime) { + psError (PM_ERR_CONFIG, false, "KH.CORRECT.APPLY requested, KH.CORRECT.ENDDATE has an invalid date/time"); + return false; + } + if (psTimeDelta (endtime, time) < 0) { + KHapply = false; + psLogMsg ("psastro", PS_LOG_INFO, "Koppenhoefer correction requested, but exposure after end date, correction skipped"); + } else { + psLogMsg ("psastro", PS_LOG_INFO, "Koppenhoefer correction requested, exposure in valid date range, correction may be applied"); + } + psFree (endtime); + } + if (KHapply) { + if (!psastroDefineFile (config, input->fpa, "PSASTRO.KH.CORRECT", "KH.CORRECT", PM_FPA_FILE_KH_CORRECT, PM_DETREND_TYPE_KH_CORRECT)) { + psError (PS_ERR_IO, false, "Can't find a Koppenhoefer Correction file"); + return NULL; + } + } + // the bool above in the recipe says "apply generally for this camera". here we are asserting that + // we can and should apply for this exposure XXX change the name? + psMetadataAddBool (recipe, PS_LIST_TAIL, "KH.CORRECT.APPLY.EXP", PS_META_REPLACE, "", KHapply); // this step is optional Index: /trunk/psastro/src/psastroLoadCrosstalk.c =================================================================== --- /trunk/psastro/src/psastroLoadCrosstalk.c (revision 36836) +++ /trunk/psastro/src/psastroLoadCrosstalk.c (revision 36837) @@ -175,5 +175,5 @@ float x_t_chip,y_t_chip; - // int faint_ct = 0; + // int faint_ct = 0; NOTE: not currently used pmCellCoordsForChip(&x_cell,&y_cell,CTsourceCell,ref->chip->x,ref->chip->y); @@ -229,5 +229,5 @@ Xt = 3; Ut = 2; - // faint_ct = 1; + // faint_ct = 1; NOTE: not currently used } } @@ -251,5 +251,5 @@ Xt = 4; Ut = 2; - // faint_ct = 1; + // faint_ct = 1; NOTE: not currently used } } Index: /trunk/psastro/src/psastroOneChipFit.c =================================================================== --- /trunk/psastro/src/psastroOneChipFit.c (revision 36836) +++ /trunk/psastro/src/psastroOneChipFit.c (revision 36837) @@ -62,4 +62,5 @@ if (match == NULL) { psLogMsg ("psastro", 3, "failed to find radius-matched sources\n"); + psastroChipFailureHeader (updates); return false; } @@ -69,4 +70,5 @@ if (!unique) { psLogMsg ("psastro", 3, "failed to generate a uniq set of matched sources\n"); + psastroChipFailureHeader (updates); return false; } @@ -99,4 +101,5 @@ psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); psFree (match); + psastroChipFailureHeader (updates); return false; } @@ -126,4 +129,5 @@ psFree (match); psFree (fitStats); + psastroChipFailureHeader (updates); return false; } @@ -185,6 +189,5 @@ psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", astNstar); } else { - psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0); - psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); + psastroChipFailureHeader (updates); } psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX", PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars @@ -222,4 +225,12 @@ } +bool psastroChipFailureHeader (psMetadata *updates) { + psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0); + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); + return true; +} + + + // psastroWriteStars ("raw.1.dat", rawstars); // psastroWriteStars ("ref.1.dat", refstars); Index: /trunk/psastro/src/psastroOneChipGrid.c =================================================================== --- /trunk/psastro/src/psastroOneChipGrid.c (revision 36836) +++ /trunk/psastro/src/psastroOneChipGrid.c (revision 36837) @@ -57,4 +57,5 @@ psFree (rawGridStars); psFree (refGridStars); + psastroChipFailureHeader (updates); return false; } @@ -68,4 +69,5 @@ psFree (rawGridStars); psFree (refGridStars); + psastroChipFailureHeader (updates); return false; }