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 GINTCO (X, RLEFT, RRIGHT, DZ, TAU, TAUL)
14 C        ********************************************
15 C        *  This subroutine finds the intersection  *
16 C        *  of a given ray (described by array X)   *
17 C        *  with a given cone (described by radii   *
18 C        *  RLEFT and RRIGHT and half-distance DZ). *
19 C        *  Output parameter is TAU, and inter-     *
20 C        *  section point is X = XP + TAU*XD, etc.  *
21 C        *                                          *
22 C        *  Called by GNOCON                        *
23 C        *  Programmed by:  Patrick Weidhaas        *
24 C        *          CERN,   March 1988              *
25 C        ********************************************
27 #include "geant321/gconsp.inc"
28       DIMENSION X(6)
29 #if !defined(CERNLIB_SINGLE)
30       DOUBLE PRECISION XP,YP,ZP,XD,YD,ZD,S,T,U,V,W
31       DOUBLE PRECISION DISCR,SQDISC
32 #endif
33 C----------------------------------------------------
35 C...... Point of origin of ray:
37       XP = X(1)
38       YP = X(2)
39       ZP = X(3)
41 C...... Direction cosines:
43       XD = X(4)
44       YD = X(5)
45       ZD = X(6)
47       TAU  = BIG
48       TAUL = BIG
50       S = 0.5 * (RLEFT + RRIGHT)
51       T = (RLEFT - RRIGHT) / DZ
53 C......  Cone equation is:    x**2 + y**2 - Az**2 + Bz + C = 0
55       A = 0.25 * T*T
56       B = S * T
57       C = -S*S
59 C......  To obtain "TAU", we must solve the quadratic equation
60 C......  Ut**2 + Vt + W = 0 .
62       U = XD**2 + YD**2 - A*ZD**2
63       V = 2.0 * (XP*XD + YP*YD - A*ZP*ZD) + B*ZD
64       W = XP**2 + YP**2 - A*ZP**2 + B*ZP + C
66       DISCR = V*V - 4.0*U*W
67       IF (DISCR .LE. 0.0) GO TO 999
68       IF(U.EQ.0.)GO TO 999
69       SQDISC = SQRT (DISCR)
70       TAU1 = (-V + SQDISC) / (2.0*U)
71       TAU2 = (-V - SQDISC) / (2.0*U)
74 C......  Set TAU to the smallest positive root;
75 C......  otherwise let TAU = BIG .
76 C
77 C......  If both roots are positive, set TAUL to
78 C......  the larger one: it may be needed in the
79 C......  case of a PHI-segmented cone.
81       IF (TAU1 .LT. 0.0) THEN
82         IF (TAU2 .LT. 0.0) GO TO 999
83         TAU = TAU2
84       ELSE
85         TAU = TAU1
86         IF (TAU2 .GT. 0.0) THEN
87           TAUL = TAU2
88           IF (TAU2.LT.TAU1) THEN
89             TAU = TAU2
90             TAUL = TAU1
91           ENDIF
92         ENDIF
93       ENDIF
95   999 CONTINUE
96       END