GEANT321/ggeom/gnoelt.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 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