Index: trunk/psastro/src/psastroLoadGhosts.c =================================================================== --- trunk/psastro/src/psastroLoadGhosts.c (revision 41290) +++ trunk/psastro/src/psastroLoadGhosts.c (revision 41338) @@ -110,5 +110,14 @@ psPolynomial1D *innerMinor = NULL; + 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) @@ -119,12 +128,4 @@ } pmFPA *fpa = astrom->fpa; - - 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); // raise an error if the config is broken @@ -202,4 +203,28 @@ ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, ghost->FP->x, ghost->FP->y); // fprintf (stderr, "-> model chip position: %f, %f\n", ghost->chip->x, ghost->chip->y); + + psPlane coordFPA,ptCH2; + psPlane coordCell,ptFP2; + coordCell.x = ghost->chip->x; + coordCell.y = ghost->chip->y; + psPlaneTransformApply(&coordFPA,ghostChip->toFPA,&coordCell); + ptCH2.x = ghost->chip->x; + ptCH2.y = ghost->chip->y; + psPlaneTransformApply(&ptFP2,ghostChip->toFPA,&ptCH2); + + psPlane ptCH, ptFP; + ptFP.x = ghost->FP->x; + ptFP.y = ghost->FP->y; + psPlaneTransformApply (&ptCH, ghostChip->fromFPA, &ptFP); + + psastroUpdateChipToFPA (fpa, ghostChip); + psPlane ptCH3, ptFP3; + ptCH3.x = ghost->chip->x; + ptCH3.y = ghost->chip->y; + psPlaneTransformApply(&ptFP3,ghostChip->toFPA,&ptCH3); + + + + psLogMsg ("psastro", 3, "Ghost refpos: %f %f %s pos: %f %f %f %f %f %f %f %f %f %f %f %f mag: %f inshape: %f %f %f outshape: %f %f %f \n", ref->FP->x,ref->FP->y,psMetadataLookupStr(NULL,chip->concepts,"CHIP.NAME"), ghost->FP->x,ghost->FP->y,coordFPA.x,coordFPA.y,ghost->chip->x, ghost->chip->y,ptCH.x,ptCH.y,ptFP2.x,ptFP2.y,ptFP3.x,ptFP3.y, ghost->Mag,ghost->inner.major,ghost->inner.minor,ghost->inner.theta/PM_RAD_DEG,ghost->outer.major,ghost->outer.minor,ghost->outer.theta/PM_RAD_DEG); if (ghostChip) { Index: trunk/psastro/src/psastroLoadGlints.c =================================================================== --- trunk/psastro/src/psastroLoadGlints.c (revision 41290) +++ trunk/psastro/src/psastroLoadGlints.c (revision 41338) @@ -43,7 +43,10 @@ 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"); double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); - // select the set of glint regions (GLINT.REGION is a MULTI of METADATA items) @@ -88,7 +91,9 @@ psastroChipBounds (fpa); - // find the possible glint star:s for each glint star, convert the position to FPA coordinates. + + // find the possible glint stars, and 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... @@ -115,11 +120,9 @@ psRegion glintRegion = psRegionFromString (glintRegionString); - // select stars that land in this regio + // 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; - - double glintLength = GLINT_LENGTH_MAG_SLOPE*(GLINT_LENGTH_MAG_ZERO - star->Mag); char *glintType = psMetadataLookupStr (&status, glintItem->data.md, "GLINT.TYPE"); @@ -129,175 +132,318 @@ } - // depending on the glint type, we need to find either the chips in the row or in the column. + double glintLength = GLINT_LENGTH_MAG_SLOPE*(GLINT_LENGTH_MAG_ZERO - star->Mag); + double glintAngle = 0; + + //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)); + } + //also 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.)); + } + + + psVector *x_glint = psVectorAlloc(2,PS_TYPE_F32); + psVector *y_glint = psVectorAlloc(2,PS_TYPE_F32); + 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); + // We want to find the coordinates of the glint end points. However, the glint is straight off the pixel focal plane and has an angle only on the focal plane. So, first find the edge + double xFPA0 = star->FP->x; + double yFPA0 = star->FP->y; + double xFPA1; + double yFPA1; + //angles for TOP and LEFT have been flipped in the fitting + if (!strcasecmp(glintType, "TOP")) { + //find the edge chip and determine the FPA coords of the edge. Then, grab the edge coord and new glint length + 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, 20000.)) { + continue; + } + if (!psastroFindChipInYrange (fpa, nChip, star->FP->x, 20000.)) { + continue; + } + + // FPA coordinates of intersections with chip edges + double yFPAs, yFPAe; + psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip); + + xFPA0 = star->FP->x; + yFPA0 = yFPAe; + glintLength-= (star->FP->y-yFPAe); + } + + xFPA1 = xFPA0 + glintLength*sin(glintAngle*-1.); + yFPA1 = yFPA0 - glintLength*cos(glintAngle*-1.); + + } else { + //find the edge chip and determine the FPA coords of the edge. Then, grab the edge coord and new glint length + 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, -20000.)) { + continue; + } + if (!psastroFindChipInYrange (fpa, nChip, star->FP->x, -20000.)) { + continue; + } + + // FPA coordinates of intersections with chip edges + double yFPAs, yFPAe; + psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip); + + xFPA0 = star->FP->x; + yFPA0 = yFPAs; + glintLength-= (yFPAs-star->FP->y); + } + + xFPA1 = xFPA0 - glintLength*sin(glintAngle); + yFPA1 = yFPA0 + glintLength*cos(glintAngle); + } + + x_glint->data.F32[0] = xFPA0; + y_glint->data.F32[0] = yFPA0; + x_glint->data.F32[1] = xFPA1; + y_glint->data.F32[1] = yFPA1; + + //we need to loop over each corner to select the chips that can have the glint on it + for (int glint_point = 0; glint_point < 2; glint_point++) { + for (int nChip = 0; nChip < fpa->chips->n; nChip++) { + + pmChip *chip = fpa->chips->data[nChip]; + if (!chip) continue; + + if (!psastroFindChipInXrange (fpa, nChip, x_glint->data.F32[glint_point], y_glint->data.F32[glint_point])) { + continue; + } + + // FPA coordinates of intersections with chip edges + double yFPAs, yFPAe; + psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip); + + if (yFPAs > yFPAe) PS_SWAP (yFPAs, yFPAe); + if (yFPA0 > yFPA1) { + PS_SWAP (xFPA0, xFPA1); + PS_SWAP (yFPA0, yFPA1); + } + + // does this glint cross this chip? + if (yFPA0 > yFPAe) continue; + if (yFPA1 < yFPAs) continue; + + //find the y-coord positions for this chip + double ycFPA0, ycFPA1; + ycFPA0 = PS_MAX (yFPA0, yFPAs); + ycFPA1 = PS_MIN (yFPA1, yFPAe); + + //now calculate the proper x-coord positions given the angle, for this chip + double xcFPA0, xcFPA1, angle; + angle = atan2(xFPA1 - xFPA0,yFPA1 - yFPA0); + xcFPA0 = xFPA0 + (ycFPA0-yFPA0)*tan(angle); + xcFPA1 = xFPA0 + (ycFPA1-yFPA0)*tan(angle); + + + double xChip0, yChip0, xChip1, yChip1, chip_angle, glint_length, xChip2,yChip2; + psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, xcFPA0, ycFPA0); + psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, xcFPA1, ycFPA1); + psastroFPAtoChip (&xChip2, &yChip2, fpa, nChip, star->FP->x,star->FP->y); + + 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"); + 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 : %d %f,%f to %f,%f (%f %f %f)\n", glintType, nChip, xChip0, yChip0, xChip1, yChip1, glint_length, glintWidth, 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; + glint->data.F32[3] = glintWidth; + glint->data.F32[4] = chip_angle; + + psArrayAdd (glints, 100, glint); + + psFree (glint); + } } } - // 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); + // We want to find the coordinates of the glint end points. However, the glint is straight off the pixel focal plane and has an angle only on the focal plane. So, first find the edge + double xFPA0 = star->FP->x; + double yFPA0 = star->FP->y; + double xFPA1; + double yFPA1; + //angles for TOP and LEFT have been flipped in the fitting + if (!strcasecmp(glintType, "RIGHT")) { + //find the edge chip and determine the FPA coords of the edge. Then, grab the edge coord and new glint length + for (int nChip = 0; nChip < fpa->chips->n; nChip++) { + pmChip *chip = fpa->chips->data[nChip]; + if (!chip) continue; + + if (!psastroFindChipInXrange (fpa, nChip, 20000.,star->FP->y)) { + continue; + } + if (!psastroFindChipInYrange (fpa, nChip, 20000.,star->FP->y)) { + continue; + } + + // FPA coordinates of intersections with chip edges + double yFPAs, yFPAe; + psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip); + + xFPA0 = yFPAe; + yFPA0 = star->FP->y; + glintLength-= (star->FP->x-yFPAe); + } + + xFPA1 = xFPA0 - glintLength*cos(glintAngle); + yFPA1 = yFPA0 - glintLength*sin(glintAngle); + + } else { + //find the edge chip and determine the FPA coords of the edge. Then, grab the edge coord and new glint length + for (int nChip = 0; nChip < fpa->chips->n; nChip++) { + pmChip *chip = fpa->chips->data[nChip]; + if (!chip) continue; + + if (!psastroFindChipInXrange (fpa, nChip, -20000.,star->FP->y)) { + continue; + } + if (!psastroFindChipInYrange (fpa, nChip, -20000.,star->FP->y)) { + continue; + } + + // FPA coordinates of intersections with chip edges + double yFPAs, yFPAe; + psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip); + + xFPA0 = yFPAs; + yFPA0 = star->FP->y; + glintLength-= (yFPAs-star->FP->x); + } + + xFPA1 = xFPA0 + glintLength*cos(glintAngle*-1.); + yFPA1 = yFPA0 + glintLength*sin(glintAngle*-1.); + } + + x_glint->data.F32[0] = xFPA0; + y_glint->data.F32[0] = yFPA0; + x_glint->data.F32[1] = xFPA1; + y_glint->data.F32[1] = yFPA1; + + //we need to loop over each corner to select the chips that can have the glint on it + for (int glint_point = 0; glint_point < 2; 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_glint->data.F32[glint_point], y_glint->data.F32[glint_point])) { + continue; + } + + // FPA coordinates of intersections with chip edges + double xFPAs, xFPAe; + psastroFindChipXedges (&xFPAs, &xFPAe, fpa, nChip); + + if (xFPAs > xFPAe) PS_SWAP (xFPAs, xFPAe); + if (xFPA0 > xFPA1) { + PS_SWAP (xFPA0, xFPA1); + PS_SWAP (yFPA0, yFPA1); + } + + // does this glint cross this chip? + if (xFPA0 > xFPAe) continue; + if (xFPA1 < xFPAs) continue; + + //find the x-coord positions for this chip + double xcFPA0, xcFPA1; + xcFPA0 = PS_MAX (xFPA0, xFPAs); + xcFPA1 = PS_MIN (xFPA1, xFPAe); + + //now calculate the proper x-coord positions given the angle, for this chip + double ycFPA0, ycFPA1, angle; + angle = atan2(xFPA1 - xFPA0,yFPA1 - yFPA0); + ycFPA0 = yFPA0 + (xcFPA0-xFPA0)/tan(angle); + ycFPA1 = yFPA0 + (xcFPA1-xFPA0)/tan(angle); + + double xChip0, yChip0, xChip1, yChip1, chip_angle, glint_length,xChip2,yChip2; + psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, xcFPA0, ycFPA0); + psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, xcFPA1, ycFPA1); + psastroFPAtoChip (&xChip2, &yChip2, fpa, nChip, star->FP->x,star->FP->y); + + 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"); + 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 : %d %f,%f to %f,%f (%f %f %f)\n", glintType, nChip, xChip0, yChip0, xChip1, yChip1, glint_length, glintWidth, 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; + glint->data.F32[3] = glintWidth; + glint->data.F32[4] = chip_angle; + + psArrayAdd (glints, 100, glint); + + psFree (glint); + } } } + + if (!strcasecmp(glintType, "HSC")) { // It's inefficient to keep looking these up. Index: trunk/psastro/src/psastroMaskUpdates.c =================================================================== --- trunk/psastro/src/psastroMaskUpdates.c (revision 41290) +++ trunk/psastro/src/psastroMaskUpdates.c (revision 41338) @@ -134,6 +134,7 @@ 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"); double REFSTAR_MASK_SATSTAR_MAG_MAX = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_MAG_MAX"); + double REFSTAR_MASK_SATSTAR_RAD_OFFSET = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_RAD_OFFSET"); + double REFSTAR_MASK_SATSTAR_EXP = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_EXP"); double REFSTAR_MASK_SATSTAR_POS_ZERO = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_POS_ZERO"); double REFSTAR_MASK_SATSPIKE_MAG_SLOPE = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSPIKE_MAG_SLOPE"); @@ -232,4 +233,14 @@ double ROTANGLE = PM_RAD_DEG * psMetadataLookupF64 (&status, fpa->concepts, "FPA.ROTANGLE"); psAssert (status, "ROTANGLE missing"); + int ROT_PARITY = psMetadataLookupS32(&status, recipe, "PSASTRO.MODEL.ROT.PARITY"); + if (!status) psAbort ("Can't find recipe option PSASTRO.MODEL.ROT.PARITY"); + + //select the astromeric rotation angle from the header + psMetadata *header = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER"); + if (!header) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO header"); + return false; + } + float AST_T0 = PM_RAD_DEG * psMetadataLookupF32(&status, header, "AST_T0"); // de-activate all files except PSASTRO.INPUT.MASK and PSASTRO.OUTPUT.MASK @@ -327,5 +338,5 @@ // CIRCLE around the stars (scaled by magnitude) - float radius = REFSTAR_MASK_SATSTAR_MAG_SLOPE * (REFSTAR_MASK_SATSTAR_MAG_MAX - ref->Mag); + float radius = REFSTAR_MASK_SATSTAR_RAD_OFFSET + pow(REFSTAR_MASK_SATSTAR_EXP,(REFSTAR_MASK_SATSTAR_MAG_MAX - ref->Mag)); // XXX for now, assume cell binning is 1x1 relative to chip @@ -333,9 +344,8 @@ for (float theta = 0.0; theta < 2*M_PI; theta += M_PI / 2.0) { - - float Theta = theta - ROTANGLE - REFSTAR_MASK_SATSTAR_POS_ZERO; + float Theta = theta - (ROT_PARITY * ROTANGLE) - AST_T0 - REFSTAR_MASK_SATSTAR_POS_ZERO; // LINE for boundaries of the saturation spikes (scaled by magnitude) - // float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2; + //float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2; float MAG_MAX = REFSTAR_MASK_SATSPIKE_MAG_MAX; @@ -400,4 +410,5 @@ } } + } @@ -412,4 +423,7 @@ for (int i = 0; i < ghosts->n; i++) { psastroGhost *ghost = ghosts->data[i]; + //offset the ghost angles with the astrometric rotation angle + ghost->inner.theta -= (ROT_PARITY * ROTANGLE) + ghost->outer.theta -= (ROT_PARITY * ROTANGLE) // XXX bright vs faint ghost bits? (OR with SUSPECT) psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer); @@ -418,5 +432,5 @@ // Select the glint mask regions for this readout (loaded in - // psastroExtractGlints.c). These glint regions are defined as rectangular + // psastroChooseGlintStars.c). These glint regions are defined as rectangular // boxes and are generated for each chip based on the position of the bright // stars beyond the edge of the focal plane. This masking is currently very @@ -426,6 +440,11 @@ // mask the glints on this readout for (int i = 0; i < glints->n; i++) { - psRegion *glint = glints->data[i]; - psastroMaskRectangle (readoutMask->mask, glintMaskValue, glint->x0, glint->y0, glint->x1, glint->y1); + psVector *glint = glints->data[i]; + //psLogMsg ("psastro", 3, "glint: %f %f %f %f\n", glint->x0, glint->y0, glint->x1, glint->y1); + //psastroMaskRectangle (readoutMask->mask, glintMaskValue, glint->x0, glint->y0, glint->x1, glint->y1); + + psastroMaskBox (readoutMask->mask, glintMaskValue, glint->data.F32[0], glint->data.F32[1], + glint->data.F32[2], glint->data.F32[3], glint->data.F32[4]); + } }