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.