GEANT321/ggeom/gnpgon.F

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