1 *
2 * \$Id\$
3 *
4 * \$Log\$
5 * Revision 1.1.1.1  1995/10/24 10:20:52  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 GNELTU(X,PAR,IACT,SNEXT,SNXT,SAFE)
13 C
14 C     ****************************************************************
15 C     *                                                              *
16 C     *     Compute distance up to intersection with 'ELTU' 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.Solano                                       *
31 C     *                                                              *
32 C     ****************************************************************
33 C
34 #include "geant321/gconsp.inc"
35 C
36       DIMENSION X(6),PAR(3)
37 #if !defined(CERNLIB_SINGLE)
38       DOUBLE PRECISION SAFZ1,SAFZ2,SAFR,A2,B2,X0,Y0,A,B,D1,D2
39       DOUBLE PRECISION X1,X2,X3,Y1,Y2,Y3
40       DOUBLE PRECISION SZ,XZ,YZ,U,V,W,DISCR,SQDISC,TAU1,TAU2
41 #endif
42 C
43       SNXT = BIG
44       SAFZ1 = PAR(3)-X(3)
45       SAFZ2 = PAR(3)+X(3)
46 C
47       A2 = PAR(1)*PAR(1)
48       B2 = PAR(2)*PAR(2)
49 C
50       IF(IACT.LT.3)THEN
51 C
52 C        -----------------------------------
53 C        |  Compute safety-distance 'SAFE' |
54 C        -----------------------------------
55 C
56          X0 = ABS(X(1))
57          Y0 = ABS(X(2))
58 C
59          A=PAR(1)
60          B=PAR(2)
61          X1=X0
62          Y1=SQRT((A-X0)*(A+X0))*B/A
63          Y2=Y0
64          X2=SQRT((B-Y0)*(B+Y0))*A/B
65          D1=(X1-X0)**2+(Y1-Y0)**2
66          D2=(X2-X0)**2+(Y2-Y0)**2
67          DO 1 I=1,8
68             IF (B.LT.A) THEN
69                X3=.5*(X1+X2)
70                Y3=SQRT((A-X3)*(A+X3))*B/A
71             ELSE
72                Y3=.5*(Y1+Y2)
73                X3=SQRT((B-Y3)*(B+Y3))*A/B
74             END IF
75             IF (D1.LT.D2) THEN
76                X2=X3
77                Y2=Y3
78                D2=(X2-X0)**2+(Y2-Y0)**2
79             ELSE
80                X1=X3
81                Y1=Y3
82                D1=(X1-X0)**2+(Y1-Y0)**2
83             END IF
84     1    CONTINUE
85     2    SAFR=SQRT(D1)-1.E-3
86 *
87          SAFE = MIN(SAFZ1,SAFZ2,SAFR)
88          IF(IACT.EQ.0)GOTO 99
89          IF(IACT.EQ.1.AND.SNEXT.LT.SAFE)GOTO 99
90 C
91       ENDIF
92 C
93 C        -----------------------------------
94 C        |  Compute vector-distance 'SNXT' |
95 C        -----------------------------------
96 C
97 C ....  First check Z
98 C
99       IF(X(6).GT.0.)THEN
100          SNXT = SAFZ1/X(6)
101       ELSEIF(X(6).LT.0.)THEN
102          SNXT = -SAFZ2/X(6)
103       ENDIF
104 C
105 C ....  Then,if necessary,find the intersection of
106 C       the given ray(described by array X) whit
107 C       the cylinder.
108 C       Ray equation : X'(1-3) = X(1-3) + TAU*X(4-6)
109 C       Cylinder equation : x**2/a**2 + y**2/b**2 = 1
110 C       To obtain TAU,solve the quadratic equation
111 C       Ut**2 + 2Vt + W = 0
112 C
113       SZ = SNXT
114       XZ = X(1)+X(4)*SZ
115       YZ = X(2)+X(5)*SZ
116       IF((XZ*XZ/A2+YZ*YZ/B2).LE.1)GOTO 99
117 C
118       U = X(4)*X(4)*B2+X(5)*X(5)*A2
119       V = X(1)*X(4)*B2+X(2)*X(5)*A2
120       W = X(1)*X(1)*B2+X(2)*X(2)*A2-A2*B2
121       DISCR = V*V-U*W
122       IF(DISCR.LT.0.)GOTO 99
123       IF(U.EQ.0.)GOTO 99
124       SQDISC = SQRT(DISCR)
125       TAU1 = (-V+SQDISC)/U
126       TAU2 = (-V-SQDISC)/U
127 C
128 C ....  Set SNXT to the smallest positive TAU
129 C
130       IF(TAU1.LT.0.)THEN
131          IF(TAU2.LT.0.)GOTO 99
132          SNXT = TAU2
133       ELSE
134          SNXT = TAU1
135          IF(TAU2.GT.0.0.AND.TAU2.LT.TAU1)THEN
136             SNXT = TAU2
137          ENDIF
138       ENDIF
139 C
140    99 END