Index: trunk/ippconfig/gpc2/ghost.model.mdc =================================================================== --- trunk/ippconfig/gpc2/ghost.model.mdc (revision 41339) +++ trunk/ippconfig/gpc2/ghost.model.mdc (revision 41339) @@ -0,0 +1,64 @@ +# TdB20200327 ghost model fit for gpc2 +GHOST.CENTER.X METADATA + NORDER_X S32 3 + NORDER_Y S32 3 + VAL_X00_Y00 F64 -1.15871916e+02 + VAL_X01_Y00 F64 4.09678844e-02 + VAL_X02_Y00 F64 2.83256020e-08 + VAL_X03_Y00 F64 5.83797628e-11 + VAL_X00_Y01 F64 2.14720353e-03 + VAL_X01_Y01 F64 -2.07766141e-07 + VAL_X02_Y01 F64 2.45261797e-12 + VAL_X00_Y02 F64 2.45187527e-07 + VAL_X01_Y02 F64 1.24779721e-10 + VAL_X00_Y03 F64 4.74932940e-12 + NELEMENTS S32 10 +END + +GHOST.CENTER.Y METADATA + NORDER_X S32 3 + NORDER_Y S32 3 + VAL_X00_Y00 F64 4.06137992e+01 + VAL_X01_Y00 F64 -3.89362940e-02 + VAL_X02_Y00 F64 -2.22548205e-07 + VAL_X03_Y00 F64 8.83360295e-11 + VAL_X00_Y01 F64 -2.50474486e-03 + VAL_X01_Y01 F64 3.37661246e-07 + VAL_X02_Y01 F64 1.76755513e-10 + VAL_X00_Y02 F64 1.92636172e-07 + VAL_X01_Y02 F64 -6.90714499e-12 + VAL_X00_Y03 F64 1.45561098e-10 + NELEMENTS S32 10 +END + + +# These are the original linear solutions +GHOST.INNER.MAJOR METADATA + NORDER_X S32 1 + VAL_X00 F64 3.926693e+01 + VAL_X01 F64 5.325759e-03 + NELEMENTS S32 2 +END + +GHOST.INNER.MINOR METADATA + NORDER_X S32 1 + VAL_X00 F64 5.287548e+01 + VAL_X01 F64 -2.191669e-03 + NELEMENTS S32 2 +END + +GHOST.OUTER.MAJOR METADATA + NORDER_X S32 2 + VAL_X00 F64 1.56552979e+02 + VAL_X01 F64 4.81347080e-04 + VAL_X02 F64 1.16115995e-06 + NELEMENTS S32 3 +END + +GHOST.OUTER.MINOR METADATA + NORDER_X S32 1 + VAL_X00 F64 1.72399795e+02 + VAL_X01 F64 -5.45331919e-03 + NELEMENTS S32 2 +END + Index: trunk/psModules/src/objects/pmSourceIO_Ghosts.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_Ghosts.c (revision 41339) +++ trunk/psModules/src/objects/pmSourceIO_Ghosts.c (revision 41339) @@ -0,0 +1,496 @@ +/** @file pmSourceIO_Ghosts.c + * + * @author TdB, IfA + * + * @version $Revision: 1.0 $ $Name: not supported by cvs2svn $ + * @date $Date: 20020-03-30 $ + * + * Copyright 2020 University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" +#include "pmModelClass.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSourceSatstar.h" +#include "pmSourceLensing.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" + +#include "pmSourceIO.h" + +#include "pmAstrometryObjects.h" +#include "pmAstrometryWCS.h" + +# include "psastroInternal.h" +static psVector *chipXmin = NULL; +static psVector *chipXmax = NULL; +static psVector *chipYmin = NULL; +static psVector *chipYmax = NULL; + +# define ESCAPE(MSG) { \ + psLogMsg ("psastro", PS_LOG_INFO, MSG); \ + return false; } + +# define GET_2D_POLY(NAME,OUT) \ + md = psMetadataLookupMetadata (&status, ghostModel, NAME); \ + if (!md) { \ + psError(PSASTRO_ERR_CONFIG, true, "Missing %s in model file %s", NAME, ghostFile); \ + goto escape; \ + } \ + OUT = psPolynomial2DfromMetadata(md); \ + if (!OUT) { \ + psError(PSASTRO_ERR_CONFIG, true, "Trouble interpretting %s in model file %s", NAME, ghostFile); \ + goto escape; \ + } + +# define GET_1D_POLY(NAME,OUT) \ + md = psMetadataLookupMetadata (&status, ghostModel, NAME); \ + if (!md) { \ + psError(PSASTRO_ERR_CONFIG, true, "Missing %s in model file %s", NAME, ghostFile); \ + goto escape; \ + } \ + OUT = psPolynomial1DfromMetadata(md); \ + if (!OUT) { \ + psError(PSASTRO_ERR_CONFIG, true, "Trouble interpretting %s in model file %s", NAME, ghostFile); \ + goto escape; \ + } + +/** + * calculate ghost FPA and Chip positions for the stars loaded on the FPA + */ + + +bool pmSourceIO_WriteGhosts (psFits *fits, pmFPA *fpa, pmConfig *config) { + + bool status; + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + float zeropt, exptime, MAX_MAG; + psMetadata *md = NULL; + psPolynomial2D *centerX = NULL; + psPolynomial2D *centerY = NULL; + psPolynomial1D *outerMajor = NULL; + psPolynomial1D *outerMinor = NULL; + psPolynomial1D *innerMajor = NULL; + psPolynomial1D *innerMinor = NULL; + psMetadata *ghostModel = NULL; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe"); + return false; + } + + bool REFSTAR_MASK_GHOST = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_GHOST"); + if (!REFSTAR_MASK_GHOST) return true; + bool GHOST_OUTPUT = psMetadataLookupBool(&status, recipe, "PSASTRO.SAVE.GHOSTS"); + if (!GHOST_OUTPUT) return true; + + char *ghostFile = psMetadataLookupStr (&status, recipe, "GHOST_MODEL"); + if (!strcasecmp(ghostFile, "NONE")) return true; + + if (!pmConfigFileRead (&ghostModel, ghostFile, "GHOST MODEL")) { + psError(PSASTRO_ERR_CONFIG, true, "Trouble loading ghost model"); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); + + GET_2D_POLY ("GHOST.CENTER.X", centerX); + GET_2D_POLY ("GHOST.CENTER.Y", centerY); + + GET_1D_POLY ("GHOST.OUTER.MAJOR", outerMajor); + GET_1D_POLY ("GHOST.OUTER.MINOR", outerMinor); + GET_1D_POLY ("GHOST.INNER.MAJOR", innerMajor); + GET_1D_POLY ("GHOST.INNER.MINOR", innerMinor); + + // select the input astrometry data (also carries the refstars) + pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!astrom) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + goto escape; + } + pmFPA *fpa_ast = astrom->fpa; + + // really error-out here? or just skip? + if (!psastroZeroPointFromRecipe (&zeropt, &exptime, &MAX_MAG, fpa_ast, recipe)) { + psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe"); + goto escape; + } + + // recipe values are given in instrumental magnitudes + // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime) + float MagOffset = zeropt + 2.5*log10(exptime); + MAX_MAG += MagOffset; + + psArray *table = psMetadataLookupPtr (&status, fpa->analysis, "GHOST_POSITIONS"); + psMemIncrRefCounter (table); + if (!table) { + table = psArrayAllocEmpty (0x1000); + + // this loop selects the matched stars for all chips + while ((chip = pmFPAviewNextChip (view, fpa_ast, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + if (!chip->fromFPA) { continue; } + + while ((cell = pmFPAviewNextCell (view, fpa_ast, 1)) != NULL) { + psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, fpa_ast, 1)) != NULL) { + if (! readout->data_exists) { continue; } + + // select the raw objects for this readout (loaded in psastroChooseRefstars.c) + // XXX : note that we place limits on the refstar sample in psastroChooseRefstars.c: + // 1) on chip and 2) < PSASTRO.MAX.NREF. magnitude limits and clump exclusion are only + psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) { continue; } + + // identify the bright stars of interest + for (int i = 0; i < refstars->n; i++) { + pmAstromObj *ref = refstars->data[i]; + psTrace("psastro.ghost",5,"Begin ghost %d/%ld: MAX_MAG: %g; ref_mag: %g @ (%.10g,%.10g)",i,refstars->n,MAX_MAG,ref->Mag,ref->sky->r * 180 / M_PI,ref->sky->d * 180.0 / M_PI); + if (ref->Mag > MAX_MAG) continue; + + psastroGhost *ghost = psastroGhostAlloc (); + ghost->srcFP->x = ref->FP->x; + ghost->srcFP->y = ref->FP->y; + + // XXX it is stupid that this takes -X_fpa,-Y_fpa --> the analysis script is reporting the guess FPA position, and the fit is using that... + ghost->FP->x = -ref->FP->x + psPolynomial2DEval(centerX, -ref->FP->x, -ref->FP->y); + ghost->FP->y = -ref->FP->y + psPolynomial2DEval(centerY, -ref->FP->x, -ref->FP->y); + + float rSrc = hypot (ref->FP->x, ref->FP->y); + + ghost->inner.major = psPolynomial1DEval (innerMajor, rSrc); + ghost->inner.minor = psPolynomial1DEval (innerMinor, rSrc); + ghost->inner.theta = atan2(ref->FP->y, ref->FP->x); + + ghost->outer.major = psPolynomial1DEval (outerMajor, rSrc); + ghost->outer.minor = psPolynomial1DEval (outerMinor, rSrc); + ghost->outer.theta = atan2(ref->FP->y, ref->FP->x); + + // report instrumental ghost star mags + ghost->Mag = ref->Mag - MagOffset; + + // XXX this code yields a single chip: we need to provide results for any chips + // which encompass the full size of the ghost + pmChip *ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, -ghost->srcFP->x, -ghost->srcFP->y); + ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, ghost->FP->x, ghost->FP->y); + + //do a rudimentary check of whether the ghost is on the pixel FPA + if (abs(ghost->FP->y) > 21000.) continue; + if (abs(ghost->FP->x) > 21000.) continue; + + //save the ghost positions + if (ghostChip) { + psMetadata *row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "X_STAR_FPA", PS_DATA_F32, "x coord on focal plane", ref->FP->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_STAR_FPA", PS_DATA_F32, "y coord on focal plane", ref->FP->y); + psMetadataAdd (row, PS_LIST_TAIL, "STAR_CHIP_NAME", PS_DATA_STRING, "star chip name", psMetadataLookupStr(NULL,chip->concepts,"CHIP.NAME")); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF", PS_DATA_F32, "reference star magnitude", ghost->Mag); + psMetadataAdd (row, PS_LIST_TAIL, "X_GHOST_FPA", PS_DATA_F32, "x coord on focal plane", ghost->FP->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_GHOST_FPA", PS_DATA_F32, "y coord on focal plane", ghost->FP->y); + psMetadataAdd (row, PS_LIST_TAIL, "X_GHOST_CHIP", PS_DATA_F32, "x coord on chip", ghost->chip->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_GHOST_CHIP", PS_DATA_F32, "y coord on chip", ghost->chip->y); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_CHIP_NAME",PS_DATA_STRING, "ghost chip name", psMetadataLookupStr(NULL,ghostChip->concepts,"CHIP.NAME")); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_INNER_MAJ",PS_DATA_F32, "inner major axis (pixels)", ghost->inner.major); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_INNER_MIN",PS_DATA_F32, "inner minor axis (pixels)", ghost->inner.minor); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_INNER_ANG",PS_DATA_F32, "inner angle (degrees)", ghost->inner.theta/PM_RAD_DEG); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_OUTER_MAJ",PS_DATA_F32, "outer major axis (pixels)", ghost->outer.major); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_OUTER_MIN",PS_DATA_F32, "outer minor axis (pixels)", ghost->outer.minor); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_OUTER_ANG",PS_DATA_F32, "outer angle (degrees)", ghost->outer.theta/PM_RAD_DEG); + + psArrayAdd (table, 100, row); + psFree (row); + + } + else { + psMetadata *row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "X_STAR_FPA", PS_DATA_F32, "x coord on focal plane", ref->FP->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_STAR_FPA", PS_DATA_F32, "y coord on focal plane", ref->FP->y); + psMetadataAdd (row, PS_LIST_TAIL, "STAR_CHIP_NAME", PS_DATA_STRING, "star chip name", psMetadataLookupStr(NULL,chip->concepts,"CHIP.NAME")); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF", PS_DATA_F32, "reference star magnitude", ghost->Mag); + psMetadataAdd (row, PS_LIST_TAIL, "X_GHOST_FPA", PS_DATA_F32, "x coord on focal plane", ghost->FP->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_GHOST_FPA", PS_DATA_F32, "y coord on focal plane", ghost->FP->y); + psMetadataAdd (row, PS_LIST_TAIL, "X_GHOST_CHIP", PS_DATA_F32, "x coord on chip", ghost->chip->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_GHOST_CHIP", PS_DATA_F32, "y coord on chip", ghost->chip->y); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_CHIP_NAME",PS_DATA_STRING, "ghost chip name", "NONE"); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_INNER_MAJ",PS_DATA_F32, "inner major axis (pixels)", ghost->inner.major); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_INNER_MIN",PS_DATA_F32, "inner minor axis (pixels)", ghost->inner.minor); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_INNER_ANG",PS_DATA_F32, "inner angle (degrees)", ghost->inner.theta/PM_RAD_DEG); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_OUTER_MAJ",PS_DATA_F32, "outer major axis (pixels)", ghost->outer.major); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_OUTER_MIN",PS_DATA_F32, "outer minor axis (pixels)", ghost->outer.minor); + psMetadataAdd (row, PS_LIST_TAIL, "GHOST_OUTER_ANG",PS_DATA_F32, "outer angle (degrees)", ghost->outer.theta/PM_RAD_DEG); + + psArrayAdd (table, 100, row); + psFree (row); + } + + } + } + } + } + + psFree (centerX); + psFree (centerY); + psFree (innerMajor); + psFree (innerMinor); + psFree (outerMajor); + psFree (outerMinor); + psFree (ghostModel); + psFree (view); + return true; + + if (table->n == 0) { + psFree(table); + return true; + } + } + + if (!psFitsWriteTable(fits, NULL, table, "GHOST_POSITIONS")) { + psError(psErrorCodeLast(), false, "writing GHOST_POSITIONS\n"); + psFree(table); + return false; + } + psFree(table); + return true; + +escape: + psFree (centerX); + psFree (centerY); + psFree (innerMajor); + psFree (innerMinor); + psFree (outerMajor); + psFree (outerMinor); + psFree (ghostModel); + psFree (view); + return false; +} + +static void psastroGhostFree (psastroGhost *ghost) { + + if (ghost == NULL) return; + + psFree (ghost->srcFP); + psFree (ghost->FP); + psFree (ghost->chip); + + return; +} + +psastroGhost *psastroGhostAlloc (void) { + + psastroGhost *ghost = (psastroGhost *) psAlloc(sizeof(psastroGhost)); + psMemSetDeallocator(ghost, (psFreeFunc) psastroGhostFree); + + ghost->srcFP = psPlaneAlloc(); + ghost->FP = psPlaneAlloc(); + ghost->chip = psPlaneAlloc(); + + ghost->Mag = 0.0; + + ghost->inner.major = 0.0; + ghost->inner.minor = 0.0; + ghost->inner.theta = 0.0; + + ghost->outer.major = 0.0; + ghost->outer.minor = 0.0; + ghost->outer.theta = 0.0; + + return ghost; +} + + +bool psastroChipBounds (pmFPA *fpa) { + + chipXmin = psVectorAlloc (fpa->chips->n, PS_TYPE_F32); + chipXmax = psVectorAlloc (fpa->chips->n, PS_TYPE_F32); + chipYmin = psVectorAlloc (fpa->chips->n, PS_TYPE_F32); + chipYmax = psVectorAlloc (fpa->chips->n, PS_TYPE_F32); + + // this loop selects the matched stars for all chips + for (int i = 0; i < fpa->chips->n; i++) { + + pmChip *chip = fpa->chips->data[i]; + if (!chip->process || !chip->file_exists) { continue; } + if (!chip->fromFPA) { continue; } + + // determine RA,DEC of 4 corners, use to find RA_MIN,MAX, DEC_MIN,MAX + psRegion *region = pmChipPixels (chip); + psPlane ptCH[4], ptFP; + + ptCH[0].x = region->x0; + ptCH[0].y = region->y0; + ptCH[1].x = region->x1; + ptCH[1].y = region->y0; + ptCH[2].x = region->x1; + ptCH[2].y = region->y1; + ptCH[3].x = region->x0; + ptCH[3].y = region->y1; + psFree (region); + + double Xmin = +FLT_MAX; + double Xmax = -FLT_MAX; + double Ymin = +FLT_MAX; + double Ymax = -FLT_MAX; + + for (int j = 0; j < 4; j++) { + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH[j]); + Xmin = PS_MIN (ptFP.x, Xmin); + Xmax = PS_MAX (ptFP.x, Xmax); + Ymin = PS_MIN (ptFP.y, Ymin); + Ymax = PS_MAX (ptFP.y, Ymax); + } + + // fpa-range for the given chip + chipXmin->data.F32[i] = Xmin; + chipXmax->data.F32[i] = Xmax; + chipYmin->data.F32[i] = Ymin; + chipYmax->data.F32[i] = Ymax; + } + + return true; +} + +pmChip *psastroFindChip (double *xChip, double *yChip, pmFPA *fpa, double xFPA, double yFPA) { + + *xChip = NAN; + *yChip = NAN; + + if (!chipXmin) { + psastroChipBounds (fpa); + } + + for (int i = 0; i < fpa->chips->n; i++) { + + if (xFPA < chipXmin->data.F32[i]) continue; + if (xFPA >= chipXmax->data.F32[i]) continue; + if (yFPA < chipYmin->data.F32[i]) continue; + if (yFPA >= chipYmax->data.F32[i]) continue; + + pmChip *chip = fpa->chips->data[i]; + psRegion *region = pmChipPixels (chip); + + psPlane ptCH, ptFP; + ptFP.x = xFPA; + ptFP.y = yFPA; + psPlaneTransformApply (&ptCH, chip->fromFPA, &ptFP); + + if (ptCH.x < region->x0) goto next_chip; + if (ptCH.x >= region->x1) goto next_chip; + if (ptCH.y < region->y0) goto next_chip; + if (ptCH.y >= region->y1) goto next_chip; + psFree (region); + + *xChip = ptCH.x; + *yChip = ptCH.y; + return chip; + + next_chip: + psFree (region); + } + + return NULL; +} + +// identify chips which land on this column (FP coords) +bool psastroFindChipInXrange (pmFPA *fpa, int nChip, double xFPA, double yFPA) { + + if (!chipXmin || !chipXmax) { + psAbort ("chip bounds not set"); + } + + if (xFPA < chipXmin->data.F32[nChip]) return false; + if (xFPA >= chipXmax->data.F32[nChip]) return false; + return true; +} + +// identify chips which land on this row (FP coords) +bool psastroFindChipInYrange (pmFPA *fpa, int nChip, double xFPA, double yFPA) { + + if (!chipYmin || !chipYmax) { + psAbort ("chip bounds not set"); + } + + if (yFPA < chipYmin->data.F32[nChip]) return false; + if (yFPA >= chipYmax->data.F32[nChip]) return false; + return true; +} + +// return the FPA coordinates of the Y edges of the chip +bool psastroFindChipYedges (double *yFPAs, double *yFPAe, pmFPA *fpa, int nChip) { + + *yFPAs = chipYmin->data.F32[nChip]; + *yFPAe = chipYmax->data.F32[nChip]; + return true; +} + +// return the FPA coordinates of the X edges of the chip +bool psastroFindChipXedges (double *yFPAs, double *yFPAe, pmFPA *fpa, int nChip) { + + *yFPAs = chipXmin->data.F32[nChip]; + *yFPAe = chipXmax->data.F32[nChip]; + return true; +} + +// convert FPA to Chip coordinates +bool psastroFPAtoChip (double *xChip, double *yChip, pmFPA *fpa, int nChip, double xFPA, double yFPA) { + + pmChip *chip = fpa->chips->data[nChip]; + + psPlane ptCH, ptFP; + ptFP.x = xFPA; + ptFP.y = yFPA; + psPlaneTransformApply (&ptCH, chip->fromFPA, &ptFP); + + *xChip = ptCH.x; + *yChip = ptCH.y; + return true; +} + +bool psastroExtractFreeChipBounds () { + + psFree (chipXmin); + psFree (chipXmax); + psFree (chipYmin); + psFree (chipYmax); + return true; +} Index: trunk/psModules/src/objects/pmSourceIO_Glints.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_Glints.c (revision 41339) +++ trunk/psModules/src/objects/pmSourceIO_Glints.c (revision 41339) @@ -0,0 +1,280 @@ +/** @file pmSourceIO_Glints.c + * + * @author TdB, IfA + * + * @version $Revision: 1.0 $ $Name: not supported by cvs2svn $ + * @date $Date: 20020-03-18 $ + * + * Copyright 2020 University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" +#include "pmModelClass.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSourceSatstar.h" +#include "pmSourceLensing.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" + +#include "pmSourceIO.h" + +#include "pmAstrometryObjects.h" +#include "pmAstrometryWCS.h" + +# include "psastroInternal.h" + +# define ESCAPE(MSG) { \ + psLogMsg ("psastro", PS_LOG_INFO, MSG); \ + return false; } + +bool pmSourceIO_WriteGlints (psFits *fits, pmFPA *fpa, pmConfig *config) { + + bool status; + float zeropt, exptime; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe"); + return false; + } + + bool REFSTAR_MASK_GLINTS = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_GLINTS"); + if (!REFSTAR_MASK_GLINTS) return true; + + // select the limiting magnitude + double GLINT_MAX_MAG = psMetadataLookupF32 (&status, recipe, "GLINT_MAX_MAG"); + double GLINT_LENGTH_MAG_SLOPE = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_MAG_SLOPE"); + double GLINT_LENGTH_MAG_ZERO = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_MAG_ZERO"); + double GLINT_LENGTH_POS_SLOPE = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_POS_SLOPE"); + double GLINT_LENGTH_POS_REF = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_POS_REF"); + double GLINT_ANGLE_POS_SLOPE = psMetadataLookupF32 (&status, recipe, "GLINT_ANGLE_POS_SLOPE"); + double GLINT_ANGLE_POS_REF = psMetadataLookupF32 (&status, recipe, "GLINT_ANGLE_POS_REF"); + double glintWidth = psMetadataLookupF32 (&status, recipe, "GLINT_WIDTH"); + bool GLINT_OUTPUT = psMetadataLookupBool(&status, recipe, "PSASTRO.SAVE.GLINTS"); + if (!GLINT_OUTPUT) return true; + + // select the set of glint regions (GLINT.REGION is a MULTI of METADATA items) + psMetadataItem *glintRegions = psMetadataLookup (recipe, "GLINT.REGION"); + if (!status) { + psWarning ("glint masking requested but glint regions are missing (GLINT.REGION)\n"); + return true; + } + if (glintRegions->type != PS_DATA_METADATA_MULTI) { + psWarning ("GLINT.REGION is not a MULTI\n"); + return true; + } + + // select the input astrometry data (also carries the glintStars) + pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!astrom) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; + } + + // recipe values are given in instrumental magnitudes + // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime) + pmFPA *fpa_ast = astrom->fpa; + if (!psastroZeroPointFromRecipe (&zeropt, &exptime, NULL, fpa_ast, recipe)) { + psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe"); + return false; + } + float MagOffset = zeropt + 2.5*log10(exptime); + GLINT_MAX_MAG += MagOffset; + GLINT_LENGTH_MAG_ZERO += MagOffset; + + // select the raw objects for this readout (loaded in psastroExtract.c) + psArray *glintStars = psMetadataLookupPtr (&status, fpa_ast->analysis, "PSASTRO.GLINT.STARS"); + if (glintStars == NULL) { + psLogMsg ("psastro", PS_LOG_INFO, "no glint stars found"); + return false; + } + + psArray *table = psMetadataLookupPtr (&status, fpa->analysis, "GLINT_POSITIONS"); + psMemIncrRefCounter (table); + if (!table) { + table = psArrayAllocEmpty (0x1000); + + for (int i = 0; i < glintStars->n; i++) { + pmAstromObj *star = glintStars->data[i]; + if (star->Mag > GLINT_MAX_MAG) continue; // XXX should not be needed... + + // project glint star to the focal-plane + psProject (star->TP, star->sky, fpa_ast->toSky); + psPlaneTransformApply (star->FP, fpa_ast->fromTPA, star->TP); + + // find the GLINT.REGION this star lands in (if any) + psListIterator *glintIter = psListIteratorAlloc(glintRegions->data.list, PS_LIST_HEAD, false); + psMetadataItem *glintItem = NULL; + while ((glintItem = psListGetAndIncrement (glintIter))) { + if (glintItem->type != PS_DATA_METADATA) { + psWarning ("GLINT.REGION entry is not a metadata folder"); + continue; + } + + char *glintRegionString = psMetadataLookupStr (&status, glintItem->data.md, "REGION"); + if (!glintRegionString) { + // psWarning ("GLINT.REGION entry is missing REGION entry"); + continue; + } + psRegion glintRegion = psRegionFromString (glintRegionString); + + // select stars that land in this region + if (star->FP->x < glintRegion.x0) continue; + if (star->FP->x > glintRegion.x1) continue; + if (star->FP->y < glintRegion.y0) continue; + if (star->FP->y > glintRegion.y1) continue; + + char *glintType = psMetadataLookupStr (&status, glintItem->data.md, "GLINT.TYPE"); + if (!status) { + psWarning ("GLINT.REGION entry is missing TYPE entry"); + continue; + } + + double glintAngle = 0; + //compute the angle of the glint, which depends on position parallel to the FPX + if (!strcasecmp(glintType, "TOP") || !strcasecmp(glintType, "BOTTOM") ){ + glintAngle = PM_RAD_DEG * (GLINT_ANGLE_POS_SLOPE*((GLINT_ANGLE_POS_REF - star->FP->x)/1000.)); + } + if (!strcasecmp(glintType, "LEFT") || !strcasecmp(glintType, "RIGHT") ) { + glintAngle = PM_RAD_DEG * (GLINT_ANGLE_POS_SLOPE*((GLINT_ANGLE_POS_REF - star->FP->y)/1000.)); + } + + double glintLength = GLINT_LENGTH_MAG_SLOPE*(GLINT_LENGTH_MAG_ZERO - star->Mag); + //Besides brightness, the length of the glints also depends on the position of the star compared to the focal plane. But, seemingly only for stars closer than 30k pixels + if ((!strcasecmp(glintType, "TOP") || !strcasecmp(glintType, "BOTTOM")) && abs(star->FP->y) < 30000 ){ + glintLength /= GLINT_LENGTH_POS_SLOPE*(GLINT_LENGTH_POS_REF - abs(star->FP->y)); + + } + if ((!strcasecmp(glintType, "LEFT") || !strcasecmp(glintType, "RIGHT")) && abs(star->FP->x) < 30000 ) { + glintLength /= GLINT_LENGTH_POS_SLOPE*(GLINT_LENGTH_POS_REF - abs(star->FP->x)); + } + + //do a rudimentary check of whether the glint enters the pixel FPA + if (!strcasecmp(glintType, "TOP") && ((star->FP->y - glintLength) > 20500.)) continue; + if (!strcasecmp(glintType, "BOTTOM") && ((star->FP->y + glintLength) < -20500.)) continue; + if (!strcasecmp(glintType, "LEFT") && ((star->FP->x + glintLength) < -20500.)) continue; + if (!strcasecmp(glintType, "RIGHT") && ((star->FP->x - glintLength) > 20500.)) continue; + + //save the glint positions + psMetadata *row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "X_FPA_START", PS_DATA_F32, "x coord on focal plane", star->FP->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA_START", PS_DATA_F32, "y coord on focal plane", star->FP->y); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF", PS_DATA_F32, "reference star magnitude", star->Mag); + psMetadataAdd (row, PS_LIST_TAIL, "GLINT_LENGTH",PS_DATA_F32, "glint length (pixels)", glintLength); + psMetadataAdd (row, PS_LIST_TAIL, "GLINT_WIDTH", PS_DATA_F32, "glint width (pixels)", glintWidth); + psMetadataAdd (row, PS_LIST_TAIL, "GLINT_ANGLE", PS_DATA_F32, "glint angle (degrees)", glintAngle/PM_RAD_DEG); + psMetadataAdd (row, PS_LIST_TAIL, "GLINT_TYPE", PS_DATA_STRING, "glint type", glintType); + + psArrayAdd (table, 100, row); + psFree (row); + } + } + + if (table->n == 0) { + psFree(table); + return true; + } + } + + if (!psFitsWriteTable(fits, NULL, table, "GLINT_POSITIONS")) { + psError(psErrorCodeLast(), false, "writing GLINT_POSITIONS\n"); + psFree(table); + return false; + } + psFree(table); + return true; +} + + +bool psastroZeroPointFromRecipe (float *zeropt, float *exptime, float *ghostMaxMag, pmFPA *fpa, psMetadata *recipe) { + + bool status; + + // select the filter; default to fixed photcode and mag limit otherwise + char *filter = psMetadataLookupStr (&status, fpa->concepts, "FPA.FILTERID"); + if (!status) ESCAPE ("missing FPA.FILTER in concepts"); + + *exptime = psMetadataLookupF32 (&status, fpa->concepts, "FPA.EXPOSURE"); + if (!status) ESCAPE ("missing FPA.EXPOSURE in concepts"); + + // we need to select the PHOTCODE.DATA folder that matches our filter + psMetadataItem *item = psMetadataLookup (recipe, "PHOTCODE.DATA"); + if (!item) ESCAPE ("PHOTCODE.DATA folders missing"); + if (item->type != PS_DATA_METADATA_MULTI) ESCAPE ("PHOTCODE.DATA not a multi"); + + // PHOTCODE.DATA is a multi of metadata items + psListIterator *iter = psListIteratorAlloc(item->data.list, PS_LIST_HEAD, false); + + psMetadataItem *refItem = NULL; + while ((refItem = psListGetAndIncrement (iter))) { + if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder"); + + char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER"); + if (!status) { + // psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); + continue; + } + + // does this entry match the current filter? + if (strcmp (refFilter, filter)) continue; + + psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter); + + *zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT"); + if (!status) { + psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing ZEROPT"); + continue; + } + if (ghostMaxMag) { + *ghostMaxMag = psMetadataLookupF32 (&status, refItem->data.md, "GHOST_MAX_MAG"); + if (!status) { + psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing GHOST_MAX_MAG"); + continue; + } + //MEH null is a pain.. so only log if set + psLogMsg ("psastro", PS_LOG_INFO, "found GHOST_MAX_MAG %f",*ghostMaxMag); + } + //MEH what zpt is set to + psLogMsg ("psastro", PS_LOG_INFO, "found ZEROPT %f",*zeropt); + psFree (iter); + return true; + } + psFree (iter); + return false; +} + + Index: trunk/psModules/src/objects/psastroInternal.h =================================================================== --- trunk/psModules/src/objects/psastroInternal.h (revision 41339) +++ trunk/psModules/src/objects/psastroInternal.h (revision 41339) @@ -0,0 +1,29 @@ +/** @file psastroInternal.h + * + * @brief + * + * @ingroup psastro + * + * @author IfA + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-07 02:03:34 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# ifdef HAVE_CONFIG_H +# include +# endif + +# ifndef PSASTRO_INTERNAL_H +# define PSASTRO_INTERNAL_H + +# include +# include +# include // for strcasecmp +# include // for unlink +# include +# include +# include +# include "psastro.h" + +#endif