[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gnoelt.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:20:52  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.29  by  S.Giani
*-- Author :
      SUBROUTINE GNOELT(X,PAR,IACT,SNEXT,SNXT,SAFE)
C
C     ****************************************************************
C     *                                                              *
C     *     Compute distance up to intersection with 'ELTU' volume,  *
C     *     from outside point X(1-3) along direction X(4-6).        *
C     *                                                              *
C     *     PAR    (input)  : volume parameters                      *
C     *     IACT   (input)  : action flag                            *
C     *       = 0   Compute SAFE only                                *
C     *       = 1   Compute SAFE, and SNXT only if SNEXT.gt.new SAFE *
C     *       = 2   Compute both SAFE and SNXT                       *
C     *       = 3   Compute SNXT only                                *
C     *     SNEXT  (input)  : see IACT = 1                           *
C     *     SNXT   (output) : distance to volume boundary            *
C     *     SAFE   (output) : shortest distance to any boundary      *
C     *                                                              *
C     *  ==>Called by : GNEXT,GTNEXT                                 *
C     *       Author  A.Solano                                       *
C     *                                                              *
C     ****************************************************************
C
#include "geant321/gconsp.inc"
C
      DIMENSION X(6),PAR(3),TAU(2)
C
#if !defined(CERNLIB_SINGLE)
      DOUBLE PRECISION SAFZ,A2,B2,X0,Y0,PHI1,PHI2,PHI3,X3,Y3
      DOUBLE PRECISION DXY2,U,V,W,DISCR,SQDISC,TAU,TAUZ,ZI,XZ,YZ
#endif
      SNXT = BIG
      SAFZ = ABS(X(3))-PAR(3)
      A2 = PAR(1)*PAR(1)
      B2 = PAR(2)*PAR(2)
C
      IF(IACT.EQ.3)GOTO 40
C
C      -----------------------------------
C      |  Compute safety-distance 'SAFE' |
C      -----------------------------------
C
C ....  First check Z
      X0 = ABS(X(1))
      Y0 = ABS(X(2))
C
      SAFE=0.
      IF(X0*X0/A2+Y0*Y0/B2.LT.1.) GO TO 30
      PHI1=0.
      PHI2=PIBY2
      DO 10    I=1,10
         PHI3=(PHI1+PHI2)*0.5
         X3=PAR(1)*COS(PHI3)
         Y3=PAR(2)*SIN(PHI3)
         D=Y3*A2*(X0-X3)-X3*B2*(Y0-Y3)
*
*        D is the (signed) distance from the point (X0,Y0)
*        to the normal to the ellipse at the point (X3,Y3).
*
         IF (D.LT.0.) THEN
            PHI1=PHI3
         ELSE
            PHI2=PHI3
         END IF
   10 CONTINUE
   20 SAFE=SQRT((X0-X3)**2+(Y0-Y3)**2)-.01
   30 IF(SAFZ.GT.0.)THEN
*
* ....   Combine the radial distance whit the Z-distance
*
         SAFE = SQRT(SAFE**2+SAFZ**2)
      ENDIF
      IF(IACT.EQ.0)GOTO 999
      IF(IACT.EQ.1.AND.SNEXT.LT.SAFE)GOTO 999
   40 CONTINUE
C
C      ---------------------------------------
C      |  Compute the vector-distance 'SNXT' |
C      ---------------------------------------
C
      IF(SAFZ.GT.0.0.AND.X(3)*X(6).GE.0.0)GOTO 999
C
      DXY2 = (1-X(6))*(1+X(6))
      IF(DXY2.LE.0.)GOTO 60
C
C ....  Find the intersection of the given ray
C       (described by array X) whit the cylider.
C       Ray equation : X'(1-3) = X(1-3) + TAU*X(4-6)
C       Cylinder equation : x**2/a**2 + y**2/b**2 = 1
C       To obtain TAU,solve the quadratic equation
C       Ut**2 + 2Vt + W = 0
C
      U = X(4)*X(4)*B2+X(5)*X(5)*A2
      V = X(1)*X(4)*B2+X(2)*X(5)*A2
      W = X(1)*X(1)*B2+X(2)*X(2)*A2-A2*B2
      DISCR = V*V-U*W
      IF(DISCR.LT.0.)GOTO 999
      IF(U.EQ.0.)GOTO 999
      SQDISC = SQRT(DISCR)
      TAU(1) = (-V+SQDISC)/U
      TAU(2) = (-V-SQDISC)/U
C
      DO 50 I=1,2
         IF(TAU(I).GE.0.)THEN
            ZI = X(3)+TAU(I)*X(6)
            IF((ABS(ZI)-PAR(3)).LT.1.E-6)THEN
C
C ....  Set SNXT to the smallest positive TAU,only if
C       the intersection point is inside the Z limits
C
               SNXT = MIN(SNXT,REAL(TAU(I)))
            ENDIF
         ENDIF
   50 CONTINUE
   60 CONTINUE
C
      IF(SAFZ.GT.0.)THEN
C
C ....  Check intersection whit Z planes
C
         IF(X(3).GT.0.) ZI = PAR(3)
         IF(X(3).LT.0.) ZI = -PAR(3)
C
         TAUZ = (ZI-X(3))/X(6)
         XZ = X(1)+X(4)*TAUZ
         YZ = X(2)+X(5)*TAUZ
         IF((XZ*XZ/A2+YZ*YZ/B2).LE.1.) SNXT = TAUZ
      ENDIF
C
  999 END
Commit	Line	Data
fe4da5cc	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(X0X0/A2+Y0Y0/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=Y3A2(X0-X3)-X3B2(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(SAFE2+SAFZ2)
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 : x2/a2 + y2/b2 = 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 = VV-UW
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((XZXZ/A2+YZYZ/B2).LE.1.) SNXT = TAUZ
138	ENDIF
139	C
140	999 END