IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41338 for trunk


Ignore:
Timestamp:
Apr 16, 2020, 1:39:13 PM (6 years ago)
Author:
tdeboer
Message:

revert to working Ohana build

Location:
trunk
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/ippconfig/gpc2/Makefile.am

    r40348 r41338  
    88        format_mef.config \
    99        format_20180131.config \
     10        ghost.model.mdc \
    1011        ppImage.config \
    1112        ppMerge.config \

  • trunk/ippconfig/gpc2/psastro.config

    r41233 r41338  
    130130DVO.GETSTAR.PHOTCODE        STR      i
    131131DVO.GETSTAR.MAX.RHO         F32      3000.0
    132 DVO.GETSTAR.MIN.MAG         F32      12.0
    133 DVO.GETSTAR.MIN.MAG.INST    F32     -25.0
     132DVO.GETSTAR.MIN.MAG         F32      2.0
     133DVO.GETSTAR.MIN.MAG.INST    F32     -35.0
    134134
    135135# rather than mean, use specified EDGE:
    … …  
    156156  ZEROPT   F32 24.611
    157157  PHOTCODE STR i
    158   GHOST_MAX_MAG                   F32 -25.0
     158  GHOST_MAX_MAG                   F32 -23.5
    159159END
    160160PHOTCODE.DATA METADATA
    … …  
    174174  ZEROPT   F32 26.3
    175175  PHOTCODE STR r
    176   GHOST_MAX_MAG                   F32 -20.0
    177 END
    178 
    179 REFSTAR_MASK                    BOOL FALSE
     176  GHOST_MAX_MAG                   F32 -23.0
     177END
     178
     179REFSTAR_MASK                    BOOL TRUE
    180180REFSTAR_MASK_MAX_MAG            F32 -15.0
    181181REFSTAR_MASK_SATSTAR_MAG_MAX    F32 -15.0
    182182REFSTAR_MASK_SATSTAR_MAG_SLOPE  F32  10.15
    183 # CZW: Original 2010-02-08 position angle offset
    184 REFSTAR_MASK_SATSTAR_POS_ZERO   F32  -0.798
    185 # CZW: Updated 2011-02-18 position angle offset. As above, offset 0.5 degree.
    186 # REFSTAR_MASK_SATSTAR_POS_ZERO   F32  -0.78927
    187 #REFSTAR_MASK_SATSTAR_POS_ZERO   F32  -0.80237
    188 #REFSTAR_MASK_SATSTAR_POS_ZERO   F32  -0.81109
    189 # Length = 10**(SLOPE*(MAG_MAX - InstMag)) - OFFSET
     183
     184#TdB20200214: Position of gpc2 spikes determined on image
     185REFSTAR_MASK_SATSTAR_POS_ZERO   F32  -0.7892
     186#TdB20200225: The rotator angle for gpc2 goes the other way around as gpc1
     187PSASTRO.MODEL.ROT.PARITY      S32      -1
     188
     189#TdB20200229: The size of saturated star circles should be an exponential for gpc2
     190#rad(mag) = REFSTAR_MASK_SATSTAR_RAD_OFFSET + REFSTAR_MASK_SATSTAR_EXP**(REFSTAR_MASK_SATSTAR_MAG_MAX-mag)
     191REFSTAR_MASK_SATSTAR_MAG_MAX    F32 -11.363
     192REFSTAR_MASK_SATSTAR_RAD_OFFSET F32  17.631
     193REFSTAR_MASK_SATSTAR_EXP        F32   1.493
    190194
    191195REFSTAR_MASK_SATSPIKE_MAG_SLOPE F32  0.096
    … …  
    204208
    205209
    206 
    207 
    208210REFSTAR_MASK_BLEED              BOOL FALSE
    209211REFSTAR_MASK_BLEED_MAG_MAX      F32 -15.0
    … …  
    216218
    217219# ghost model file is relative to PATH (eg, ippconfig directory)
    218 #GHOST_MODEL                     STR  gpc1/ghost.model.mdc
     220GHOST_MODEL                     STR  gpc2/ghost.model.mdc
    219221# GHOST_MAX_MAG                   F32 -20.0
    220222
    221 GLINT_LENGTH_MAG_SLOPE          F32 2500.0
    222 GLINT_LENGTH_MAG_ZERO           F32 -20.0
    223 GLINT_WIDTH                     F32  150.0
    224223REFSTAR_COUNT_GHOSTS            BOOL TRUE
    225224REFSTAR_MASK_GHOST              BOOL TRUE
    226 
    227 # GHOST_CENTER_X_0              F32   66.2
    228 # GHOST_CENTER_X_SLOPE          F32   -0.02454
    229 # GHOST_CENTER_Y_0              F32   16.6
    230 # GHOST_CENTER_Y_SLOPE          F32   +0.01565
    231 #
    232 # GHOST_INNER_MAJOR_0           F32   32.9
    233 # GHOST_INNER_MAJOR_SLOPE         F32   +0.00612
    234 # GHOST_INNER_MINOR_0           F32   58.7
    235 # GHOST_INNER_MINOR_SLOPE         F32   -0.00258
    236 #
    237 # GHOST_OUTER_MAJOR_0           F32   90.6
    238 # GHOST_OUTER_MAJOR_SLOPE         F32   +0.01548
    239 # GHOST_OUTER_MINOR_0           F32  133.9
    240 # GHOST_OUTER_MINOR_SLOPE         F32   -0.00288
     225PSASTRO.SAVE.GHOSTS             BOOL TRUE  # save ghost positions as table in output smf file
     226
     227#TdB20200305: The length of glints as function of mag is recomputed for gpc2
     228#glint_length(mag) = GLINT_LENGTH_MAG_SLOPE*(GLINT_LENGTH_MAG_ZERO-mag)
     229GLINT_LENGTH_MAG_SLOPE          F32 1901.0
     230GLINT_LENGTH_MAG_ZERO           F32 -16.8
     231#TdB20200308: The length of glints also depends on position off the focal plane (perpendicular to pixel FPA)
     232#glint_length(pos) = GLINT_LENGTH_POS_SLOPE*(GLINT_LENGTH_POS_REF-pos)
     233GLINT_LENGTH_POS_SLOPE          F32 0.00019
     234GLINT_LENGTH_POS_REF            F32 33821
     235#TdB20200312: The angle that glints have depends on position along the focal plane (parallel to pixel FPA)
     236#glint_angle(pos) = GLINT_ANGLE_SLOPE*(GLINT_ANGLE_POS_REF-(pos/1000))
     237GLINT_ANGLE_POS_SLOPE           F32  0.07903 
     238GLINT_ANGLE_POS_REF             F32 -575.5
     239GLINT_WIDTH                     F32  30.0
     240
     241# save glint positions?
     242PSASTRO.SAVE.GLINTS             BOOL TRUE  # save glint positions as table in output smf file
    241243
    242244PR_RECIPE METADATA
    … …  
    271273REFSTAR_MASK_CROSSTALK_MAG_MAX   F32 -14.47
    272274CROSSTALK_MAX_MAG                F32  20.0
    273 CROSSTALK_MASK                   BOOL TRUE
     275CROSSTALK_MASK                   BOOL FALSE
    274276
    275277PS1_REFERENCE METADATA

  • trunk/ippconfig/recipes/psastro.config

    r41287 r41338  
    150150REFSTAR_MASK_REGIONS            BOOL FALSE
    151151REFSTAR_MASK_MAX_MAG            F32 -15.0
    152 REFSTAR_MASK_SATSTAR_MAG_MAX    F32 -17.0
    153152REFSTAR_MASK_SATSTAR_MAG_SLOPE  F32  10.0
    154153REFSTAR_MASK_SATSTAR_POS_ZERO   F32   0.0
    … …  
    161160REFSTAR_MASK_SATSPIKE_WIDTH_SLOPE F32 0.00
    162161
     162REFSTAR_MASK_SATSTAR_MAG_MAX    F32 -11.363
     163REFSTAR_MASK_SATSTAR_RAD_OFFSET F32  17.631
     164REFSTAR_MASK_SATSTAR_EXP        F32   1.493
     165
     166PSASTRO.SAVE.GHOSTS             BOOL FALSE
    163167REFSTAR_MASK_GHOST              BOOL FALSE
    164168REFSTAR_COUNT_GHOSTS            BOOL FALSE
    … …  
    184188GLINT_LENGTH_MAG_SLOPE          F32 0.0
    185189GLINT_LENGTH_MAG_ZERO           F32 0.0
     190GLINT_LENGTH_POS_SLOPE          F32 0.0
     191GLINT_LENGTH_POS_REF            F32 0.0
     192GLINT_ANGLE_POS_SLOPE           F32 0.0
     193GLINT_ANGLE_POS_REF             F32 0.0
     194# save glint positions?
     195PSASTRO.SAVE.GLINTS             BOOL FALSE
     196
    186197GLINT_WIDTH                     F32 0.0
    187198GLINT.REGION                    MULTI

  • trunk/psModules/src/objects/Makefile.am

    r40491 r41338  
    6060        pmSourceIO_CMF_PS1_DV5.c \
    6161        pmSourceIO_MatchedRefs.c \
     62        pmSourceIO_Glints.c \
     63        pmSourceIO_Ghosts.c \
    6264        pmSourcePlots.c \
    6365        pmSourcePlotPSFModel.c \

  • trunk/psModules/src/objects/pmSourceIO.c

    r38376 r41338  
    701701    if (file->fpa->chips->n == 1) {
    702702        pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config);
     703        pmSourceIO_WriteGlints (file->fits, file->fpa, config);
    703704        return true;
    704705    }
    … …  
    808809
    809810    pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config);
     811    pmSourceIO_WriteGlints (file->fits, file->fpa, config);
    810812    return true;
    811813}

  • trunk/psModules/src/objects/pmSourceIO.h

    r37964 r41338  
    9191bool pmSourceIO_ReadMatchedRefs (psFits *fits, pmFPA *fpa, const pmConfig *config);
    9292
     93bool pmSourceIO_WriteGlints (psFits *fits, pmFPA *fpa, pmConfig *config);
     94
    9395/// @}
    9496# endif /* PM_SOURCE_IO_H */

  • trunk/psastro/src/psastroLoadGhosts.c

    r41290 r41338  
    110110    psPolynomial1D *innerMinor = NULL;
    111111
     112
    112113    pmFPAview *view = pmFPAviewAlloc (0);
     114
     115    GET_2D_POLY ("GHOST.CENTER.X", centerX);
     116    GET_2D_POLY ("GHOST.CENTER.Y", centerY);
     117
     118    GET_1D_POLY ("GHOST.OUTER.MAJOR", outerMajor);
     119    GET_1D_POLY ("GHOST.OUTER.MINOR", outerMinor);
     120    GET_1D_POLY ("GHOST.INNER.MAJOR", innerMajor);
     121    GET_1D_POLY ("GHOST.INNER.MINOR", innerMinor);
    113122
    114123    // select the input astrometry data (also carries the refstars)
    … …  
    119128    }
    120129    pmFPA *fpa = astrom->fpa;
    121 
    122     GET_2D_POLY ("GHOST.CENTER.X", centerX);
    123     GET_2D_POLY ("GHOST.CENTER.Y", centerY);
    124 
    125     GET_1D_POLY ("GHOST.OUTER.MAJOR", outerMajor);
    126     GET_1D_POLY ("GHOST.OUTER.MINOR", outerMinor);
    127     GET_1D_POLY ("GHOST.INNER.MAJOR", innerMajor);
    128     GET_1D_POLY ("GHOST.INNER.MINOR", innerMinor);
    129130
    130131    // raise an error if the config is broken
    … …  
    202203                    ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, ghost->FP->x, ghost->FP->y);
    203204                    // fprintf (stderr, "-> model chip position: %f, %f\n", ghost->chip->x, ghost->chip->y);
     205
     206                    psPlane coordFPA,ptCH2;
     207                    psPlane coordCell,ptFP2;
     208                    coordCell.x = ghost->chip->x;
     209                    coordCell.y = ghost->chip->y;
     210                    psPlaneTransformApply(&coordFPA,ghostChip->toFPA,&coordCell);
     211                    ptCH2.x = ghost->chip->x;
     212                    ptCH2.y = ghost->chip->y;
     213                    psPlaneTransformApply(&ptFP2,ghostChip->toFPA,&ptCH2);
     214
     215                    psPlane ptCH, ptFP;
     216                    ptFP.x = ghost->FP->x;
     217                    ptFP.y = ghost->FP->y;
     218                    psPlaneTransformApply (&ptCH, ghostChip->fromFPA, &ptFP);
     219
     220                    psastroUpdateChipToFPA (fpa, ghostChip);
     221                    psPlane ptCH3, ptFP3;
     222                    ptCH3.x = ghost->chip->x;
     223                    ptCH3.y = ghost->chip->y;
     224                    psPlaneTransformApply(&ptFP3,ghostChip->toFPA,&ptCH3);
     225
     226
     227
     228                    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);
    204229
    205230                    if (ghostChip) {

  • trunk/psastro/src/psastroLoadGlints.c

    r41290 r41338  
    4343    double GLINT_LENGTH_MAG_SLOPE = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_MAG_SLOPE");
    4444    double GLINT_LENGTH_MAG_ZERO = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_MAG_ZERO");
     45    double GLINT_LENGTH_POS_SLOPE = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_POS_SLOPE");
     46    double GLINT_LENGTH_POS_REF = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_POS_REF");
     47    double GLINT_ANGLE_POS_SLOPE = psMetadataLookupF32 (&status, recipe, "GLINT_ANGLE_POS_SLOPE");
     48    double GLINT_ANGLE_POS_REF = psMetadataLookupF32 (&status, recipe, "GLINT_ANGLE_POS_REF");
    4549    double glintWidth = psMetadataLookupF32 (&status, recipe, "GLINT_WIDTH");
    4650    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
    47    
    4851
    4952    // select the set of glint regions (GLINT.REGION is a MULTI of METADATA items)
    … …  
    8891    psastroChipBounds (fpa);
    8992
    90     // find the possible glint star:s for each glint star, convert the position to FPA coordinates.
     93
     94    // find the possible glint stars, and convert the position to FPA coordinates.
    9195    // search for stars within the glint regions
    9296    for (int i = 0; i < glintStars->n; i++) {
     97
    9398        pmAstromObj *star = glintStars->data[i];
    9499        if (star->Mag > GLINT_MAX_MAG) continue; // XXX should not be needed...
    … …  
    115120            psRegion glintRegion = psRegionFromString (glintRegionString);
    116121
    117             // select stars that land in this regio
     122            // select stars that land in this region
    118123            if (star->FP->x < glintRegion.x0) continue;
    119124            if (star->FP->x > glintRegion.x1) continue;
    120125            if (star->FP->y < glintRegion.y0) continue;
    121126            if (star->FP->y > glintRegion.y1) continue;
    122 
    123             double glintLength = GLINT_LENGTH_MAG_SLOPE*(GLINT_LENGTH_MAG_ZERO - star->Mag);
    124127
    125128            char *glintType = psMetadataLookupStr (&status, glintItem->data.md, "GLINT.TYPE");
    … …  
    129132            }
    130133
    131             // depending on the glint type, we need to find either the chips in the row or in the column.
     134            double glintLength = GLINT_LENGTH_MAG_SLOPE*(GLINT_LENGTH_MAG_ZERO - star->Mag);
     135            double glintAngle = 0;
     136
     137            //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     
     138            if ((!strcasecmp(glintType, "TOP") || !strcasecmp(glintType, "BOTTOM")) && abs(star->FP->y) < 30000 ){
     139              glintLength /= GLINT_LENGTH_POS_SLOPE*(GLINT_LENGTH_POS_REF - abs(star->FP->y));
     140
     141            }
     142            if ((!strcasecmp(glintType, "LEFT") || !strcasecmp(glintType, "RIGHT")) && abs(star->FP->x) < 30000 ) {
     143              glintLength /= GLINT_LENGTH_POS_SLOPE*(GLINT_LENGTH_POS_REF - abs(star->FP->x));
     144            }
     145            //also compute the angle of the glint, which depends on position parallel to the FPX
     146            if (!strcasecmp(glintType, "TOP") || !strcasecmp(glintType, "BOTTOM") ){
     147              glintAngle = PM_RAD_DEG * (GLINT_ANGLE_POS_SLOPE*((GLINT_ANGLE_POS_REF - star->FP->x)/1000.));
     148
     149            }
     150            if (!strcasecmp(glintType, "LEFT") || !strcasecmp(glintType, "RIGHT") ) {
     151              glintAngle = PM_RAD_DEG * (GLINT_ANGLE_POS_SLOPE*((GLINT_ANGLE_POS_REF - star->FP->y)/1000.));
     152            }
     153
     154
     155            psVector *x_glint = psVectorAlloc(2,PS_TYPE_F32);
     156            psVector *y_glint = psVectorAlloc(2,PS_TYPE_F32);
     157
    132158            if (!strcasecmp(glintType, "TOP") || !strcasecmp(glintType, "BOTTOM")) {
    133                 for (int nChip = 0; nChip < fpa->chips->n; nChip++) {
    134 
    135                     pmChip *chip = fpa->chips->data[nChip];
    136                     if (!chip) continue;
    137 
    138                     if (!psastroFindChipInXrange (fpa, nChip, star->FP->x, star->FP->y)) {
    139                         continue;
    140                     }
    141 
    142                     // find the coordinate of the end point
    143                     double xFPA  = star->FP->x;
    144                     double yFPA0 = star->FP->y;
    145                     double yFPA1;
    146                     if (!strcasecmp(glintType, "TOP")) {
    147                         yFPA1 = yFPA0 - glintLength;
    148                     } else {
    149                         yFPA1 = yFPA0 + glintLength;
    150                     }
    151 
    152                     // FPA coordinates of intersections with chip edges
    153                     double yFPAs, yFPAe;
    154                     psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip);
    155 
    156                     if (yFPA0 > yFPA1) PS_SWAP (yFPA0, yFPA1);
    157                     if (yFPAs > yFPAe) PS_SWAP (yFPAs, yFPAe);
    158 
    159                     // does this glint cross this chip?
    160                     if (yFPA0 > yFPAe) continue;
    161                     if (yFPA1 < yFPAs) continue;
    162 
    163                     yFPA0 = PS_MAX (yFPA0, yFPAs);
    164                     yFPA1 = PS_MIN (yFPA1, yFPAe);
    165 
    166                     double xChip0, yChip0, xChip1, yChip1;
    167                     psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, xFPA, yFPA0);
    168                     psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, xFPA, yFPA1);
    169 
    170                     // we now have the location of the glint ends or edge-intersections in chip coordinates
    171                     // double-check if this chip overlaps this glint
    172                     if (xChip0 > xChip1) PS_SWAP (xChip0, xChip1);
    173                     if (yChip0 > yChip1) PS_SWAP (yChip0, yChip1);
    174 
    175                     // bounds of this chip
    176                     psRegion *region = pmChipPixels (chip);
    177                     if (yChip1 < region->y0) continue;
    178                     if (yChip0 > region->y1) continue;
    179 
    180                     // this glint touches this chip. calculate the start and end
    181                     // coordinates on this chip
    182                     double xChipS = PS_MAX(xChip0 - 0.5*glintWidth, region->x0);
    183                     double xChipE = PS_MIN(xChip1 + 0.5*glintWidth, region->x1);
    184                     double yChipS = PS_MAX (yChip0, region->y0);
    185                     double yChipE = PS_MIN (yChip1, region->y1);
    186 
    187                     // select the 0th readout of the 0th cell for this chip
    188                     if (!chip->cells) continue;
    189                     if (!chip->cells->n) continue;
    190                     pmCell *glintCell = chip->cells->data[0];
    191                     if (!glintCell) continue;
    192                     if (!glintCell->readouts) continue;
    193                     if (!glintCell->readouts->n) continue;
    194                     pmReadout *glintReadout = glintCell->readouts->data[0];
    195                     if (!glintReadout) continue;
    196 
    197                     // save the glints on the readout->analysis metadata, creating if needed
    198                     psArray *glints = psMetadataLookupPtr (&status, glintReadout->analysis, "PSASTRO.GLINTS");
    199                     if (glints == NULL) {
    200                         glints = psArrayAllocEmpty (100);
    201                         if (!psMetadataAdd (glintReadout->analysis, PS_LIST_TAIL, "PSASTRO.GLINTS", PS_DATA_ARRAY, "astrometry matches", glints)) {
    202                             psWarning("failure to add glints to readout");
    203                             psFree (glints);
    204                             continue;
    205                         }
    206                         psFree (glints);
    207                     }
    208 
    209                     fprintf (stderr, "glint %s : %f,%f to %f,%f (%f,%f to %f,%f)\n", glintType, xChip0, yChip0, xChip1, yChip1, xChipS, yChipS, xChipE, yChipE);
    210                     psRegion *glint = psRegionAlloc(xChipS, xChipE, yChipS, yChipE);
    211                     psArrayAdd (glints, 100, glint);
    212                     psFree (glint);
    213                     psFree (region);
     159                // 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
     160                double xFPA0 = star->FP->x;
     161                double yFPA0 = star->FP->y;
     162                double xFPA1;
     163                double yFPA1;
     164                //angles for TOP and LEFT have been flipped in the fitting
     165                if (!strcasecmp(glintType, "TOP")) {
     166                    //find the edge chip and determine the FPA coords of the edge. Then, grab the edge coord and new glint length
     167                    for (int nChip = 0; nChip < fpa->chips->n; nChip++) {
     168                        pmChip *chip = fpa->chips->data[nChip];
     169                        if (!chip) continue;
     170
     171                        if (!psastroFindChipInXrange (fpa, nChip, star->FP->x, 20000.)) {
     172                          continue;
     173                        }
     174                        if (!psastroFindChipInYrange (fpa, nChip, star->FP->x, 20000.)) {
     175                          continue;
     176                        }
     177
     178                        // FPA coordinates of intersections with chip edges
     179                        double yFPAs, yFPAe;
     180                        psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip);
     181
     182                        xFPA0 = star->FP->x;
     183                        yFPA0 = yFPAe;
     184                        glintLength-= (star->FP->y-yFPAe);
     185                    }
     186
     187                    xFPA1 = xFPA0 + glintLength*sin(glintAngle*-1.);
     188                    yFPA1 = yFPA0 - glintLength*cos(glintAngle*-1.);
     189
     190                } else {
     191                    //find the edge chip and determine the FPA coords of the edge. Then, grab the edge coord and new glint length
     192                    for (int nChip = 0; nChip < fpa->chips->n; nChip++) {
     193                        pmChip *chip = fpa->chips->data[nChip];
     194                        if (!chip) continue;
     195
     196                        if (!psastroFindChipInXrange (fpa, nChip, star->FP->x, -20000.)) {
     197                          continue;
     198                        }
     199                        if (!psastroFindChipInYrange (fpa, nChip, star->FP->x, -20000.)) {
     200                          continue;
     201                        }
     202
     203                        // FPA coordinates of intersections with chip edges
     204                        double yFPAs, yFPAe;
     205                        psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip);
     206
     207                        xFPA0 = star->FP->x;
     208                        yFPA0 = yFPAs;
     209                        glintLength-= (yFPAs-star->FP->y);
     210                    }
     211
     212                    xFPA1 = xFPA0 - glintLength*sin(glintAngle);
     213                    yFPA1 = yFPA0 + glintLength*cos(glintAngle);
     214                }
     215
     216                x_glint->data.F32[0] = xFPA0;
     217                y_glint->data.F32[0] = yFPA0;
     218                x_glint->data.F32[1] = xFPA1;
     219                y_glint->data.F32[1] = yFPA1;
     220
     221                //we need to loop over each corner to select the chips that can have the glint on it
     222                for (int glint_point = 0; glint_point < 2; glint_point++) {
     223                    for (int nChip = 0; nChip < fpa->chips->n; nChip++) {
     224
     225                        pmChip *chip = fpa->chips->data[nChip];
     226                        if (!chip) continue;
     227
     228                        if (!psastroFindChipInXrange (fpa, nChip, x_glint->data.F32[glint_point], y_glint->data.F32[glint_point])) {
     229                          continue;
     230                        }
     231
     232                        // FPA coordinates of intersections with chip edges
     233                        double yFPAs, yFPAe;
     234                        psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip);
     235
     236                        if (yFPAs > yFPAe) PS_SWAP (yFPAs, yFPAe);
     237                        if (yFPA0 > yFPA1) {
     238                            PS_SWAP (xFPA0, xFPA1);
     239                            PS_SWAP (yFPA0, yFPA1);
     240                        }
     241
     242                        // does this glint cross this chip?
     243                        if (yFPA0 > yFPAe) continue;
     244                        if (yFPA1 < yFPAs) continue;
     245
     246                        //find the y-coord positions for this chip
     247                        double ycFPA0, ycFPA1;
     248                        ycFPA0 = PS_MAX (yFPA0, yFPAs);
     249                        ycFPA1 = PS_MIN (yFPA1, yFPAe);
     250
     251                        //now calculate the proper x-coord positions given the angle, for this chip
     252                        double xcFPA0, xcFPA1, angle;
     253                        angle = atan2(xFPA1 - xFPA0,yFPA1 - yFPA0);
     254                        xcFPA0 = xFPA0 + (ycFPA0-yFPA0)*tan(angle);
     255                        xcFPA1 = xFPA0 + (ycFPA1-yFPA0)*tan(angle);
     256                       
     257
     258                        double xChip0, yChip0, xChip1, yChip1, chip_angle, glint_length, xChip2,yChip2;
     259                        psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, xcFPA0, ycFPA0);
     260                        psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, xcFPA1, ycFPA1);
     261                        psastroFPAtoChip (&xChip2, &yChip2, fpa, nChip, star->FP->x,star->FP->y);
     262
     263                        chip_angle = atan2(yChip1 - yChip0, xChip1 - xChip0);
     264                        glint_length = sqrt(pow(yChip1 - yChip0,2) + pow(xChip1 - xChip0,2));
     265
     266                        // select the 0th readout of the 0th cell for this chip
     267                        if (!chip->cells) continue;
     268                        if (!chip->cells->n) continue;
     269                        pmCell *glintCell = chip->cells->data[0];
     270                        if (!glintCell) continue;
     271                        if (!glintCell->readouts) continue;
     272                        if (!glintCell->readouts->n) continue;
     273                        pmReadout *glintReadout = glintCell->readouts->data[0];
     274                        if (!glintReadout) continue;
     275                       
     276                        // save the glints on the readout->analysis metadata, creating if needed
     277                        psArray *glints = psMetadataLookupPtr (&status, glintReadout->analysis, "PSASTRO.GLINTS");
     278                        if (glints == NULL) {
     279                          glints = psArrayAllocEmpty (100);
     280                          if (!psMetadataAdd (glintReadout->analysis, PS_LIST_TAIL, "PSASTRO.GLINTS", PS_DATA_ARRAY, "astrometry matches", glints)) {
     281                            psWarning("failure to add glints to readout");
     282                            psFree (glints);
     283                            continue;
     284                        }
     285                          psFree (glints);
     286                        }
     287                       
     288                        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);
     289                        psVector *glint = psVectorAlloc(5,PS_TYPE_F32);
     290                        glint->data.F32[0] = xChip0;
     291                        glint->data.F32[1] = yChip0;
     292                        glint->data.F32[2] = glint_length;
     293                        glint->data.F32[3] = glintWidth;
     294                        glint->data.F32[4] = chip_angle;
     295
     296                        psArrayAdd (glints, 100, glint);
     297
     298                        psFree (glint);
     299                    }
    214300                }
    215301            }
    216302
    217             // depending on the glint type, we need to find either the chips in the row or in the column.
    218303            if (!strcasecmp(glintType, "LEFT") || !strcasecmp(glintType, "RIGHT")) {
    219                 for (int nChip = 0; nChip < fpa->chips->n; nChip++) {
    220 
    221                     pmChip *chip = fpa->chips->data[nChip];
    222                     if (!chip) continue;
    223 
    224                     if (!psastroFindChipInYrange (fpa, nChip, star->FP->x, star->FP->y)) {
    225                         continue;
    226                     }
    227 
    228                     // find the coordinate of the end point
    229                     double yFPA  = star->FP->y;
    230                     double xFPA0 = star->FP->x;
    231                     double xFPA1;
    232                     if (!strcasecmp(glintType, "RIGHT")) {
    233                         xFPA1 = xFPA0 - glintLength;
    234                     } else {
    235                         xFPA1 = xFPA0 + glintLength;
    236                     }
    237 
    238                     // FPA coordinates of intersections with chip edges
    239                     double xFPAs, xFPAe;
    240                     psastroFindChipXedges (&xFPAs, &xFPAe, fpa, nChip);
    241 
    242                     if (xFPA0 > xFPA1) PS_SWAP (xFPA0, xFPA1);
    243                     if (xFPAs > xFPAe) PS_SWAP (xFPAs, xFPAe);
    244 
    245                     // does this glint cross this chip?
    246                     if (xFPA0 > xFPAe) continue;
    247                     if (xFPA1 < xFPAs) continue;
    248 
    249                     xFPA0 = PS_MAX (xFPA0, xFPAs);
    250                     xFPA1 = PS_MIN (xFPA1, xFPAe);
    251 
    252                     double xChip0, yChip0, xChip1, yChip1;
    253                     psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, xFPA0, yFPA);
    254                     psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, xFPA1, yFPA);
    255 
    256                     // we now have the location of the glint ends or edge-intersections in chip coordinates
    257                     // double-check if this chip overlaps this glint
    258                     if (xChip0 > xChip1) PS_SWAP (xChip0, xChip1);
    259                     if (yChip0 > yChip1) PS_SWAP (yChip0, yChip1);
    260 
    261                     // bounds of this chip
    262                     psRegion *region = pmChipPixels (chip);
    263                     if (xChip1 < region->x0) continue;
    264                     if (xChip0 > region->x1) continue;
    265 
    266                     // this glint touches this chip. calculate the start and end
    267                     // coordinates on this chip
    268                     double yChipS = PS_MAX (yChip0 - 0.5*glintWidth, region->y0);
    269                     double yChipE = PS_MIN (yChip1 + 0.5*glintWidth, region->y1);
    270                     double xChipS = PS_MAX (xChip0, region->x0);
    271                     double xChipE = PS_MIN (xChip1, region->x1);
    272 
    273                     // select the 0th readout of the 0th cell for this chip
    274                     if (!chip->cells) continue;
    275                     if (!chip->cells->n) continue;
    276                     pmCell *glintCell = chip->cells->data[0];
    277                     if (!glintCell) continue;
    278                     if (!glintCell->readouts) continue;
    279                     if (!glintCell->readouts->n) continue;
    280                     pmReadout *glintReadout = glintCell->readouts->data[0];
    281                     if (!glintReadout) continue;
    282 
    283                     // save the glints on the readout->analysis metadata, creating if needed
    284                     psArray *glints = psMetadataLookupPtr (&status, glintReadout->analysis, "PSASTRO.GLINTS");
    285                     if (glints == NULL) {
    286                         glints = psArrayAllocEmpty (100);
    287                         if (!psMetadataAdd (glintReadout->analysis, PS_LIST_TAIL, "PSASTRO.GLINTS", PS_DATA_ARRAY, "astrometry matches", glints)) {
    288                             psWarning("failure to add glints to readout");
    289                             psFree (glints);
    290                             continue;
    291                         }
    292                         psFree (glints);
    293                     }
    294 
    295                     fprintf (stderr, "glint %s : %f,%f to %f,%f (%f,%f to %f,%f)\n", glintType, xChip0, yChip0, xChip1, yChip1, xChipS, yChipS, xChipE, yChipE);
    296                     psRegion *glint = psRegionAlloc(xChipS, xChipE, yChipS, yChipE);
    297                     psArrayAdd (glints, 100, glint);
    298                     psFree (glint);
    299                     psFree (region);
     304                // 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
     305                double xFPA0 = star->FP->x;
     306                double yFPA0 = star->FP->y;
     307                double xFPA1;
     308                double yFPA1;
     309                //angles for TOP and LEFT have been flipped in the fitting
     310                if (!strcasecmp(glintType, "RIGHT")) {
     311                    //find the edge chip and determine the FPA coords of the edge. Then, grab the edge coord and new glint length
     312                    for (int nChip = 0; nChip < fpa->chips->n; nChip++) {
     313                        pmChip *chip = fpa->chips->data[nChip];
     314                        if (!chip) continue;
     315
     316                        if (!psastroFindChipInXrange (fpa, nChip, 20000.,star->FP->y)) {
     317                          continue;
     318                        }
     319                        if (!psastroFindChipInYrange (fpa, nChip, 20000.,star->FP->y)) {
     320                          continue;
     321                        }
     322
     323                        // FPA coordinates of intersections with chip edges
     324                        double yFPAs, yFPAe;
     325                        psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip);
     326
     327                        xFPA0 = yFPAe;
     328                        yFPA0 = star->FP->y;
     329                        glintLength-= (star->FP->x-yFPAe);
     330                    }
     331
     332                    xFPA1 = xFPA0 - glintLength*cos(glintAngle);
     333                    yFPA1 = yFPA0 - glintLength*sin(glintAngle);
     334
     335                } else {
     336                    //find the edge chip and determine the FPA coords of the edge. Then, grab the edge coord and new glint length
     337                    for (int nChip = 0; nChip < fpa->chips->n; nChip++) {
     338                        pmChip *chip = fpa->chips->data[nChip];
     339                        if (!chip) continue;
     340
     341                        if (!psastroFindChipInXrange (fpa, nChip, -20000.,star->FP->y)) {
     342                          continue;
     343                        }
     344                        if (!psastroFindChipInYrange (fpa, nChip, -20000.,star->FP->y)) {
     345                          continue;
     346                        }
     347
     348                        // FPA coordinates of intersections with chip edges
     349                        double yFPAs, yFPAe;
     350                        psastroFindChipYedges (&yFPAs, &yFPAe, fpa, nChip);
     351
     352                        xFPA0 = yFPAs;
     353                        yFPA0 = star->FP->y;
     354                        glintLength-= (yFPAs-star->FP->x);
     355                    }
     356
     357                    xFPA1 = xFPA0 + glintLength*cos(glintAngle*-1.);
     358                    yFPA1 = yFPA0 + glintLength*sin(glintAngle*-1.);
     359                }
     360
     361                x_glint->data.F32[0] = xFPA0;
     362                y_glint->data.F32[0] = yFPA0;
     363                x_glint->data.F32[1] = xFPA1;
     364                y_glint->data.F32[1] = yFPA1;
     365
     366                //we need to loop over each corner to select the chips that can have the glint on it
     367                for (int glint_point = 0; glint_point < 2; glint_point++) {
     368                    for (int nChip = 0; nChip < fpa->chips->n; nChip++) {
     369
     370                        pmChip *chip = fpa->chips->data[nChip];
     371                        if (!chip) continue;
     372
     373                        if (!psastroFindChipInYrange (fpa, nChip, x_glint->data.F32[glint_point], y_glint->data.F32[glint_point])) {
     374                          continue;
     375                        }
     376
     377                        // FPA coordinates of intersections with chip edges
     378                        double xFPAs, xFPAe;
     379                        psastroFindChipXedges (&xFPAs, &xFPAe, fpa, nChip);
     380
     381                        if (xFPAs > xFPAe) PS_SWAP (xFPAs, xFPAe);
     382                        if (xFPA0 > xFPA1) {
     383                            PS_SWAP (xFPA0, xFPA1);
     384                            PS_SWAP (yFPA0, yFPA1);
     385                        }
     386
     387                        // does this glint cross this chip?
     388                        if (xFPA0 > xFPAe) continue;
     389                        if (xFPA1 < xFPAs) continue;
     390
     391                        //find the x-coord positions for this chip
     392                        double xcFPA0, xcFPA1;
     393                        xcFPA0 = PS_MAX (xFPA0, xFPAs);
     394                        xcFPA1 = PS_MIN (xFPA1, xFPAe);
     395
     396                        //now calculate the proper x-coord positions given the angle, for this chip
     397                        double ycFPA0, ycFPA1, angle;
     398                        angle = atan2(xFPA1 - xFPA0,yFPA1 - yFPA0);
     399                        ycFPA0 = yFPA0 + (xcFPA0-xFPA0)/tan(angle);
     400                        ycFPA1 = yFPA0 + (xcFPA1-xFPA0)/tan(angle);
     401
     402                        double xChip0, yChip0, xChip1, yChip1, chip_angle, glint_length,xChip2,yChip2;
     403                        psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, xcFPA0, ycFPA0);
     404                        psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, xcFPA1, ycFPA1);
     405                        psastroFPAtoChip (&xChip2, &yChip2, fpa, nChip, star->FP->x,star->FP->y);
     406
     407                        chip_angle = atan2(yChip1 - yChip0, xChip1 - xChip0);
     408                        glint_length = sqrt(pow(yChip1 - yChip0,2) + pow(xChip1 - xChip0,2));
     409
     410                        // select the 0th readout of the 0th cell for this chip
     411                        if (!chip->cells) continue;
     412                        if (!chip->cells->n) continue;
     413                        pmCell *glintCell = chip->cells->data[0];
     414                        if (!glintCell) continue;
     415                        if (!glintCell->readouts) continue;
     416                        if (!glintCell->readouts->n) continue;
     417                        pmReadout *glintReadout = glintCell->readouts->data[0];
     418                        if (!glintReadout) continue;
     419                       
     420                        // save the glints on the readout->analysis metadata, creating if needed
     421                        psArray *glints = psMetadataLookupPtr (&status, glintReadout->analysis, "PSASTRO.GLINTS");
     422                        if (glints == NULL) {
     423                          glints = psArrayAllocEmpty (100);
     424                          if (!psMetadataAdd (glintReadout->analysis, PS_LIST_TAIL, "PSASTRO.GLINTS", PS_DATA_ARRAY, "astrometry matches", glints)) {
     425                            psWarning("failure to add glints to readout");
     426                            psFree (glints);
     427                            continue;
     428                        }
     429                          psFree (glints);
     430                        }
     431                       
     432                        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);
     433                        psVector *glint = psVectorAlloc(5,PS_TYPE_F32);
     434                        glint->data.F32[0] = xChip0;
     435                        glint->data.F32[1] = yChip0;
     436                        glint->data.F32[2] = glint_length;
     437                        glint->data.F32[3] = glintWidth;
     438                        glint->data.F32[4] = chip_angle;
     439
     440                        psArrayAdd (glints, 100, glint);
     441
     442                        psFree (glint);
     443                    }
    300444                }
    301445            }
     446
     447
    302448            if (!strcasecmp(glintType, "HSC")) {
    303449              // It's inefficient to keep looking these up.

  • trunk/psastro/src/psastroMaskUpdates.c

    r40490 r41338  
    134134   
    135135    double REFSTAR_MASK_MAX_MAG            = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_MAX_MAG");
    136     double REFSTAR_MASK_SATSTAR_MAG_SLOPE  = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_MAG_SLOPE");
    137136    double REFSTAR_MASK_SATSTAR_MAG_MAX    = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_MAG_MAX");
     137    double REFSTAR_MASK_SATSTAR_RAD_OFFSET = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_RAD_OFFSET");
     138    double REFSTAR_MASK_SATSTAR_EXP        = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_EXP");
    138139    double REFSTAR_MASK_SATSTAR_POS_ZERO   = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_POS_ZERO");
    139140    double REFSTAR_MASK_SATSPIKE_MAG_SLOPE = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSPIKE_MAG_SLOPE");
    … …  
    232233    double ROTANGLE = PM_RAD_DEG * psMetadataLookupF64 (&status, fpa->concepts, "FPA.ROTANGLE");
    233234    psAssert (status, "ROTANGLE missing");
     235    int ROT_PARITY = psMetadataLookupS32(&status, recipe, "PSASTRO.MODEL.ROT.PARITY");
     236    if (!status) psAbort ("Can't find recipe option PSASTRO.MODEL.ROT.PARITY");
     237
     238    //select the astromeric rotation angle from the header
     239    psMetadata *header = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
     240    if (!header) {
     241        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO header");
     242        return false;
     243    }
     244    float AST_T0 = PM_RAD_DEG * psMetadataLookupF32(&status, header, "AST_T0");
    234245
    235246    // de-activate all files except PSASTRO.INPUT.MASK and PSASTRO.OUTPUT.MASK
    … …  
    327338
    328339                    // CIRCLE around the stars (scaled by magnitude)
    329                     float radius = REFSTAR_MASK_SATSTAR_MAG_SLOPE * (REFSTAR_MASK_SATSTAR_MAG_MAX - ref->Mag);
     340                    float radius = REFSTAR_MASK_SATSTAR_RAD_OFFSET + pow(REFSTAR_MASK_SATSTAR_EXP,(REFSTAR_MASK_SATSTAR_MAG_MAX - ref->Mag));
    330341
    331342                    // XXX for now, assume cell binning is 1x1 relative to chip
    … …  
    333344
    334345                    for (float theta = 0.0; theta < 2*M_PI; theta += M_PI / 2.0) {
    335 
    336                         float Theta = theta - ROTANGLE - REFSTAR_MASK_SATSTAR_POS_ZERO;
     346                        float Theta = theta - (ROT_PARITY * ROTANGLE) - AST_T0 - REFSTAR_MASK_SATSTAR_POS_ZERO;
    337347
    338348                        // LINE for boundaries of the saturation spikes (scaled by magnitude)
    339                         // float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2;
     349                        //float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2;
    340350
    341351                        float MAG_MAX = REFSTAR_MASK_SATSPIKE_MAG_MAX;
    … …  
    400410                        }
    401411                    }
     412
    402413                }
    403414
    … …  
    412423                    for (int i = 0; i < ghosts->n; i++) {
    413424                        psastroGhost *ghost = ghosts->data[i];
     425                        //offset the ghost angles with the astrometric rotation angle
     426                        ghost->inner.theta -= (ROT_PARITY * ROTANGLE)
     427                        ghost->outer.theta -= (ROT_PARITY * ROTANGLE)
    414428                        // XXX bright vs faint ghost bits? (OR with SUSPECT)
    415429                        psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer);
    … …  
    418432
    419433                // Select the glint mask regions for this readout (loaded in
    420                 // psastroExtractGlints.c).  These glint regions are defined as rectangular
     434                // psastroChooseGlintStars.c).  These glint regions are defined as rectangular
    421435                // boxes and are generated for each chip based on the position of the bright
    422436                // stars beyond the edge of the focal plane.  This masking is currently very
    … …  
    426440                    // mask the glints on this readout
    427441                    for (int i = 0; i < glints->n; i++) {
    428                         psRegion *glint = glints->data[i];
    429                         psastroMaskRectangle (readoutMask->mask, glintMaskValue, glint->x0, glint->y0, glint->x1, glint->y1);
     442                        psVector *glint = glints->data[i];
     443                        //psLogMsg ("psastro", 3, "glint: %f %f %f %f\n", glint->x0, glint->y0, glint->x1, glint->y1);
     444                        //psastroMaskRectangle (readoutMask->mask, glintMaskValue, glint->x0, glint->y0, glint->x1, glint->y1);
     445
     446                        psastroMaskBox (readoutMask->mask, glintMaskValue, glint->data.F32[0], glint->data.F32[1],
     447                                    glint->data.F32[2], glint->data.F32[3], glint->data.F32[4]);
     448
    430449                    }
    431450                }
Note: See TracChangeset for help on using the changeset viewer.