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 GNOELT(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 outside 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),TAU(2)
37 C
38 #if !defined(CERNLIB_SINGLE)
39       DOUBLE PRECISION SAFZ,A2,B2,X0,Y0,PHI1,PHI2,PHI3,X3,Y3
40       DOUBLE PRECISION DXY2,U,V,W,DISCR,SQDISC,TAU,TAUZ,ZI,XZ,YZ
41 #endif
42       SNXT = BIG
43       SAFZ = ABS(X(3))-PAR(3)
44       A2 = PAR(1)*PAR(1)
45       B2 = PAR(2)*PAR(2)
46 C
47       IF(IACT.EQ.3)GOTO 40
48 C
49 C      -----------------------------------
50 C      |  Compute safety-distance 'SAFE' |
51 C      -----------------------------------
52 C
53 C ....  First check Z
54       X0 = ABS(X(1))
55       Y0 = ABS(X(2))
56 C
57       SAFE=0.
58       IF(X0*X0/A2+Y0*Y0/B2.LT.1.) GO TO 30
59       PHI1=0.
60       PHI2=PIBY2
61       DO 10    I=1,10
62          PHI3=(PHI1+PHI2)*0.5
63          X3=PAR(1)*COS(PHI3)
64          Y3=PAR(2)*SIN(PHI3)
65          D=Y3*A2*(X0-X3)-X3*B2*(Y0-Y3)
66 *
67 *        D is the (signed) distance from the point (X0,Y0)
68 *        to the normal to the ellipse at the point (X3,Y3).
69 *
70          IF (D.LT.0.) THEN
71             PHI1=PHI3
72          ELSE
73             PHI2=PHI3
74          END IF
75    10 CONTINUE
76    20 SAFE=SQRT((X0-X3)**2+(Y0-Y3)**2)-.01
77    30 IF(SAFZ.GT.0.)THEN
78 *
79 * ....   Combine the radial distance whit the Z-distance
80 *
81          SAFE = SQRT(SAFE**2+SAFZ**2)
82       ENDIF
83       IF(IACT.EQ.0)GOTO 999
84       IF(IACT.EQ.1.AND.SNEXT.LT.SAFE)GOTO 999
85    40 CONTINUE
86 C
87 C      ---------------------------------------
88 C      |  Compute the vector-distance 'SNXT' |
89 C      ---------------------------------------
90 C
91       IF(SAFZ.GT.0.0.AND.X(3)*X(6).GE.0.0)GOTO 999
92 C
93       DXY2 = (1-X(6))*(1+X(6))
94       IF(DXY2.LE.0.)GOTO 60
95 C
96 C ....  Find the intersection of the given ray
97 C       (described by array X) whit the cylider.
98 C       Ray equation : X'(1-3) = X(1-3) + TAU*X(4-6)
99 C       Cylinder equation : x**2/a**2 + y**2/b**2 = 1
100 C       To obtain TAU,solve the quadratic equation
101 C       Ut**2 + 2Vt + W = 0
102 C
103       U = X(4)*X(4)*B2+X(5)*X(5)*A2
104       V = X(1)*X(4)*B2+X(2)*X(5)*A2
105       W = X(1)*X(1)*B2+X(2)*X(2)*A2-A2*B2
106       DISCR = V*V-U*W
107       IF(DISCR.LT.0.)GOTO 999
108       IF(U.EQ.0.)GOTO 999
109       SQDISC = SQRT(DISCR)
110       TAU(1) = (-V+SQDISC)/U
111       TAU(2) = (-V-SQDISC)/U
112 C
113       DO 50 I=1,2
114          IF(TAU(I).GE.0.)THEN
115             ZI = X(3)+TAU(I)*X(6)
116             IF((ABS(ZI)-PAR(3)).LT.1.E-6)THEN
117 C
118 C ....  Set SNXT to the smallest positive TAU,only if
119 C       the intersection point is inside the Z limits
120 C
121                SNXT = MIN(SNXT,REAL(TAU(I)))
122             ENDIF
123          ENDIF
124    50 CONTINUE
125    60 CONTINUE
126 C
127       IF(SAFZ.GT.0.)THEN
128 C
129 C ....  Check intersection whit Z planes
130 C
131          IF(X(3).GT.0.) ZI = PAR(3)
132          IF(X(3).LT.0.) ZI = -PAR(3)
133 C
134          TAUZ = (ZI-X(3))/X(6)
135          XZ = X(1)+X(4)*TAUZ
136          YZ = X(2)+X(5)*TAUZ
137          IF((XZ*XZ/A2+YZ*YZ/B2).LE.1.) SNXT = TAUZ
138       ENDIF
139 C
140   999 END