Index: /branches/eam_branches/ipp-20140423/psModules/src/astrom/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/astrom/Makefile.am (revision 36817) +++ /branches/eam_branches/ipp-20140423/psModules/src/astrom/Makefile.am (revision 36818) @@ -11,5 +11,6 @@ pmAstrometryRefstars.c \ pmAstrometryWCS.c \ - pmAstrometryVisual.c + pmAstrometryVisual.c \ + pmKHcorrect.c pkginclude_HEADERS = \ @@ -21,5 +22,6 @@ pmAstrometryRefstars.h \ pmAstrometryWCS.h \ - pmAstrometryVisual.h + pmAstrometryVisual.h \ + pmKHcorrect.h CLEANFILES = *~ Index: /branches/eam_branches/ipp-20140423/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/astrom/pmAstrometryObjects.h (revision 36817) +++ /branches/eam_branches/ipp-20140423/psModules/src/astrom/pmAstrometryObjects.h (revision 36818) @@ -42,4 +42,5 @@ float Color; ///< object color float dMag; ///< error on object magnitude + float SBinst; ///< surface brightness, used for Koppenhoefer correction } pmAstromObj; Index: /branches/eam_branches/ipp-20140423/psModules/src/astrom/pmKHcorrect.c =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/astrom/pmKHcorrect.c (revision 36818) +++ /branches/eam_branches/ipp-20140423/psModules/src/astrom/pmKHcorrect.c (revision 36818) @@ -0,0 +1,252 @@ +/** @file pmKHcorrect.c + * @brief Functions to read (and write?) Koppenhoefer correction file + * + * The Koppenhoefer correction is needed for some chips of gpc1 before the camera voltages were adjusted 2011/05/11. + * The correction is a modification of the X (and possibly Y) coordinate of a star which depends on the instrumental + * surface brightness, defined as -2.5 log_10 (DN) + 5.0 log_10 fwhm_maj [XXX be careful about the definition of fwhm_maj] + * + * @ingroup AstroImage + * @author EAM, IfA + * + * Copyright 2014 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ +#include +#include +#include +#include +#include +#include // for unlink +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" +#include "pmFPAExtent.h" +#include "pmFPAfileFitsIO.h" +#include "pmConcepts.h" +#include "pmKHcorrect.h" + +static void KHcorrectDataFree (KHcorrectData *spline) { + + if (!spline) return; + + psFree (spline->xk); + psFree (spline->yk); + psFree (spline->y2); + return; +} + +KHcorrectData *KHcorrectDataAlloc (int Nrow) { + + // allocate xk[Nrow], etc + + KHcorrectData *spline = (KHcorrectData *) psAlloc(sizeof(KHcorrectData)); + psMemSetDeallocator(spline, (psFreeFunc) KHcorrectDataFree); + + spline->N = Nrow; + spline->xk = (float *) psAlloc(Nrow*sizeof(float)); + spline->yk = (float *) psAlloc(Nrow*sizeof(float)); + spline->y2 = (float *) psAlloc(Nrow*sizeof(float)); + + return spline; +} + +// the KH correction is a function of the instrumental surface brightness, SBinst +float KHcorrectApply (KHcorrectData *spline, float X) { + + int N = spline->N; + + // saturate correction at high and low ends + if (X < spline->xk[ 0]) return spline->yk[ 0]; + if (X > spline->xk[N-1]) return spline->yk[N-1]; + + float *xk = spline->xk; + float *yk = spline->yk; + float *y2 = spline->y2; + + /* find correct element in array (x must be sorted) */ + int lo = 0; + int hi = N-1; + while (hi - lo > 1) { + int i = 0.5*(hi+lo); + if (xk[i] > X) { + hi = i; + } else { + lo = i; + } + } + + /* error condition: duplicate abssisca */ + float dx = xk[hi] - xk[lo]; + if (dx == 0.0) { + return (0.0); + } + + /* evaluate spline */ + float a = (xk[hi] - X) / dx; + float b = (X - xk[lo]) / dx; + + float value = a*yk[lo] + b*yk[hi] + ((a*a - 1.0)*a*y2[lo] + (b*b - 1.0)*b*y2[hi])*(dx*dx) / 6.0; + return (value); +} + +/********************* CheckDataStatus functions *****************************/ + +bool pmKHcorrectCheckDataStatusForView (const pmFPAview *view, pmFPAfile *file) { + psError(PS_ERR_IO, false, "Check Data Status not defined"); + return false; +} + +bool pmKHcorrectCheckDataStatusForFPA (const pmFPA *fpa) { + psError(PS_ERR_IO, false, "Check Data Status not defined"); + return false; +} + +bool pmKHcorrectCheckDataStatusForChip (const pmChip *chip) { + psError(PS_ERR_IO, false, "Check Data Status not defined"); + return false; +} + +/********************* Write Data functions *****************************/ + +// NOTE : these are not exposed because I don't think we need them (we do not create KHcorrect +// in psModules based tool, but in DVO. + +bool pmKHcorrectWriteFPA (pmFPAfile *file, const pmFPA *fpa); +bool pmKHcorrectWriteForView (const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + // write the full model in one pass: require the level to be FPA + if (view->chip != -1) { + psError(PS_ERR_IO, false, "Koppenhoefer Correction must be written at the FPA level"); + return false; + } + + pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); // Suitable FPA for writing + + if (!pmKHcorrectWriteFPA(file, fpa)) { + psError(PS_ERR_IO, false, "Failed to write KH Correction for fpa"); + psFree(fpa); + return false; + } + + psFree(fpa); + + return true; +} + +// write out all chip-level KH Correction data for this FPA +bool pmKHcorrectWriteFPA (pmFPAfile *file, const pmFPA *fpa) +{ + psError(PS_ERR_IO, false, "output for KH Correction is not defined"); + return false; +} + +/********************* Read Data functions *****************************/ + +bool pmKHcorrectReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) + { + // read the full model in one pass: require the level to be FPA + if (view->chip != -1) { + psError(PS_ERR_IO, false, "KH Correction must be read at the FPA level"); + return false; + } + + if (!pmKHcorrectReadFPA (file)) { + psError(PS_ERR_IO, false, "Failed to read KH Correction for fpa"); + return false; + } + return true; + } + +// read in all chip-level KH Correction data for this FPA +bool pmKHcorrectReadFPA (pmFPAfile *file) { + + if (!pmKHcorrectReadChips (file)) { + psError(PS_ERR_IO, false, "Failed to read KH Correction for chips"); + return false; + } + + return true; +} + +// read the set of tables, one for each chip +bool pmKHcorrectReadChips (pmFPAfile *file) { + + bool haveData, status; + + // loop over the extensions + // for each extension, use the extname (eg, XY01.DX.T0) to assign to a chip + + // move to the start of the file + haveData = psFitsMoveExtNum (file->fits, 1, false); + if (!haveData) { + psError(PS_ERR_IO, false, "Failed to read even the first extension?"); + return false; + } + + while (haveData) { + + // load the header + psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header + if (!header) psAbort("cannot read model header"); + + // load the full model in one shot + psArray *model = psFitsReadTable (file->fits); + if (!model) psAbort("cannot read model"); + + // determine the chip: + char *extname = psMetadataLookupStr (&status, header, "EXTNAME"); + psLogMsg ("psModules.astrom", 4, "read %ld rows from Koppenhoefer correction file, extname %s\n", model->n, extname); + + // I expect to find a name of the form: chipName.dir.tset (eg, XY01.DX.T0) + // where chipName like 'XY01' + // dir = 'DX' + // tset = 'T0' + psAssert (strlen(extname) == 10, "invalid extension %s", extname); + psAssert (extname[5] == 'D', "invalid extension %s", extname); + psAssert (extname[6] == 'X', "invalid extension %s", extname); + psAssert (extname[8] == 'T', "invalid extension %s", extname); + psAssert (extname[9] == '0', "invalid extension %s", extname); + + char chipName[5]; + strncpy (chipName, extname, 4); + chipName[4] = 0; + + pmChip *chip = pmConceptsChipFromName (file->fpa, chipName); + if (!chip) psAbort ("invalid chip?"); + + KHcorrectData *spline = KHcorrectDataAlloc (model->n); + + // parse the model entries + for (int i = 0; i < model->n; i++) { + psMetadata *row = model->data[i]; + + spline->xk[i] = psMetadataLookupF32(&status, row, "X_KNOT"); + spline->yk[i] = psMetadataLookupF32(&status, row, "Y_KNOT"); + spline->y2[i] = psMetadataLookupF32(&status, row, "DY2_DX"); + } + psMetadataAddUnknown (chip->analysis, PS_LIST_TAIL, "KH.CORRECT", PS_META_REPLACE, "", spline); + psFree (spline); + + psFree (model); + psFree (header); + + // move to the next extension + haveData = psFitsMoveExtNum (file->fits, 1, true); + } + + return true; +} + Index: /branches/eam_branches/ipp-20140423/psModules/src/astrom/pmKHcorrect.h =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/astrom/pmKHcorrect.h (revision 36818) +++ /branches/eam_branches/ipp-20140423/psModules/src/astrom/pmKHcorrect.h (revision 36818) @@ -0,0 +1,35 @@ +/** @file pmKHcorrect.h + * @brief Functions to read (and write?) Koppenhoefer correction file + * + * The Koppenhoefer correction is needed for some chips of gpc1 before the camera voltages were adjusted 2011/05/11. + * The correction is a modification of the X (and possibly Y) coordinate of a star which depends on the instrumental + * surface brightness, defined as -2.5 log_10 (DN) + 5.0 log_10 fwhm_maj [XXX be careful about the definition of fwhm_maj] + * + * @ingroup AstroImage + * @author EAM, IfA + * + * Copyright 2014 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_KH_CORRECT_H +#define PM_KH_CORRECT_H + +/// @addtogroup Astrometry +/// @{ + +typedef struct { + int N; + float *xk; + float *yk; + float *y2; +} KHcorrectData; + +KHcorrectData *KHcorrectDataAlloc (int Nrow); +float KHcorrectApply (KHcorrectData *spline, float X); + +bool pmKHcorrectReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config); +bool pmKHcorrectReadFPA (pmFPAfile *file); +bool pmKHcorrectReadChips (pmFPAfile *file); + +/// @} +#endif Index: /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfile.c =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfile.c (revision 36817) +++ /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfile.c (revision 36818) @@ -531,7 +531,6 @@ return PM_FPA_FILE_LINEARITY; } - // deprecate this? if (!strcasecmp(type, "ASTROM")) { - return PM_FPA_FILE_ASTROM_MODEL; + return PM_FPA_FILE_ASTROM_MODEL; } if (!strcasecmp(type, "ASTROM.MODEL")) { @@ -540,4 +539,7 @@ if (!strcasecmp(type, "ASTROM.REFSTARS")) { return PM_FPA_FILE_ASTROM_REFSTARS; + } + if (!strcasecmp(type, "KH.CORRECT")) { + return PM_FPA_FILE_KH_CORRECT; } if (!strcasecmp(type, "SUBKERNEL")) { @@ -593,4 +595,6 @@ case PM_FPA_FILE_ASTROM_REFSTARS: return ("ASTROM.REFSTARS"); + case PM_FPA_FILE_KH_CORRECT: + return ("KH.CORRECT"); case PM_FPA_FILE_SUBKERNEL: return ("SUBKERNEL"); Index: /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfile.h =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfile.h (revision 36817) +++ /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfile.h (revision 36818) @@ -48,4 +48,5 @@ PM_FPA_FILE_ASTROM_MODEL, PM_FPA_FILE_ASTROM_REFSTARS, + PM_FPA_FILE_KH_CORRECT, PM_FPA_FILE_SUBKERNEL, PM_FPA_FILE_SRCTEXT, Index: /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfileFitsIO.c =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfileFitsIO.c (revision 36817) +++ /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfileFitsIO.c (revision 36818) @@ -150,5 +150,7 @@ case PM_FPA_FILE_PSF: case PM_FPA_FILE_ASTROM_MODEL: - case PM_FPA_FILE_ASTROM_REFSTARS: { + case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: + { pmHDU *hdu = pmFPAviewThisHDU(view, fpa); if (hdu) { Index: /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfileIO.c =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfileIO.c (revision 36817) +++ /branches/eam_branches/ipp-20140423/psModules/src/camera/pmFPAfileIO.c (revision 36818) @@ -48,4 +48,5 @@ #include "pmPSF_IO.h" +#include "pmKHcorrect.h" #include "pmAstrometryModel.h" #include "pmAstrometryRefstars.h" @@ -234,4 +235,7 @@ status = pmAstromModelReadForView (view, file, config); break; + case PM_FPA_FILE_KH_CORRECT: + status = pmKHcorrectReadForView (view, file, config); + break; case PM_FPA_FILE_EXPNUM: status = pmExpNumRead(view, file, config); @@ -324,4 +328,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_JPEG: case PM_FPA_FILE_KAPA: @@ -408,4 +413,8 @@ } } + if (file->type == PM_FPA_FILE_KH_CORRECT) { + psTrace("psModules.camera", 6, "skip write for %s, no write function defined", file->name); + return true; + } // open the file if not yet opened @@ -500,4 +509,8 @@ case PM_FPA_FILE_ASTROM_REFSTARS: status = pmAstromRefstarsWriteForView (view, file, config); + break; + + case PM_FPA_FILE_KH_CORRECT: + psError(PS_ERR_IO, true, "cannot write type KH.CORRECT (%s)", file->name); break; @@ -576,4 +589,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_LINEARITY: case PM_FPA_FILE_EXPNUM: @@ -653,4 +667,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_EXPNUM: psTrace ("psModules.camera", 6, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name); @@ -815,4 +830,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_LINEARITY: case PM_FPA_FILE_EXPNUM: @@ -1018,4 +1034,5 @@ case PM_FPA_FILE_EXPNUM: case PM_FPA_FILE_ASTROM_MODEL: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_SX: case PM_FPA_FILE_RAW: Index: /branches/eam_branches/ipp-20140423/psModules/src/detrend/pmDetrendDB.c =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/detrend/pmDetrendDB.c (revision 36817) +++ /branches/eam_branches/ipp-20140423/psModules/src/detrend/pmDetrendDB.c (revision 36818) @@ -111,4 +111,5 @@ DETREND_STRING_CASE(LINEARITY); DETREND_STRING_CASE(AUXMASK); + DETREND_STRING_CASE(KH_CORRECT); default: return NULL; Index: /branches/eam_branches/ipp-20140423/psModules/src/detrend/pmDetrendDB.h =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/detrend/pmDetrendDB.h (revision 36817) +++ /branches/eam_branches/ipp-20140423/psModules/src/detrend/pmDetrendDB.h (revision 36818) @@ -40,4 +40,5 @@ PM_DETREND_TYPE_LINEARITY, PM_DETREND_TYPE_AUXMASK, + PM_DETREND_TYPE_KH_CORRECT, } pmDetrendType; Index: /branches/eam_branches/ipp-20140423/psModules/src/objects/pmModel.c =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/objects/pmModel.c (revision 36817) +++ /branches/eam_branches/ipp-20140423/psModules/src/objects/pmModel.c (revision 36818) @@ -39,4 +39,6 @@ { psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__); + if (!tmp) return; + psFree(tmp->params); psFree(tmp->dparams); Index: /branches/eam_branches/ipp-20140423/psModules/src/psmodules.h =================================================================== --- /branches/eam_branches/ipp-20140423/psModules/src/psmodules.h (revision 36817) +++ /branches/eam_branches/ipp-20140423/psModules/src/psmodules.h (revision 36818) @@ -95,4 +95,5 @@ #include #include +#include // the following headers are from psModule:imcombine