Index: /trunk/psastro/src/Makefile.am =================================================================== --- /trunk/psastro/src/Makefile.am (revision 24644) +++ /trunk/psastro/src/Makefile.am (revision 24645) @@ -89,4 +89,6 @@ psastroMaskUpdates.c \ psastroMaskUtils.c \ + psastroChooseGlintStars.c \ + psastroLoadGlints.c \ psastroLoadGhosts.c \ psastroGhostUtils.c \ Index: /trunk/psastro/src/psastroChooseGlintStars.c =================================================================== --- /trunk/psastro/src/psastroChooseGlintStars.c (revision 24645) +++ /trunk/psastro/src/psastroChooseGlintStars.c (revision 24645) @@ -0,0 +1,73 @@ +/** @file psastroChooseGlintStars.c + * + * @brief: Select stars by magnitude which are likely glint sources + * + * @ingroup libpsastro + * + * @author IfA + * @version $Revision: 1.21 $ + * @date $Date: 2009-02-07 02:03:34 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "psastroInternal.h" + +bool psastroChooseGlintStars (pmConfig *config, psArray *refs, const char *source) { + + bool status; + float zeropt, exptime; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!\n"); + return false; + } + + bool REFSTAR_MASK_GLINTS = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_GLINTS"); + if (!REFSTAR_MASK_GLINTS) return true; + + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, source); + if (!input) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data!\n"); + return false; + } + pmFPA *fpa = input->fpa; + + // really error-out here? or just skip? + if (!psastroZeroPointFromRecipe (&zeropt, &exptime, NULL, fpa, recipe)) { + psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe"); + return false; + } + + // select the limiting magnitude + double GLINT_MAX_MAG = psMetadataLookupF32 (&status, recipe, "GLINT_MAX_MAG"); + + // 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); + GLINT_MAX_MAG += MagOffset; + + // There is no selection based on the location of the star. We need the full astrometry + // solution before we can trust the star positions. + + // the refstars is a subset within range of this chip + psArray *glintStars = psArrayAllocEmpty (100); + + // select the reference objects brighter than the cutoff magnitude within range of the FPA + for (int i = 0; i < refs->n; i++) { + pmAstromObj *ref = refs->data[i]; + if (ref->Mag > GLINT_MAX_MAG) continue; + + pmAstromObj *glint = pmAstromObjCopy(ref); + psArrayAdd (glintStars, 100, glint); + psFree (glint); + } + psTrace ("psastro", 4, "Added %ld glint stars\n", glintStars->n); + + psMetadataAdd (fpa->analysis, PS_LIST_TAIL, "PSASTRO.GLINT.STARS", PS_DATA_ARRAY, "possible glint stars", glintStars); + psFree (glintStars); + + return true; +} Index: /trunk/psastro/src/psastroLoadGlints.c =================================================================== --- /trunk/psastro/src/psastroLoadGlints.c (revision 24645) +++ /trunk/psastro/src/psastroLoadGlints.c (revision 24645) @@ -0,0 +1,391 @@ +/** @file psastroLoadGlints.c + * + * @brief calculate glint FPA and Chip positions for the stars loaded on the FPA + * + * @ingroup libpsastro + * + * @author IfA + * @version $Revision: 1.7 $ + * @date $Date: 2009-02-07 02:03:34 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "psastroInternal.h" + +# define ESCAPE(MSG) { \ + psError(PS_ERR_UNKNOWN, false, "I/O failure in psastroMaskUpdate: %s", MSG); \ + psFree (view); \ + return false; \ + } + +/** + * calculate glint FPA and Chip positions for the stars loaded on the FPA + */ +bool psastroLoadGlints (pmConfig *config) { + + bool status; + float zeropt, exptime; + + psLogMsg ("psastro", PS_LOG_INFO, "determine glint positions"); + + // 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 glintWidth = psMetadataLookupF32 (&status, recipe, "GLINT_WIDTH"); + + // 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; + } + pmFPA *fpa = astrom->fpa; + + // really error-out here? or just skip? + if (!psastroZeroPointFromRecipe (&zeropt, &exptime, NULL, fpa, recipe)) { + psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe"); + 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) + 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->analysis, "PSASTRO.GLINT.STARS"); + if (glintStars == NULL) { + psLogMsg ("psastro", PS_LOG_INFO, "no glint stars found"); + return false; + } + + // set up the chip boundary vectors. + psastroChipBounds (fpa); + + // find the possible glint star:s for each glint star, convert the position to FPA coordinates. + // search for stars within the glint regions + 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->toSky); + psPlaneTransformApply (star->FP, fpa->fromTPA, star->TP); + fprintf (stderr, "glint: %7.2f @ %8.1f, %8.1f\n", star->Mag, star->FP->x, star->FP->y); + + // 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 regio + 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; + + double glintLength = GLINT_LENGTH_MAG_SLOPE*(GLINT_LENGTH_MAG_ZERO - star->Mag); + + char *glintType = psMetadataLookupStr (&status, glintItem->data.md, "GLINT.TYPE"); + if (!status) { + psWarning ("GLINT.REGION entry is missing TYPE entry"); + continue; + } + + // depending on the glint type, we need to find either the chips in the row or in the column. + if (!strcasecmp(glintType, "TOP") || !strcasecmp(glintType, "BOTTOM")) { + for (int nChip = 0; nChip < fpa->chips->n; nChip++) { + + pmChip *chip = fpa->chips->data[nChip]; + if (!chip) continue; + + if (!psastroFindChipInXrange (fpa, nChip, star->FP->x, star->FP->y)) { + continue; + } + + // find the coordinate of the end point + double xFPA = star->FP->x; + double yFPA0 = star->FP->y; + double yFPA1; + if (!strcasecmp(glintType, "TOP")) { + yFPA1 = yFPA0 - glintLength; + } else { + yFPA1 = yFPA0 + glintLength; + } + + // FPA coordinates of intersections with chip edges + double yFPAs, yFPAe; + psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip); + + if (yFPA0 > yFPA1) PS_SWAP (yFPA0, yFPA1); + if (yFPAs > yFPAe) PS_SWAP (yFPAs, yFPAe); + + // does this glint cross this chip? + if (yFPA0 > yFPAe) continue; + if (yFPA1 < yFPAs) continue; + + yFPA0 = PS_MAX (yFPA0, yFPAs); + yFPA1 = PS_MIN (yFPA1, yFPAe); + + double xChip0, yChip0, xChip1, yChip1; + psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, xFPA, yFPA0); + psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, xFPA, yFPA1); + + // we now have the location of the glint ends or edge-intersections in chip coordinates + // double-check if this chip overlaps this glint + if (xChip0 > xChip1) PS_SWAP (xChip0, xChip1); + if (yChip0 > yChip1) PS_SWAP (yChip0, yChip1); + + // bounds of this chip + psRegion *region = pmChipPixels (chip); + if (yChip1 < region->y0) continue; + if (yChip0 > region->y1) continue; + + // this glint touches this chip. calculate the start and end + // coordinates on this chip + double xChipS = PS_MAX(xChip0 - 0.5*glintWidth, region->x0); + double xChipE = PS_MIN(xChip1 + 0.5*glintWidth, region->x1); + double yChipS = PS_MAX (yChip0, region->y0); + double yChipE = PS_MIN (yChip1, region->y1); + + // select the 0th readout of the 0th cell for this chip + if (!chip->cells) continue; + if (!chip->cells->n) continue; + pmCell *glintCell = chip->cells->data[0]; + if (!glintCell) continue; + if (!glintCell->readouts) continue; + if (!glintCell->readouts->n) continue; + pmReadout *glintReadout = glintCell->readouts->data[0]; + if (!glintReadout) continue; + + // save the glints on the readout->analysis metadata, creating if needed + psArray *glints = psMetadataLookupPtr (&status, glintReadout->analysis, "PSASTRO.GLINTS"); + if (glints == NULL) { + glints = psArrayAllocEmpty (100); + if (!psMetadataAdd (glintReadout->analysis, PS_LIST_TAIL, "PSASTRO.GLINTS", PS_DATA_ARRAY, "astrometry matches", glints)) { + psWarning("failure to add glints to readout"); + psFree (glints); + continue; + } + psFree (glints); + } + + fprintf (stderr, "glint %s : %f,%f to %f,%f (%f,%f to %f,%f)\n", glintType, xChip0, yChip0, xChip1, yChip1, xChipS, yChipS, xChipE, yChipE); + psRegion *glint = psRegionAlloc(xChipS, xChipE, yChipS, yChipE); + psArrayAdd (glints, 100, glint); + psFree (glint); + psFree (region); + } + } + + // depending on the glint type, we need to find either the chips in the row or in the column. + if (!strcasecmp(glintType, "LEFT") || !strcasecmp(glintType, "RIGHT")) { + for (int nChip = 0; nChip < fpa->chips->n; nChip++) { + + pmChip *chip = fpa->chips->data[nChip]; + if (!chip) continue; + + if (!psastroFindChipInYrange (fpa, nChip, star->FP->x, star->FP->y)) { + continue; + } + + // find the coordinate of the end point + double yFPA = star->FP->y; + double xFPA0 = star->FP->x; + double xFPA1; + if (!strcasecmp(glintType, "RIGHT")) { + xFPA1 = xFPA0 - glintLength; + } else { + xFPA1 = xFPA0 + glintLength; + } + + // FPA coordinates of intersections with chip edges + double xFPAs, xFPAe; + psastroFindChipXedges (&xFPAs, &xFPAe, fpa, nChip); + + if (xFPA0 > xFPA1) PS_SWAP (xFPA0, xFPA1); + if (xFPAs > xFPAe) PS_SWAP (xFPAs, xFPAe); + + // does this glint cross this chip? + if (xFPA0 > xFPAe) continue; + if (xFPA1 < xFPAs) continue; + + xFPA0 = PS_MAX (xFPA0, xFPAs); + xFPA1 = PS_MIN (xFPA1, xFPAe); + + double xChip0, yChip0, xChip1, yChip1; + psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, xFPA0, yFPA); + psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, xFPA1, yFPA); + + // we now have the location of the glint ends or edge-intersections in chip coordinates + // double-check if this chip overlaps this glint + if (xChip0 > xChip1) PS_SWAP (xChip0, xChip1); + if (yChip0 > yChip1) PS_SWAP (yChip0, yChip1); + + // bounds of this chip + psRegion *region = pmChipPixels (chip); + if (xChip1 < region->x0) continue; + if (xChip0 > region->x1) continue; + + // this glint touches this chip. calculate the start and end + // coordinates on this chip + double yChipS = PS_MAX (yChip0 - 0.5*glintWidth, region->y0); + double yChipE = PS_MIN (yChip1 + 0.5*glintWidth, region->y1); + double xChipS = PS_MAX (xChip0, region->x0); + double xChipE = PS_MIN (xChip1, region->x1); + + // select the 0th readout of the 0th cell for this chip + if (!chip->cells) continue; + if (!chip->cells->n) continue; + pmCell *glintCell = chip->cells->data[0]; + if (!glintCell) continue; + if (!glintCell->readouts) continue; + if (!glintCell->readouts->n) continue; + pmReadout *glintReadout = glintCell->readouts->data[0]; + if (!glintReadout) continue; + + // save the glints on the readout->analysis metadata, creating if needed + psArray *glints = psMetadataLookupPtr (&status, glintReadout->analysis, "PSASTRO.GLINTS"); + if (glints == NULL) { + glints = psArrayAllocEmpty (100); + if (!psMetadataAdd (glintReadout->analysis, PS_LIST_TAIL, "PSASTRO.GLINTS", PS_DATA_ARRAY, "astrometry matches", glints)) { + psWarning("failure to add glints to readout"); + psFree (glints); + continue; + } + psFree (glints); + } + + fprintf (stderr, "glint %s : %f,%f to %f,%f (%f,%f to %f,%f)\n", glintType, xChip0, yChip0, xChip1, yChip1, xChipS, yChipS, xChipE, yChipE); + psRegion *glint = psRegionAlloc(xChipS, xChipE, yChipS, yChipE); + psArrayAdd (glints, 100, glint); + psFree (glint); + psFree (region); + } + } + } + } + psastroExtractFreeChipBounds(); + return true; +} + +// XXX need to place the glints on the right analysis... +// psMetadataAdd (fpa->analysis, PS_LIST_TAIL, "PSASTRO.GLINT.STARS", PS_DATA_ARRAY, "possible glint stars", glintStars); +// psFree (glintStars); + + + +# if (0) + + // depending on the glint type, we need to find either the chips in the row or in the column. + if (!strcasecmp(glintType, "LEFT") || !strcasecmp(glintType, "RIGHT")) { + for (int nChip = 0; nChip < fpa->chips->n; nChip++) { + + double xChip0, yChip0; + if (!psastroFindChipInYrange (&xChip0, &yChip0, fpa, nChip, star->FP->x, star->FP->y)) { + continue; + } + + pmChip *chip = fpa->chips->data[nChip]; + + // select the 0th readout of the 0th cell for this chip + if (!chip) continue; + if (!chip->cells) continue; + if (!chip->cells->n) continue; + pmCell *glintCell = chip->cells->data[0]; + if (!glintCell) continue; + if (!glintCell->readouts) continue; + if (!glintCell->readouts->n) continue; + pmReadout *glintReadout = glintCell->readouts->data[0]; + if (!glintReadout) continue; + + // save the glints on the readout->analysis metadata, creating if needed + psArray *glints = psMetadataLookupPtr (&status, glintReadout->analysis, "PSASTRO.GLINTS"); + if (glints == NULL) { + glints = psArrayAllocEmpty (100); + if (!psMetadataAdd (glintReadout->analysis, PS_LIST_TAIL, "PSASTRO.GLINTS", PS_DATA_ARRAY, "astrometry matches", glints)) { + psWarning("failure to add glints to readout"); + continue; + } + psFree (glints); + } + + // bounds of this chip + psRegion *region = pmChipPixels (chip); + + // find the coordinate of the end point + double xChip1, yChip1; + if (!strcasecmp(glintType, "RIGHT")) { + if (!psastroFindChipInYrange (&xChip1, &yChip1, fpa, nChip, star->FP->x - glintLength, star->FP->y)) { + psAbort ("inconsistent chip position result"); + } + } else { + // find the coordinate of the end point + if (!psastroFindChipInYrange (&xChip1, &yChip1, fpa, nChip, star->FP->x + glintLength, star->FP->y)) { + psAbort ("inconsistent chip position result"); + } + } + + // we have the location in chip coordinates of the two glint end-points. + // check if this chip overlaps this glint + if (xChip0 > xChip1) PS_SWAP (xChip0, xChip1); + if (xChip1 < region->x0) continue; + if (xChip0 > region->x1) continue; + + // this glint touches this chip. calculate the start and end + // coordinates on this chip + double yChip; + double xChipS = PS_MAX (xChip0, region->x0); + double xChipE = PS_MIN (xChip1, region->x1); + + // if the line has any tilt (in chip coordinates), interpolate for Y: + if (fabs(region->y1 - region->y0) > 1.0) { + double yChipS = yChip0 + (yChip1 - yChip0) * (xChipS - xChip0) / (xChip1 - xChip0); + double yChipE = yChip0 + (yChip1 - yChip0) * (xChipE - xChip0) / (xChip1 - xChip0); + yChip = 0.5*(yChipS + yChipE); + } else { + yChip = 0.5*(yChip0 + yChip1); + } + + fprintf (stderr, "glint %s : %f,%f to %f,%f (%f - %f @ %f)\n", glintType, xChip0, yChip0, xChip1, yChip1, xChipS, xChipE, yChip); + psRegion *glint = psRegionAlloc(xChipS, xChipE, yChip - 0.5*glintWidth, yChip + 0.5*glintWidth); + psArrayAdd (glints, 100, glint); + psFree (glint); + psFree (region); + } + } +# endif