GEANT321/ggeom/gnctub.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:51  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.29  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GNCTUB (X, PAR, IACT, SNEXT, SNXT, SAFE)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *      COMPUTE DISTANCE UP TO INTERSECTION WITH 'CTUB'           *
  17 C.    *      VOLUME 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, GNPCON, GTNEXT                        *
  30 C.    *         Authors A.McPherson    ********                        *
  31 C.    *       MODIFICATION LOG.                                        *
  32 C.    *       18-July-89 M.Guckes modifications due to GEANG 3.13      *
  33 C.    *                                                                *
  34 C.    ******************************************************************
  35 C.
  36 #include "geant321/gconsp.inc"
  37       REAL X(6),PAR(11)
  38 C
  39 C-------------------------------------------------------------
  40 C
  41       SNXT = BIG
  42       R2   = X(1)*X(1) + X(2)*X(2)
  43       R    = SQRT (R2)
  44       IF(PAR(5).GE.360.) THEN
  45          IFLAG = 1
  46       ELSE
  47          IFLAG = 2
  48       ENDIF
  49 *
  50       SAFZ1  = (-PAR(3)-X(3) )*PAR(8)-X(1)*PAR(6)-X(2)*PAR(7)
  51       SAFZ2  = (PAR(3)-X(3) )*PAR(11)-X(1)*PAR(9)-X(2)*PAR(10)
  52       IF (PAR(1).NE.0.) THEN
  53          SAFR1  = R - PAR(1)
  54       ELSE
  55          SAFR1  = BIG
  56       ENDIF
  57       SAFR2  = PAR(2) - R
  58 *
  59       IF (IACT .LT. 3) THEN
  60 *
  61 * *** Compute safety-distance 'SAFE' (P.Weidhaas)
  62 *
  63          SAFSEG = BIG
  64          IF (IFLAG .EQ. 2) THEN
  65 *
  66 *     In addition to the radial distances (SAFR1 and SAFR2) and the
  67 *     axial distances (SAFZ1 and SAFZ2) we compute here the distance
  68 *     to the PHI-segment boundary that is closest to the point:
  69 *
  70 *     For each PHI-boundary we find the distance from the given
  71 *     point to the outer (at RMAX) point of the segment boundary
  72 *     (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
  73 *     SAFSEG to be the distance to segment PHI1. Else we set
  74 *     SAFSEG to be the distance to segment PHI2.
  75 *
  76             PHI1 = PAR(4) * DEGRAD
  77             PHI2 = PAR(5) * DEGRAD
  78 *
  79             COSPH1 = COS (PHI1)
  80             COSPH2 = COS (PHI2)
  81             SINPH1 = SIN (PHI1)
  82             SINPH2 = SIN (PHI2)
  83 *
  84 * ***   Get coordinates of outer endpoints (at ROUT) of both segments.
  85 *
  86             XS1 = PAR(2) * COSPH1
  87             YS1 = PAR(2) * SINPH1
  88             XS2 = PAR(2) * COSPH2
  89             YS2 = PAR(2) * SINPH2
  90 *
  91 * ***   Get distances (squared) from the given point to each endpoint.
  92 *
  93             DISTS1 = (X(1) - XS1)**2 + (X(2) - YS1)**2
  94             DISTS2 = (X(1) - XS2)**2 + (X(2) - YS2)**2
  95 *
  96 * ***   Get distance to that PHI-segment whose endpoint
  97 * ***   is closest to the given point.
  98 *
  99             IF (DISTS1 .LE. DISTS2) THEN
 100                SAFSEG = ABS(X(1) * SINPH1 - X(2) * COSPH1)
 101             ELSE
 102                SAFSEG = ABS(X(1) * SINPH2 - X(2) * COSPH2)
 103             ENDIF
 104          ENDIF
 105 *
 106          SAFE = MIN (SAFZ1, SAFZ2, SAFR1, SAFR2, SAFSEG)
 107          IF (IACT .EQ. 0) GO TO 999
 108          IF (IACT .EQ. 1 .AND. SNEXT .LT. SAFE) GO TO 999
 109       ENDIF
 110 *
 111 * ***   Compute intersection with z-boundaries
 112 *
 113       V1 = X(4)*PAR(6)+X(5)*PAR(7)+X(6)*PAR(8)
 114       SZ1 = SAFZ1/V1
 115       V2 = X(4)*PAR(9)+X(5)*PAR(10)+X(6)*PAR(11)
 116       SZ2 = SAFZ2/V2
 117       IF( SZ1 .GT. 0. ) THEN
 118          SNXT = SZ1
 119       ELSE
 120          SNXT = BIG
 121       ENDIF
 122       IF( SZ2 .GT. 0.0 .AND. SZ2 .LT. SNXT) SNXT = SZ2
 123       SZ=SNXT
 124 *
 125       IF (ABS(X(6)).LT.1.)THEN
 126 *
 127 * ***      Compute z-intercept with inner cylinder.
 128 *
 129          BA2=-1.
 130          IF(PAR(1).GT.0.)THEN
 131             RMIN2  = PAR(1)*PAR(1)
 132             ZP2    = 1./(1.-X(6)*X(6))
 133             BA     = (X(4)*X(1)+X(5)*X(2))*ZP2
 134             BA2    = BA*BA
 135             CA     = (R2-RMIN2)*ZP2
 136             DIS2   = BA2-CA
 137             IF (DIS2.GE.0.)THEN
 138                XSIN    = -BA-SQRT(DIS2)
 139                IF (XSIN.GE.0.)THEN
 140                   IF(XSIN.LT.SNXT)SNXT = XSIN
 141                   GO TO 10
 142                ENDIF
 143             ENDIF
 144          ENDIF
 145 *
 146 *  ***    Compute z-intercept with outer cylinder.
 147 *
 148          RMAX2  = PAR(2)*PAR(2)
 149          XZ     = X(1) + X(4)*SZ
 150          YZ     = X(2) + X(5)*SZ
 151          IF (XZ*XZ+YZ*YZ.GT.RMAX2)THEN
 152             IF(BA2.LT.0.)THEN
 153                ZP2    = 1./(1.-X(6)*X(6))
 154                BA     = (X(4)*X(1)+X(5)*X(2))*ZP2
 155                BA2    = BA*BA
 156             ENDIF
 157             CA     = (R2-RMAX2)*ZP2
 158             DIS2   = BA2-CA
 159             IF (DIS2.GE.0.)THEN
 160                SRMAX = -BA+SQRT(DIS2)
 161                IF(SRMAX.LT.SNXT)SNXT=SRMAX
 162             ENDIF
 163          ENDIF
 164       ENDIF
 165 *
 166    10 IF(IFLAG.EQ.2) THEN
 167 *
 168 *     =======>PHI segmented tube
 169 *             We have checked the radius and Z.
 170 *             now check PHI.
 171 *
 172          DPSGN=X(1)*X(5)-X(2)*X(4)
 173 *
 174 *             Tells us which way its going.
 175 *
 176          PHI2=PAR(5)
 177          IF(DPSGN.LT.0.0) PHI2=PAR(4)
 178          PH2=PHI2*DEGRAD
 179 *
 180 *             Have set up the limit.
 181 *
 182          TSGN=SIN(PH2)
 183          TCSG=COS(PH2)
 184          DX45=X(4)*TSGN-X(5)*TCSG
 185          IF(DX45.EQ.0.)GO TO 999
 186          SN1=(X(2)*TCSG-X(1)*TSGN)/DX45
 187 *
 188 *             Distance until tangents are right.
 189 *
 190          IF(SN1.LT.0.0) GO TO 999
 191          IF(ABS(TSGN).GT.1.E-6.AND.(X(2)+SN1*X(5))*TSGN .LT.0.)GO TO
 192      +   999
 193          IF(ABS(TCSG).GT.1.E-6.AND.(X(1)+SN1*X(4))*TCSG .LT.0.)GO TO
 194      +   999
 195 *
 196 *             Have checked that the distance is +VE and that the
 197 *             SINE is the right sign.
 198 *
 199          IF(SN1.LT.SNXT) SNXT=SN1
 200       END IF
 201 *
 202   999 END