Index: /trunk/psastro/src/psastro.h =================================================================== --- /trunk/psastro/src/psastro.h (revision 40489) +++ /trunk/psastro/src/psastro.h (revision 40490) @@ -60,4 +60,5 @@ psastroGhost *psastroGhostAlloc (void); bool psastroLoadGhosts (pmConfig *config); +bool psastroLoadGhostsHSC (pmConfig *config); bool psastroDataSave (pmConfig *config, psMetadata *stats); Index: /trunk/psastro/src/psastroLoadGhosts.c =================================================================== --- /trunk/psastro/src/psastroLoadGhosts.c (revision 40489) +++ /trunk/psastro/src/psastroLoadGhosts.c (revision 40490) @@ -40,4 +40,7 @@ * calculate ghost FPA and Chip positions for the stars loaded on the FPA */ + + + bool psastroLoadGhosts (pmConfig *config) { @@ -58,4 +61,12 @@ psLogMsg ("psastro", PS_LOG_INFO, "determine ghost positions"); + // XXX this needs to be camera-dependent + if (0) { + psWarning("Doing the hard coded switch to HSC ghosts still."); + return (psastroLoadGhostsHSC(config)); + } + + psWarning("hard coded switch to GPC ghosts."); + // select the current recipe psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); @@ -68,4 +79,5 @@ if (!REFSTAR_MASK_GHOST) return true; + char *ghostFile = psMetadataLookupStr (&status, recipe, "GHOST_MODEL"); if (!strcasecmp(ghostFile, "NONE")) return true; @@ -256,2 +268,229 @@ return false; } + + + +bool psastroLoadGhostsHSC (pmConfig *config) { + bool status; + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + float zeropt, exptime, MAX_MAG; + psMetadata *ghostModel = NULL; + psVector *C_terms; + psVector *R_terms; + float glintPixelScale; + + psLogMsg ("psastro", PS_LOG_INFO, "determine ghost 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_GHOST = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_GHOST"); + if (!REFSTAR_MASK_GHOST) return true; + + // Load model from file + char *ghostFile = psMetadataLookupStr (&status, recipe, "GHOST_MODEL"); + if (!strcasecmp(ghostFile, "NONE")) return true; + psLogMsg("psastro", PS_LOG_INFO, "ghost file %s",ghostFile); + if (!pmConfigFileRead (&ghostModel, ghostFile, "GHOST MODEL")) { + psError(PSASTRO_ERR_CONFIG, true, "Trouble loading ghost model"); + return false; + } + + // Allocate FPA and refstars + pmFPAview *view = pmFPAviewAlloc (0); + // 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 = astrom->fpa; + + char *filter = psMetadataLookupStr (&status, fpa->concepts, "FPA.FILTERID"); + psMetadataItem *item = psMetadataLookup(ghostModel, "GHOST.MODEL.HSC"); + if (!item) { + psLogMsg ("psastro", PS_LOG_INFO, "GHOST.MODEL.HSC data missing"); + return false; + } + if (item->type != PS_DATA_METADATA_MULTI) { + psLogMsg ("psastro", PS_LOG_INFO, "GHOST.MODEL.HSC not multi"); + return false; + } + + psListIterator *iter = psListIteratorAlloc(item->data.list,PS_LIST_HEAD, false); + psMetadataItem *refItem = NULL; + while ((refItem = psListGetAndIncrement(iter))) { + if (refItem->type != PS_DATA_METADATA) { + psLogMsg ("psastro", PS_LOG_INFO, "GHOST.MODEL.HSC entry not metadata"); + return false; + } + char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER"); + if (!status) { + // psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); + continue; + } + if (strcmp(refFilter, filter)) continue; + + C_terms = psMetadataLookupPtr(&status, refItem->data.md, "C_TERMS"); + R_terms = psMetadataLookupPtr(&status, refItem->data.md, "R_TERMS"); + glintPixelScale = psMetadataLookupF32(&status, refItem->data.md, "SCALE"); + } + + // really error-out here? or just skip? + if (!psastroZeroPointFromRecipe (&zeropt, &exptime, &MAX_MAG, fpa, 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; + + // this loop selects the matched stars for all chips + while ((chip = pmFPAviewNextChip (view, fpa, 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, 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, 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]; + if (ref->Mag > MAX_MAG) continue; + + psastroGhost *ghost = psastroGhostAlloc (); + + double R_TPA = sqrt(pow(ref->TP->y,2) + pow(ref->TP->x,2)); + double theta0 = atan2(ref->TP->y,ref->TP->x); + double star_radius_deg = R_TPA * glintPixelScale; + + // psTrace("psastro.ghost",5, + psLogMsg("psastro",PS_LOG_INFO, + "Begin ghost %d/%ld: MAX_MAG: %g; ref_mag: %g @ (%.10g,%.10g) TPA: (%f %f) R: %f %f %f", + i,refstars->n,MAX_MAG,ref->Mag,ref->sky->r * 180 / M_PI,ref->sky->d * 180.0 / M_PI,ref->TP->x,ref->TP->y, + R_TPA,star_radius_deg,theta0); + + if (star_radius_deg > 0.5) { continue; } + + // Calculate expected position. + double C = C_terms->data.F32[0] + C_terms->data.F32[1] * star_radius_deg + C_terms->data.F32[2] * pow(star_radius_deg,2) + + C_terms->data.F32[3] * pow(star_radius_deg,3) + C_terms->data.F32[4] * pow(star_radius_deg,4) + + C_terms->data.F32[5] * pow(star_radius_deg,5) + C_terms->data.F32[6] * pow(star_radius_deg,6) + + C_terms->data.F32[7] * pow(star_radius_deg,7); + double R = R_terms->data.F32[0] + R_terms->data.F32[1] * star_radius_deg + R_terms->data.F32[2] * pow(star_radius_deg,2) + + R_terms->data.F32[3] * pow(star_radius_deg,3) + R_terms->data.F32[4] * pow(star_radius_deg,4) + + R_terms->data.F32[5] * pow(star_radius_deg,5) + R_terms->data.F32[6] * pow(star_radius_deg,6) + + R_terms->data.F32[7] * pow(star_radius_deg,7); + R = R / 0.015; + + psPlane *fp = psPlaneAlloc(); + psPlane *tp = psPlaneAlloc(); + + tp->x = 13.5 * (C * cos(theta0) / 0.015 + 12.78); + tp->y = 13.5 * (C * sin(theta0) / 0.015 + 57.74); + psPlaneTransformApply(fp,fpa->fromTPA, tp); + + ghost->srcFP->x = ref->FP->x; + ghost->srcFP->y = ref->FP->y; + ghost->FP->x = fp->x; + ghost->FP->y = fp->y; + + ghost->inner.major = 0.0; + ghost->inner.minor = 0.0; + ghost->outer.major = R; + ghost->outer.minor = R; + + pmChip *ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, ghost->FP->x, ghost->FP->y); + + if (ghostChip) { + psLogMsg("psastro",PS_LOG_INFO, + // psTrace("psastro.ghost",5, + " in ghost: %d/%ld: ref: (%f,%f) or (%s) ghost: (%f,%f) or (%f,%f,%s) %f inner: (%f,%f,%f) outer: (%f,%f,%f)", + i,refstars->n, + ref->FP->x,ref->FP->y, + psMetadataLookupStr(NULL,chip->concepts,"CHIP.NAME"), + ghost->FP->x,ghost->FP->y, + ghost->chip->x,ghost->chip->y,psMetadataLookupStr(NULL,ghostChip->concepts,"CHIP.NAME"), + C, + ghost->inner.major,ghost->inner.minor,ghost->inner.theta, + ghost->outer.major,ghost->outer.minor,ghost->outer.theta); + } + else { + psTrace("psastro.ghost",5," in ghost: %d/%ld: ref: (%f,%f) ghost: (%f,%f) or (%f,%f,%s) inner: (%f,%f,%f) outer: (%f,%f,%f)", + i,refstars->n, + ref->FP->x,ref->FP->y, + ghost->FP->x,ghost->FP->y, + ghost->chip->x,ghost->chip->y,"NONE", + ghost->inner.major,ghost->inner.minor,ghost->inner.theta, + ghost->outer.major,ghost->outer.minor,ghost->outer.theta); + } + + + + if (!ghostChip) goto skip; + if (!ghostChip->cells) goto skip; + if (!ghostChip->cells->n) goto skip; + + + pmCell *ghostCell = ghostChip->cells->data[0]; + if (!ghostCell) goto skip; + if (!ghostCell->readouts) goto skip; + if (!ghostCell->readouts->n) goto skip; + pmReadout *ghostReadout = ghostCell->readouts->data[0]; + if (!ghostReadout) goto skip; + + psArray *ghosts = psMetadataLookupPtr (&status, ghostReadout->analysis, "PSASTRO.GHOSTS"); + if (ghosts == NULL) { + ghosts = psArrayAllocEmpty (100); + if (!psMetadataAdd (ghostReadout->analysis, PS_LIST_TAIL, "PSASTRO.GHOSTS", PS_DATA_ARRAY, "astrometry matches", ghosts)) { + psError(PSASTRO_ERR_CONFIG, false, "failure to add ghosts to readout"); + goto escape; + } + psFree (ghosts); + } + + psArrayAdd (ghosts, 100, ghost); + + skip: + + psFree (ghost); + } + } + } + } + + psastroExtractFreeChipBounds(); + + psFree (C_terms); + psFree (R_terms); + // psFree (ghostModel); + psFree (view); + return true; + + escape: + psFree (C_terms); + psFree (R_terms); + // psFree (ghostModel); + psFree (view); + return false; +} Index: /trunk/psastro/src/psastroLoadGlints.c =================================================================== --- /trunk/psastro/src/psastroLoadGlints.c (revision 40489) +++ /trunk/psastro/src/psastroLoadGlints.c (revision 40490) @@ -44,4 +44,6 @@ double GLINT_LENGTH_MAG_ZERO = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_MAG_ZERO"); double glintWidth = psMetadataLookupF32 (&status, recipe, "GLINT_WIDTH"); + double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); + // select the set of glint regions (GLINT.REGION is a MULTI of METADATA items) @@ -95,5 +97,5 @@ 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); + fprintf (stderr, "glint: %7.2f @ %8.1f, %8.1f (%f %f) %8.1f %8.1f\n", star->Mag, star->FP->x, star->FP->y, star->sky->r * PS_DEG_RAD, star->sky->d * PS_DEG_RAD, star->TP->x,star->TP->y); // find the GLINT.REGION this star lands in (if any) @@ -298,4 +300,232 @@ } } + if (!strcasecmp(glintType, "HSC")) { + // It's inefficient to keep looking these up. + double GLINT_RADIUS_INNER = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.RADIUS.INNER"); + double GLINT_RADIUS_OUTER = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.RADIUS.OUTER"); + + // double R_FPA = sqrt(pow(star->FP->y,2) + pow(star->FP->x,2)); + double R_FPA = sqrt(pow(star->TP->y,2) + pow(star->TP->x,2)); + // R_FPA *= 1.30; + if (R_FPA < GLINT_RADIUS_INNER) { continue; } + if (R_FPA > GLINT_RADIUS_OUTER) { continue; } + psTrace("psastro.masks",4,"HSC_GLINT_STARS: %f %f %f\n", + // star->FP->x,star->FP->y,star->Mag); + star->TP->x,star->TP->y,star->Mag); + psVector *C_terms; + psVector *R_terms; + double scale_factor = 1.0; + + char *filter = psMetadataLookupStr (&status, fpa->concepts, "FPA.FILTERID"); + psMetadataItem *item = psMetadataLookup(glintItem->data.md, "GLINT.FILTER.TERM"); + if (!item) { + psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM data missing"); + return false; + } + if (item->type != PS_DATA_METADATA_MULTI) { + psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM not multi"); + return false; + } + psListIterator *iter = psListIteratorAlloc(item->data.list,PS_LIST_HEAD, false); + psMetadataItem *refItem = NULL; + while ((refItem = psListGetAndIncrement(iter))) { + if (refItem->type != PS_DATA_METADATA) { + psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM entry not metadata"); + return false; + } + char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER"); + if (!status) { + // psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); + continue; + } + if (strcmp(refFilter, filter)) continue; + + C_terms = psMetadataLookupPtr(&status, refItem->data.md, "C_TERMS"); + R_terms = psMetadataLookupPtr(&status, refItem->data.md, "R_TERMS"); + scale_factor = psMetadataLookupF32(&status, refItem->data.md, "FACTOR"); + double limit = psMetadataLookupF32(&status, refItem->data.md, "LIMIT"); + limit += MagOffset; + if (star->Mag > limit) { continue; } + } + // check + double glintPixelScale = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.PIXEL.SCALE"); + // psMetadata *GLINT_PARAMETER_SET = psMetadataLookupPtr(&status, glintItem->data.md, "GLINT.FILTER.TERM"); + // psMetadata *GLINT_PARAMETERS = psMetadataLookupPtr(&status, GLINT_PARAMETER_SET, filter); + + + // double theta0 = atan2(star->FP->y,star->FP->x); + double theta0 = atan2(star->TP->y,star->TP->x); + double star_radius_deg = R_FPA * glintPixelScale; + + double C = C_terms->data.F32[0] + C_terms->data.F32[1] * star_radius_deg + C_terms->data.F32[2] * pow(star_radius_deg,2); + double R = R_terms->data.F32[0] + R_terms->data.F32[1] * star_radius_deg + R_terms->data.F32[2] * pow(star_radius_deg,2); + + double inner_edge_angle, outer_edge_angle; + if (star_radius_deg < 0.8831) { + outer_edge_angle = 0.0; + } + else if (star_radius_deg < 0.9368) { + outer_edge_angle = 46.2985 * sqrt(star_radius_deg - 0.8831); + } + else { + outer_edge_angle = -2.67 + 14.3 * star_radius_deg; + } + if (star_radius_deg < 0.964) { + inner_edge_angle = -243.866 * pow(star_radius_deg - 0.932,2) + 2.4636; + } + else { + inner_edge_angle = 66.2061 * sqrt(star_radius_deg - 0.9635); + } + + inner_edge_angle = inner_edge_angle * PS_RAD_DEG; + outer_edge_angle = outer_edge_angle * PS_RAD_DEG; + + // These define the ends of a single arc in arc-centric coordinates + double x1 = C - R * cos(inner_edge_angle); + double y1 = R * sin(inner_edge_angle); + double x2 = C - R * cos(outer_edge_angle); + double y2 = R * sin(outer_edge_angle); + + // Create the endpoints for the pair of arcs in device coordinates (in millimeters) + double theta1 = theta0 + M_PI / 2.0; + double x1s = (x1 * cos(theta1) - y1 * sin(theta1)) * scale_factor; + double y1s = (y1 * cos(theta1) + x1 * sin(theta1)) * scale_factor; + double x1e = (x2 * cos(theta1) - y2 * sin(theta1)) * scale_factor; + double y1e = (y2 * cos(theta1) + x2 * sin(theta1)) * scale_factor; + + double x2s = (x1 * cos(theta1) + y1 * sin(theta1)) * scale_factor; + double y2s = (-1.0 * y1 * cos(theta1) + x1 * sin(theta1)) * scale_factor; + double x2e = (x2 * cos(theta1) + y2 * sin(theta1)) * scale_factor; + double y2e = (-1.0 * y2 * cos(theta1) + x2 * sin(theta1)) * scale_factor; + + psVector *x_start = psVectorAlloc(4,PS_TYPE_F32); + psVector *y_start = psVectorAlloc(4,PS_TYPE_F32); + psVector *x_end = psVectorAlloc(4,PS_TYPE_F32); + psVector *y_end = psVectorAlloc(4,PS_TYPE_F32); + + // CZW: I know this looks like a typo, but trust me, it's not a typo. + psPlane *fp = psPlaneAlloc(); + psPlane *tp = psPlaneAlloc(); + + fprintf(stderr,"HSC GLINT %f: x1x2(%f %f) (%f %f) x1sx1e (%f %f %f %f) x2sx2e (%f %f %f %f)\n", + star->Mag, + x1,y1,x2,y2, + x1s,y1s,x1e,y1e, + x2s,y2s,x2e,y2e); + + + tp->x = 13.5 * (y1s / 0.015 + 12.78); + tp->y = 13.5 * (x1s / -0.015 + 57.74); + psPlaneTransformApply(fp,fpa->fromTPA, tp); + x_start->data.F32[0] = fp->x; + y_start->data.F32[0] = fp->y; + // x_start->data.F32[0] = tp->x; + // y_start->data.F32[0] = tp->y; + + tp->x = 13.5 * (y1e / 0.015 + 12.78); + tp->y = 13.5 * (x1e / -0.015 + 57.74); + psPlaneTransformApply(fp,fpa->fromTPA, tp); + x_end->data.F32[0] = fp->x; + y_end->data.F32[0] = fp->y; + // x_end->data.F32[0] = tp->x; + // y_end->data.F32[0] = tp->y; + + x_start->data.F32[1] = x_end->data.F32[0]; + y_start->data.F32[1] = y_end->data.F32[0]; + x_end->data.F32[1] = x_start->data.F32[0]; + y_end->data.F32[1] = y_start->data.F32[0]; + + tp->x = 13.5 * (y2s / 0.015 + 12.78); + tp->y = 13.5 * (x2s / -0.015 + 57.74); + psPlaneTransformApply(fp,fpa->fromTPA, tp); + x_start->data.F32[2] = fp->x; + y_start->data.F32[2] = fp->y; + // x_start->data.F32[2] = tp->x; + // y_start->data.F32[2] = tp->y; + + tp->x = 13.5 * (y2e / 0.015 + 12.78); + tp->y = 13.5 * (x2e / -0.015 + 57.74); + psPlaneTransformApply(fp,fpa->fromTPA, tp); + x_end->data.F32[2] = fp->x; + y_end->data.F32[2] = fp->y; + // x_end->data.F32[2] = tp->x; + // y_end->data.F32[2] = tp->y; + + x_start->data.F32[3] = x_end->data.F32[2]; + y_start->data.F32[3] = y_end->data.F32[2]; + x_end->data.F32[3] = x_start->data.F32[2]; + y_end->data.F32[3] = y_start->data.F32[2]; + + psFree(fp); + psFree(tp); + + for (int glint_point = 0; glint_point < 4; glint_point++) { + for (int nChip = 0; nChip < fpa->chips->n; nChip++) { + pmChip *chip = fpa->chips->data[nChip]; + if (!chip) continue; + + if (!psastroFindChipInYrange (fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point])) { + continue; + } + if (!psastroFindChipInXrange (fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point])) { + continue; + } + + double xChip0, yChip0, xChip1, yChip1, chip_angle, glint_length; + psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point]); + psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, x_end->data.F32[glint_point], y_end->data.F32[glint_point]); + + chip_angle = atan2(yChip1 - yChip0, xChip1 - xChip0); + glint_length = sqrt(pow(yChip1 - yChip0,2) + pow(xChip1 - xChip0,2)); + + // 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.HSC"); + if (glints == NULL) { + glints = psArrayAllocEmpty (100); + if (!psMetadataAdd (glintReadout->analysis, PS_LIST_TAIL, "PSASTRO.GLINTS.HSC", PS_DATA_ARRAY, "astrometry matches", glints)) { + psWarning("failure to add glints to readout"); + psFree (glints); + continue; + } + psFree (glints); + } + + fprintf (stderr, "glint %s : %d %f,%f to %f,%f (%f %f %f)\n", glintType, nChip, xChip0, yChip0, xChip1, yChip1, glint_length, glintWidth, chip_angle); + psTrace("psastro.masks",4,"HSC_GLINT: Star: %f %f Glint CR %f %f Parameters: %f %f %f Ends: %f %f -> %f %f Chip: %f %f -> %f %f @ %s %f\n", + star->FP->x,star->FP->y, + C,R, + inner_edge_angle,outer_edge_angle,theta0, + x_start->data.F32[glint_point],y_start->data.F32[glint_point],x_end->data.F32[glint_point],y_end->data.F32[glint_point], + xChip0,yChip0,xChip1,yChip1, + psMetadataLookupStr(&status,glintReadout->parent->parent->concepts,"CHIP.NAME"),chip_angle); + psVector *glint = psVectorAlloc(5,PS_TYPE_F32); + glint->data.F32[0] = xChip0; + glint->data.F32[1] = yChip0; + glint->data.F32[2] = glint_length * pixelScale / pixelScale; + glint->data.F32[3] = glintWidth; + glint->data.F32[4] = chip_angle; + + psArrayAdd (glints, 100, glint); + + psFree (glint); + } // End loop over chips + } // End loop over glint endpoints + + psFree(x_start); + psFree(x_end); + psFree(y_start); + psFree(y_end); + } // End HSC glint block + } } Index: /trunk/psastro/src/psastroMaskUpdates.c =================================================================== --- /trunk/psastro/src/psastroMaskUpdates.c (revision 40489) +++ /trunk/psastro/src/psastroMaskUpdates.c (revision 40490) @@ -129,5 +129,8 @@ bool REFSTAR_MASK_BLEED = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED"); - + bool REFSTAR_MASK_BLEED_ORIENTATION_X = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED_ORIENTATION_X"); + bool REFSTAR_MASK_BLEED_ORIENTATION_Y = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED_ORIENTATION_Y"); + + double REFSTAR_MASK_MAX_MAG = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_MAX_MAG"); double REFSTAR_MASK_SATSTAR_MAG_SLOPE = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_MAG_SLOPE"); @@ -143,4 +146,9 @@ double REFSTAR_MASK_BLEED_MAG_SLOPE = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE"); + double REFSTAR_MASK_BLEED_MAG_MAX_X = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_MAX_X"); + double REFSTAR_MASK_BLEED_MAG_SLOPE_X = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE_X"); + double REFSTAR_MASK_BLEED_MAG_MAX_Y = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_MAX_Y"); + double REFSTAR_MASK_BLEED_MAG_SLOPE_Y = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE_Y"); + double REFSTAR_MASK_CROSSTALK_MAG_MAX = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_MAX"); double REFSTAR_MASK_CROSSTALK_MAG_SLOPE= psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_SLOPE"); @@ -180,4 +188,6 @@ } + psTrace("psastro.masks",2,"Configuration loaded for masking."); + // 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) @@ -187,6 +197,12 @@ REFSTAR_MASK_SATSPIKE_MAG_MAX += MagOffset; REFSTAR_MASK_BLEED_MAG_MAX += MagOffset; + REFSTAR_MASK_BLEED_MAG_MAX_X += MagOffset; + REFSTAR_MASK_BLEED_MAG_MAX_Y += MagOffset; + GHOST_MAX_MAG += MagOffset; + psTrace("psastro.masks",2,"Magnitudes: max_mag: %f satstar: %f satspike: %f bleed: %f ghost %f\n", + REFSTAR_MASK_MAX_MAG,REFSTAR_MASK_SATSTAR_MAG_MAX,REFSTAR_MASK_SATSPIKE_MAG_MAX,REFSTAR_MASK_BLEED_MAG_MAX,GHOST_MAX_MAG); + // select the output mask image :: we mosaic to chip mosaic format pmFPAfile *outMask = psMetadataLookupPtr (NULL, config->files, "PSASTRO.OUTPUT.MASK"); @@ -287,4 +303,5 @@ if (!refstars) continue; + psTrace("psastro.masks",2,"In Readout loop."); // we need to generate the following masks regions: // 1) circle around the saturated stars (scaled by magnitude) @@ -301,5 +318,6 @@ } } - + psTrace("psastro.masks",4,"In refstar loop: %d/%ld %f %f\n", + i,refstars->n,ref->Mag,REFSTAR_MASK_MAX_MAG); if (ref->Mag > REFSTAR_MASK_MAX_MAG) continue; @@ -334,5 +352,7 @@ // The length should also be a function of the image background level - + psTrace("psastro.masks",4,"Masking: Radius: %f Theta: %f Length: %f Width: %f\n", + radius,Theta,spikeLength,spikeWidth); + psTrace("psastro.masks",4,"Also: %f %f %f %f\n",ref->chip->x,ref->chip->y,ROTANGLE,REFSTAR_MASK_SATSTAR_POS_ZERO); psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta); } @@ -349,15 +369,33 @@ float xCell = 0.0; float yCell = 0.0; - float xEnd = 0.0; - float yEnd = 0.0; + float xStart = 0.0; + float xEnd = 0.0; + float yStart = 0.0; + float yEnd = 0.0; + // float xEnd = 0.0; + // float yEnd = 0.0; + float width = 0.0; + float length = 0.0; // find coordinate of star on cell pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y); // find coordinate of end-point on chip - int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE"); - pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize); - - float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag); - psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd); + if (REFSTAR_MASK_BLEED_ORIENTATION_X) { + length = pow(10,REFSTAR_MASK_BLEED_MAG_SLOPE_X*(REFSTAR_MASK_BLEED_MAG_MAX_X - ref->Mag)); + width = REFSTAR_MASK_BLEED_MAG_SLOPE *(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag); + + pmChipCoordsForCell (&xStart, &yStart, refCell, xCell - length, yCell - 0.5 * width); + pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell + length, yCell + 0.5 * width); + psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) xStart, (int) yStart, (int) xEnd, (int) yEnd + 1); + } + if (REFSTAR_MASK_BLEED_ORIENTATION_Y) { + length = pow(10,REFSTAR_MASK_BLEED_MAG_SLOPE_Y*(REFSTAR_MASK_BLEED_MAG_MAX_Y - ref->Mag)); + width = REFSTAR_MASK_BLEED_MAG_SLOPE *(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag); + + pmChipCoordsForCell (&xStart, &yStart, refCell, xCell - 0.5 * width, yCell - length); + pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell + 0.5 * width, yCell + length); + psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) xStart, (int) yStart, (int) xEnd + 1, (int) yEnd); + } + } } @@ -392,5 +430,21 @@ } } - + // Because it's very GPC1-specific, handle alternate glint types here + psArray *hsc_glints = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GLINTS.HSC"); + if (hsc_glints) { + for (int i = 0; i < hsc_glints->n; i++) { + // Contains (start_x, start_y, length, width, theta) + psVector *glint = hsc_glints->data[i]; + + psTrace("psastro.masks", 4, "Masking glint %d on chip %s with params (%f,%f,%f,%f,%f)\n", + i, psMetadataLookupStr(&status,readout->parent->parent->concepts,"CHIP.NAME"), + glint->data.F32[0], glint->data.F32[1], + glint->data.F32[2], glint->data.F32[3], glint->data.F32[4]); + + psastroMaskBox (readoutMask->mask, glintMaskValue, glint->data.F32[0], glint->data.F32[1], + glint->data.F32[2], glint->data.F32[3], glint->data.F32[4]); + } + } + // Select the crosstalk artifact regions for this readout, and mask a circular region // corresponding to the source star