This commit was generated by cvs2svn to compensate for changes in r2,

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

diff --git a/GEANT321/ggeom/gneltu.F b/GEANT321/ggeom/gneltu.F
new file mode 100644 (file)
index 0000000..094e056
--- /dev/null
+++ b/GEANT321/ggeom/gneltu.F
@@ -0,0 +1,140 @@
+*
+* $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 GNELTU(X,PAR,IACT,SNEXT,SNXT,SAFE)
+C
+C     ****************************************************************
+C     *                                                              *
+C     *     Compute distance up to intersection with 'ELTU' volume,  *
+C     *     from inside 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)
+#if !defined(CERNLIB_SINGLE)
+      DOUBLE PRECISION SAFZ1,SAFZ2,SAFR,A2,B2,X0,Y0,A,B,D1,D2
+      DOUBLE PRECISION X1,X2,X3,Y1,Y2,Y3
+      DOUBLE PRECISION SZ,XZ,YZ,U,V,W,DISCR,SQDISC,TAU1,TAU2
+#endif
+C
+      SNXT = BIG
+      SAFZ1 = PAR(3)-X(3)
+      SAFZ2 = PAR(3)+X(3)
+C
+      A2 = PAR(1)*PAR(1)
+      B2 = PAR(2)*PAR(2)
+C
+      IF(IACT.LT.3)THEN
+C
+C        -----------------------------------
+C        |  Compute safety-distance 'SAFE' |
+C        -----------------------------------
+C
+         X0 = ABS(X(1))
+         Y0 = ABS(X(2))
+C
+         A=PAR(1)
+         B=PAR(2)
+         X1=X0
+         Y1=SQRT((A-X0)*(A+X0))*B/A
+         Y2=Y0
+         X2=SQRT((B-Y0)*(B+Y0))*A/B
+         D1=(X1-X0)**2+(Y1-Y0)**2
+         D2=(X2-X0)**2+(Y2-Y0)**2
+         DO 1 I=1,8
+            IF (B.LT.A) THEN
+               X3=.5*(X1+X2)
+               Y3=SQRT((A-X3)*(A+X3))*B/A
+            ELSE
+               Y3=.5*(Y1+Y2)
+               X3=SQRT((B-Y3)*(B+Y3))*A/B
+            END IF
+            IF (D1.LT.D2) THEN
+               X2=X3
+               Y2=Y3
+               D2=(X2-X0)**2+(Y2-Y0)**2
+            ELSE
+               X1=X3
+               Y1=Y3
+               D1=(X1-X0)**2+(Y1-Y0)**2
+            END IF
+    1    CONTINUE
+    2    SAFR=SQRT(D1)-1.E-3
+*
+         SAFE = MIN(SAFZ1,SAFZ2,SAFR)
+         IF(IACT.EQ.0)GOTO 99
+         IF(IACT.EQ.1.AND.SNEXT.LT.SAFE)GOTO 99
+C
+      ENDIF
+C
+C        -----------------------------------
+C        |  Compute vector-distance 'SNXT' |
+C        -----------------------------------
+C
+C ....  First check Z
+C
+      IF(X(6).GT.0.)THEN
+         SNXT = SAFZ1/X(6)
+      ELSEIF(X(6).LT.0.)THEN
+         SNXT = -SAFZ2/X(6)
+      ENDIF
+C
+C ....  Then,if necessary,find the intersection of
+C       the given ray(described by array X) whit
+C       the cylinder.
+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
+      SZ = SNXT
+      XZ = X(1)+X(4)*SZ
+      YZ = X(2)+X(5)*SZ
+      IF((XZ*XZ/A2+YZ*YZ/B2).LE.1)GOTO 99
+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 99
+      IF(U.EQ.0.)GOTO 99
+      SQDISC = SQRT(DISCR)
+      TAU1 = (-V+SQDISC)/U
+      TAU2 = (-V-SQDISC)/U
+C
+C ....  Set SNXT to the smallest positive TAU
+C
+      IF(TAU1.LT.0.)THEN
+         IF(TAU2.LT.0.)GOTO 99
+         SNXT = TAU2
+      ELSE
+         SNXT = TAU1
+         IF(TAU2.GT.0.0.AND.TAU2.LT.TAU1)THEN
+            SNXT = TAU2
+         ENDIF
+      ENDIF
+C
+   99 END