GEANT321/ggeom/gnpgon.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:54  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.30  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GNPGON (X, PAR, IACT, SNEXT, SNXT, SAFE)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'PGON' VOLUME,  *
  17 C.    *       FROM INSIDE POINT X(1-3) ALONG DIRECTION X(4-6)          *
  18 C.    *                                                                *
  19 C.    *       PAR   (input)  : volume parameters                       *
  20 C.    *       IACT  (input)  : action flag                             *
  21 C.    *         = 0  Compute SAFE only                                 *
  22 C.    *         = 1  Compute SAFE, and SNXT only if SNEXT .GT.new SAFE *
  23 C.    *         = 2  Compute both SAFE and SNXT                        *
  24 C.    *         = 3  Compute SNXT only                                 *
  25 C.    *       SNEXT (input)  : see IACT = 1                            *
  26 C.    *       SNXT  (output) : distance to volume boundary             *
  27 C.    *       SAFE  (output) : shortest distance to any boundary       *
  28 C.    *                                                                *
  29 C.    *    ==>Called by : GNEXT, GTNEXT                                *
  30 C.    *         Author  A.McPherson,  P.Weidhaas  *********            *
  31 C.    *                                                                *
  32 C.    ******************************************************************
  33 C.
  34 #include "geant321/gconsp.inc"
  35
  36       DIMENSION X(6), PAR(50), TPAR(50)
  37
  38 C-----------------------------------------------------------------
  39
  40       SNXT = BIG
  41       SAFSEG = BIG
  42       R2 = X(1)*X(1) + X(2)*X(2)
  43       R  = SQRT (R2)
  44
  45       PDIV  = PAR(2) / PAR(3)
  46       DELPHI  = PDIV * DEGRAD
  47       DPHI2 = 0.5 * DELPHI
  48       CSDPH2 = COS (DPHI2)
  49
  50       TPAR(1)  = PAR(1)
  51       TPAR(2)  = PAR(2)
  52       TPAR(3)  = PAR(3)
  53       TPAR(4)  = PAR(4)
  54       NZ  = PAR(4)
  55
  56       DO 5 I=1, NZ
  57         I3  = 3*(I-1)
  58         TPAR(5+I3) = PAR(5+I3)
  59         TPAR(6+I3) = PAR(6+I3) / CSDPH2
  60         TPAR(7+I3) = PAR(7+I3)
  61     5 CONTINUE
  62
  63 C**************************************************************************
  64 C
  65 C......  Here we start the logic from "GNPCON" (which for reasons of
  66 C......  efficiency and clarity has been implemented inline).
  67 C
  68 C**************************************************************************
  69
  70       ZMIN  = TPAR(5)
  71       ZMAX  = TPAR(3*NZ+2)
  72       SAFZ1 = X(3) - ZMIN
  73       SAFZ2 = ZMAX - X(3)
  74       SAFEZ = MIN (SAFZ1, SAFZ2)
  75
  76 C
  77 C......  First determine in which segment the particle is located.
  78 C
  79       DO 10 JPH=8, 3*NZ-1,  3
  80         IF (X(3) .LT. TPAR(JPH)) THEN
  81            IPH=JPH
  82            GO TO 20
  83         ENDIF
  84    10 CONTINUE
  85       IPH = 3*NZ+2
  86
  87    20 CONTINUE
  88 C
  89 C......  The particle is in the segment bounded by z-planes at
  90 C......  Z1=PAR(IPL) and Z2=PAR(IPH), i.e.,  Z1 < X(3) < Z2.
  91 C
  92 C......  Set parameters for this segment and translate z-coordinate
  93 C......  of point relative to center of this segment. this is done in
  94 C......  preparation of invoking the algorithms used in "GNTUBE" and
  95 C......  "GNCONE" (which for reasons of efficiency and clarity are
  96 C......  implemented inline).
  97 C
  98       IPL  = IPH - 3
  99       DZ    = 0.5 * (TPAR(IPH) - TPAR(IPL))
 100       PT2 = TPAR(IPL+1)
 101       PT3 = TPAR(IPL+2)
 102       PT4 = TPAR(IPH+1)
 103       PT5 = TPAR(IPH+2)
 104       PT6 = TPAR(1)
 105       PT7 = TPAR(1) + TPAR(2)
 106       IF (PT7 .GT. 360.0)  PT7  = PT7 - 360.0
 107
 108       XT3 = X(3) - 0.5 * (TPAR(IPL) + TPAR(IPH))
 109
 110       SAFZ2  = DZ + XT3
 111       ZLENI  = 0.5 / DZ
 112
 113       PHI1  = PT6 * DEGRAD
 114       PHI2  = PT7 * DEGRAD
 115       IF (PHI2.LE.PHI1) PHI2=PHI2+TWOPI
 116
 117       IF (IACT .LT. 3) THEN
 118
 119 C       -------------------------------------------------
 120 C       |  Compute safety-distance 'SAFE' (P.Weidhaas)  |
 121 C       -------------------------------------------------
 122
 123         IFLAG = 2
 124         IF (TPAR(2) .EQ. 360.0)  IFLAG = 1
 125         SAFZ1 = DZ - ABS(XT3)
 126         SAFEZ = MIN (SAFEZ,SAFZ1)
 127 C
 128 C......  Next determine whether the segment is a tube or a cone.
 129 C
 130         IF (PT2 .NE. PT4) GO TO 50
 131         IF (PT3 .NE. PT5) GO TO 50
 132
 133 C*********************************************************
 134 C
 135 C......  The segment is a tube; invoke the algorithm
 136 C......  from routine "GNTUBE" inline to get "SAFER".
 137 C
 138 C*********************************************************
 139
 140         SAFR1  = BIG
 141         IF(PT2.GT.0.) SAFR1  = R - PT2
 142         SAFR2  = PT3 - R
 143         SAFER  = MIN (SAFR1, SAFR2)
 144
 145         IF (IFLAG .EQ. 2) GO TO 70
 146
 147         GO TO 100
 148
 149
 150    50   CONTINUE
 151
 152 C*********************************************************
 153 C
 154 C......  The segment is a cone; invoke the algorithm
 155 C......  from routine "GNCONE" inline to get "SAFER".
 156 C
 157 C*********************************************************
 158
 159
 160 C......  Compute radial distance to inner wall.
 161
 162         IF(PT2+PT4.GT.0.) THEN
 163         FACT  = (PT4 - PT2) * ZLENI
 164         RAD1  = PT2 + FACT * SAFZ2
 165         SAFR1 = (R - RAD1) / SQRT(1.0 + FACT*FACT)
 166         ELSE
 167         SAFR1 = BIG
 168         ENDIF
 169
 170 C......  Compute radial distance to outer wall.
 171
 172         FACT  = (PT5 - PT3) * ZLENI
 173         RAD2  = PT3 + FACT * SAFZ2
 174         SAFR2 = (RAD2 - R) / SQRT(1.0 + FACT*FACT)
 175
 176         SAFER  = MIN (SAFR1, SAFR2)
 177
 178         IF (IFLAG .EQ. 1) GO TO 100
 179
 180    70   CONTINUE
 181
 182 C********************************************************************
 183 C......  Here we handle the case of a phi-segment of a tube or cone.
 184 C......  in addition to the radial distances (SAFR1, SAFR2) and the
 185 C......  axial distances (SAFZ1, SAFZ2) we compute here the distance
 186 C......  to the phi-segment boundary that is closest to the point:
 187 C
 188 C......  For each phi-boundary we find the distance from the given
 189 C......  point to the outer (at R2) point of the segment boundary
 190 C......  (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
 191 C......  "SAFSEG" to be the distance to segment PHI1, else we set
 192 C......  "SAFSEG" to be the distance to segment PHI2.
 193 C*********************************************************************
 194
 195
 196           COSPH1  = COS (PHI1)
 197           SINPH1  = SIN (PHI1)
 198           COSPH2  = COS (PHI2)
 199           SINPH2  = SIN (PHI2)
 200
 201 C......  Get coordinates of outer endpoints (at R2) of both phi-segments.
 202
 203           XS1  = R * COSPH1
 204           YS1  = R * SINPH1
 205           XS2  = R * COSPH2
 206           YS2  = R * SINPH2
 207
 208 C......  Get distances (squared) from given point to each endpoint.
 209
 210           DISTS1 = (X(1) - XS1)**2  +  (X(2) - YS1)**2
 211           DISTS2 = (X(1) - XS2)**2  +  (X(2) - YS2)**2
 212
 213 C......  Get distance to that phi-segment whose endpoint
 214 C......  is closest to the given point.
 215
 216           IF (DISTS1 .LE. DISTS2) THEN
 217             SAFSEG = ABS(SINPH1 * X(1) - COSPH1 * X(2))
 218           ELSE
 219             SAFSEG = ABS(SINPH2 * X(1) - COSPH2 * X(2))
 220           ENDIF
 221
 222   100   CONTINUE
 223
 224
 225         IF (SAFER .LE. 0.0)  THEN
 226
 227 C---------------------------------------------------------------------------
 228 C
 229 C......  Here we handle the case in which  SAFER < 0, i.e., the point is
 230 C......  inside the polygon but outside the inscribed polycone. We must
 231 C......  do an accurate calculation of "SAFER".
 232 C
 233 C---------------------------------------------------------------------------
 234
 235           FACT = SAFZ2 * ZLENI
 236           RAD1 = PT2 + FACT * (PT4 - PT2)
 237           RAD2 = PT3 + FACT * (PT5 - PT3)
 238           RR1 = RAD1 * RAD1
 239           RR2 = RAD2 * RAD2
 240
 241           IF(X(1).EQ.0.)THEN
 242              PHI   = ATAN2(X(2), 1.E-8)
 243           ELSE
 244              PHI   = ATAN2(X(2), X(1))
 245           ENDIF
 246           IF (PHI .LT. PHI1)  PHI = PHI + TWOPI
 247
 248           DIST=0.
 249           IF (PHI .GE. PHI1  .AND.  PHI .LE. PHI2)  THEN
 250             PHIREL = PHI - PHI1
 251             NSECTR = INT(PHIREL / DELPHI) + 1
 252             PHICTR  = PHI1 + (2.0*NSECTR - 1.0) * DPHI2
 253             COSPHC = COS (PHICTR)
 254             SINPHC = SIN (PHICTR)
 255
 256
 257             IF (R2 .GE. RR2) THEN
 258               SR2 = RAD2
 259               DIST = ABS (COSPHC * X(1) + SINPHC * X(2) - SR2)
 260               FACTC  = (PT5 - PT3) * ZLENI
 261             ELSEIF (R2 .LE. RR1) THEN
 262               RIN  = RAD1 * CSDPH2
 263               SRIN = RIN
 264               DIST = ABS (COSPHC * X(1) + SINPHC * X(2) - SRIN)
 265               FACTC = (PT4 - PT2) * ZLENI
 266             ENDIF
 267
 268           ENDIF
 269
 270           SAFER  = DIST / SQRT(1.0 + FACTC*FACTC)
 271
 272         ENDIF
 273
 274         SAFE  = MIN (SAFEZ, SAFER, SAFSEG)
 275
 276
 277         IF (IACT .EQ. 0) GO TO 999
 278         IF (IACT .EQ. 1) THEN
 279           IF (SNEXT .LT. SAFE) GO TO 999
 280         ENDIF
 281       ENDIF
 282
 283
 284 C     ------------------------------------------------
 285 C     |  Compute vector-distance 'SNXT' (Nierhaus)  |
 286 C     ------------------------------------------------
 287
 288
 289       CALL GNPGO1(X,PAR,SNXT)
 290 C
 291   999 CONTINUE
 292       END
 293