GEANT321/ggeom/gneltu.F

   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