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)
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