GEANT321/ggeom/gnpcon.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:53  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 GNPCON (X, PAR, IACT, SNEXT, SNXT, SAFE)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'PCON' 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       DIMENSION X(6), PAR(9),  XT(6), PT(7)
  36
  37       EQUIVALENCE (PT(1), DZ),  (PT(2), PT2), (PT(3), PT3)
  38       EQUIVALENCE (PT(4), PT4), (PT(5), PT5), (PT(6), PT6)
  39       EQUIVALENCE (PT(7), PT7)
  40
  41       EQUIVALENCE (XT(3), XT3)
  42 C.
  43 C.   ---------------------------------------------------
  44 C.
  45
  46       SNXT = BIG
  47       R2  = X(1)*X(1) + X(2)*X(2)
  48       R   = SQRT (R2)
  49
  50       NZ  = PAR(3)
  51       ZMIN  = PAR(4)
  52       ZMAX  = PAR(3*NZ+1)
  53       SAFZ1 = X(3) - ZMIN
  54       SAFZ2 = ZMAX - X(3)
  55       SAFEZ = MIN (SAFZ1, SAFZ2)
  56
  57 C
  58 C......  First determine in which z-segment the particle is located.
  59 C
  60       DO 10 JPH=7, 3*(NZ-1)+1,  3
  61         IF (X(3) .LE. PAR(JPH)) THEN
  62            IPH=JPH
  63            GO TO 20
  64         ENDIF
  65    10 CONTINUE
  66       IPH=3*NZ+1
  67
  68    20 CONTINUE
  69 C
  70 C......  The particle is in the segment bounded by z-planes at
  71 C......  Z1=PAR(IPL) and Z2=PAR(IPH), i.e.,  Z1 < X(3) < Z2.
  72 C
  73 C......  Set parameters for this segment and translate z-coordinate
  74 C......  of point relative to center of this segment.This is done in
  75 C......  preparation of invoking the algorithms used in "GNTUBE" and
  76 C......  "GNCONE" (which for reasons of efficiency and clarity are
  77 C......  implemented inline).
  78 C
  79       IPL  = IPH - 3
  80       DZ    = 0.5 * (PAR(IPH) - PAR(IPL))
  81       PT2 = PAR(IPL+1)
  82       PT3 = PAR(IPL+2)
  83       PT4 = PAR(IPH+1)
  84       PT5 = PAR(IPH+2)
  85       PT6 = PAR(1)
  86       PT7 = PAR(1) + PAR(2)
  87       IF (PT7 .GT. 360.0)  PT7  = PT7 - 360.0
  88
  89       XT3 = X(3) - 0.5 * (PAR(IPL) + PAR(IPH))
  90
  91       XT(1)  = X(1)
  92       XT(2)  = X(2)
  93       XT(4)  = X(4)
  94       XT(5)  = X(5)
  95       XT(6)  = X(6)
  96
  97       IND = 2
  98       IF (PAR(2) .EQ. 360.0)  IND = 1
  99
 100       IF (IACT .LT. 3) THEN
 101
 102 C       -------------------------------------------------
 103 C       |  Compute safety-distance 'SAFE' (P.Weidhaas)  |
 104 C       -------------------------------------------------
 105         SAFZ1   = DZ - ABS(XT3)
 106         SAFEZ   = MIN (SAFEZ,SAFZ1)
 107         SAFSEG  = BIG
 108 C
 109 C......  Determine whether the segment is a tube or a cone.
 110 C
 111
 112         IF (PT2 .NE. PT4) GO TO 50
 113         IF (PT3 .NE. PT5) GO TO 50
 114
 115 C*********************************************************
 116 C
 117 C......  The segment is a tube: invoke the algorithm
 118 C......  from routine "GNTUBE" inline to get "SAFER".
 119 C
 120 C*********************************************************
 121
 122         SAFR1  = R - PT2
 123         SAFR2  = PT3 - R
 124         SAFER  = MIN (SAFR1, SAFR2)
 125
 126         IF (IND .EQ. 2) GO TO 70
 127
 128         GO TO 100
 129
 130
 131    50   CONTINUE
 132
 133 C*********************************************************
 134 C
 135 C......  The segment is a cone: invoke the algorithm
 136 C......  from routine "GNCONE" inline to get "SAFER".
 137 C
 138 C*********************************************************
 139
 140         SAFZ2  = DZ + XT3
 141         ZLENI  = 0.5 / DZ
 142
 143 C......  Compute radial distance to inner wall.
 144
 145         FACT  = (PT4 - PT2) * ZLENI
 146         R1  = PT2 + FACT * SAFZ2
 147         SAFR1 = (R - R1) / SQRT(1.0 + FACT*FACT)
 148
 149 C......  Compute radial distance to outer wall.
 150
 151         FACT  = (PT5 - PT3) * ZLENI
 152         R2  = PT3 + FACT * SAFZ2
 153         SAFR2 = (R2 - R) / SQRT(1.0 + FACT*FACT)
 154         SAFER = MIN (SAFR1, SAFR2)
 155
 156         IF (IND .EQ. 1) GO TO 100
 157
 158    70   CONTINUE
 159
 160 C********************************************************************
 161 C......  Here we handle the case of a PHI-segment of a tube or cone.
 162 C......  In addition to the radial distances (SAFR1, SAFR2) and the
 163 C......  axial distances (SAFZ1, SAFZ2) we compute here the distance
 164 C......  to the PHI-segment boundary that is closest to the point.
 165 C
 166 C......  For each PHI-boundary we find the distance from the given
 167 C......  point to the outer (at R2) point of the segment boundary
 168 C......  (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
 169 C......  "SAFSEG" to be the distance to segment PHI1, else we set
 170 C......  "SAFSEG" to be the distance to segment PHI2.
 171 C*********************************************************************
 172
 173           PHI1  = PT6 * DEGRAD
 174           PHI2  = PT7 * DEGRAD
 175           IF (PHI2 .LT. PHI1)  PHI2 = PHI2 + TWOPI
 176
 177           COSPH1  = COS (PHI1)
 178           COSPH2  = COS (PHI2)
 179           SINPH1  = SIN (PHI1)
 180           SINPH2  = SIN (PHI2)
 181
 182 C......  Get coordinates of outer endpoints (at R2) of both PHI-segments.
 183
 184           XS1  = R2 * COSPH1
 185           YS1  = R2 * SINPH1
 186           XS2  = R2 * COSPH2
 187           YS2  = R2 * SINPH2
 188
 189 C......  Get distances (squared) from given point to each endpoint.
 190
 191           DISTS1 = (X(1) - XS1)**2  +  (X(2) - YS1)**2
 192           DISTS2 = (X(1) - XS2)**2  +  (X(2) - YS2)**2
 193
 194 C......  Get distance to that PHI-segment whose endpoint
 195 C......  is closest to the given point.
 196
 197           IF (DISTS1 .LE. DISTS2) THEN
 198             SAFSEG = ABS(X(1) * SINPH1 - X(2) * COSPH1)
 199           ELSE
 200             SAFSEG = ABS(X(1) * SINPH2 - X(2) * COSPH2)
 201           ENDIF
 202
 203
 204   100   CONTINUE
 205
 206         SAFE  = MIN (SAFEZ, SAFER, SAFSEG)
 207
 208         IF (IACT .EQ. 0) GO TO 999
 209         IF (IACT .EQ. 1) THEN
 210           IF (SNEXT .LT. SAFE) GO TO 999
 211         ENDIF
 212       ENDIF
 213
 214
 215 C     ------------------------------------------------
 216 C     |  Compute vector-distance 'SNXT' (McPherson)  |
 217 C     ------------------------------------------------
 218 *
 219 *     Avoid rounding effects induced by translation ********************
 220 *
 221       IF (ABS(XT(3)).GE.DZ) XT(3) = (1.-0.000001)*XT(3)
 222 *
 223       IF (PT2 .NE. PT4) GO TO 200
 224       IF (PT3 .NE. PT5) GO TO 200
 225
 226       DELZ  = DZ
 227       PT(1) = PT2
 228       PT(2) = PT3
 229       PT(3) = DELZ
 230       PT(4) = PT6
 231       PT(5) = PT7
 232
 233       CALL GNTUBE (XT, PT, 3, IND, SNEXT, SNXT, TSAFE)
 234       GO TO 999
 235
 236   200 CONTINUE
 237       CALL GNCONE (XT, PT, 3, IND, SNEXT, SNXT, TSAFE)
 238
 239   999 CONTINUE
 240       END
 241