GEANT321/ggeom/gnotub.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:53  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.30  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GNOTUB(X,P,IACT,IFL,SNEXT,SNXT,SAFE)
  13 C.    ******************************************************************
  14 C.    *                                                                *
  15 C.    *      Compute distance to intersection with boundary surface of *
  16 C     *      volume TUBE or TUBS, from point X(1),X(2),X(3) outside    *
  17 C     *      the volume along track with direction cosines X(4),X(5),  *
  18 C     *      X(6)                                                      *
  19 C.    *      P     (input)  : volume parameters                        *
  20 C.    *      IACT  (input)  : action flag                              *
  21 C.    *         = 0  Compute SAFE only                                 *
  22 C.    *         = 1  Compute SAFE, compute SNXT only if SAFE.LT.SNEXT  *
  23 C.    *         = 2  Compute both SAFE and SNXT                        *
  24 C.    *         = 3  Compute SNXT only                                 *
  25 C.    *      IFL   (input)  : 1 for TUBE, 2 for PHI segmented TUBE     *
  26 C.    *      SNEXT (input)  : see IACT = 1                             *
  27 C.    *      SNXT  (output) : distance to volume boundary along track  *
  28 C.    *      SAFE  (output) : not larger than scalar distance to       *
  29 C.    *                       volume boundaray                         *
  30 C.    *      Called by : GNEXT, GNOPCO, GTNEXT                         *
  31 C.    *                                                                *
  32 C.    *      Authors   : Michel Maire and Rolf Nierhaus   21-JUN-1990  *
  33 C.    *                                                                *
  34 C.    ******************************************************************
  35 C.    *                                                                *
  36 C.    * 'TUBE'    is a tube. It has 3 parameters, the inside radius,   *
  37 C.    *           the outside radius and the half length in z.         *
  38 C.    * 'TUBS'    is a phi segment of a  tube.  It has 5 parameters,   *
  39 C.    *           the same  3 as  'TUBE' plus  the phi  limits.  The   *
  40 C.    *           segment  starts at  the first  limit and  includes   *
  41 C.    *           increasing phi  value up  to the  second limit  or   *
  42 C.    *           that plus 360 degrees.                               *
  43 C.    *                                                                *
  44 C.    ******************************************************************
  45 #if !defined(CERNLIB_SINGLE)
  46       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  47 #endif
  48 #include "geant321/gconsp.inc"
  49       REAL X(6),P(5),SNEXT,SNXT,SAFE
  50 *
  51 *     this part has to be moved outside the routine
  52       IF (IFL.EQ.2) THEN
  53          P4=P(4)*DEGRAD
  54          P5=P(5)*DEGRAD
  55          IF (P5.LT.P4) P5=P5+TWOPI
  56          C1=COS(P4)
  57          S1=SIN(P4)
  58          C2=COS(P5)
  59          S2=SIN(P5)
  60          FIO=0.5*(P5+P4)
  61          CFIO=COS(FIO)
  62          SFIO=SIN(FIO)
  63          DFI=0.5*(P5-P4)
  64          CDFI=COS(DFI)
  65 *        SDFI=SIN(DFI)
  66       END IF
  67 *
  68       SNXT=1.E10
  69       R=SQRT(X(1)**2+X(2)**2)
  70 *
  71 *     Compute SAFE radius
  72       IF (IACT.LT.3) THEN
  73          SAF1=P(1)-R
  74          SAF2=R-P(2)
  75          SAF3=ABS(X(3))-P(3)
  76          SAF4=0.
  77          IF (IFL.EQ.2.AND.R.GT.0.) THEN
  78             CPSI=(X(1)*CFIO+X(2)*SFIO)/R
  79             IF (CPSI.LT.CDFI) THEN
  80                IF ((X(2)*CFIO-X(1)*SFIO).LE.0.) THEN
  81                   SAF4=ABS(X(1)*S1-X(2)*C1)
  82                ELSE
  83                   SAF4=ABS(X(1)*S2-X(2)*C2)
  84                END IF
  85             END IF
  86          ENDIF
  87          SAFE=MAX(SAF1,SAF2,SAF3,SAF4)
  88          IF (IACT.EQ.0) GO TO 999
  89          IF (IACT.EQ.1.AND.SNEXT.LE.SAFE) GO TO 999
  90       END IF
  91 *
  92 *     Intersection with z-plane
  93 *     only points outside the z range need to be considered
  94       IF (ABS(X(3)).GE.P(3)) THEN
  95          IF (X(3)*X(6).LT.0.) THEN
  96             S=(ABS(X(3))-P(3))/ABS(X(6))
  97             XI=X(1)+S*X(4)
  98             YI=X(2)+S*X(5)
  99             R2=XI**2+YI**2
 100             IF (P(1)**2.LE.R2.AND.R2.LE.P(2)**2) THEN
 101                IF (IFL.EQ.1) GO TO 101
 102                IF (R2.EQ.0.) GO TO 101
 103                CPSI=(XI*CFIO+YI*SFIO)/SQRT(R2)
 104                IF (CPSI.GE.CDFI) GO TO 101
 105             END IF
 106          END IF
 107       END IF
 108 *
 109 *     Intersection with cylinders
 110 *     Intersection point (x,y,z)
 111 *     (x,y,z) is on track :    x=X(1)+t*X(4)
 112 *                              y=X(2)+t*X(5)
 113 *                              z=X(3)+t*X(6)
 114 *     (x,y,z) is on cylinder : x**2 + y**2 = R**2
 115 *
 116 *     (X(4)**2+X(5)**2)*t**2
 117 *     +2.*(X(1)*X(4)+X(2)*X(5))*t
 118 *     +X(1)**2+X(2)**2-R**2=0
 119 *
 120       T1=X(4)**2+X(5)**2
 121       T2=(X(1)*X(4)+X(2)*X(5))
 122       T3=X(1)**2+X(2)**2
 123 *     track parallel to the z axis ?
 124 ***** 21-JUN-1990
 125 *     IF (T1.EQ.0.) GO TO 999
 126       IF (ABS(T1).LT.1.E-32) GO TO 999
 127       B =T2/T1
 128 *
 129 *     Intersection with outer cylinder
 130 *     only points outside the outer cylinder need to be considered
 131       IF (R.GE.P(2)) THEN
 132          C=(T3-P(2)**2)/T1
 133          D=B**2-C
 134          IF (D.GE.0.) THEN
 135             S=-B-SQRT(D)
 136             IF (S.GE.0.) THEN
 137                ZI=X(3)+S*X(6)
 138                IF (ABS(ZI).LE.P(3)) THEN
 139                   IF (IFL.EQ.1) GO TO 101
 140                   XI=X(1)+S*X(4)
 141                   YI=X(2)+S*X(5)
 142                   CPSI=(XI*CFIO+YI*SFIO)/P(2)
 143                   IF (CPSI.GE.CDFI) GO TO 101
 144                END IF
 145             END IF
 146          END IF
 147       END IF
 148 *     Intersection with inner cylinder
 149       IF (P(1).GT.0.) THEN
 150          C=(T3-P(1)**2)/T1
 151          D=B**2-C
 152          IF (D.GE.0.) THEN
 153             S=-B+SQRT(D)
 154             IF (S.GE.0.) THEN
 155                ZI=X(3)+S*X(6)
 156                IF (ABS(ZI).LE.P(3)) THEN
 157                   IF (IFL.EQ.1) GO TO 101
 158                   XI=X(1)+S*X(4)
 159                   YI=X(2)+S*X(5)
 160                   CPSI=(XI*CFIO+YI*SFIO)/P(1)
 161                   IF (CPSI.GE.CDFI) SNXT=S
 162                END IF
 163             END IF
 164          END IF
 165       END IF
 166 *
 167 *     Intersection with phi-planes
 168 *     x=r*cos(phi)=X(1)+t*X(4)
 169 *     y=r*sin(phi)=X(2)+t*X(5)
 170 *     z           =X(3)+t*X(6)
 171 *     t=(X(2)*cos(phi)-X(1)*sin(phi))/(X(4)*sin(phi)-X(5)*cos(phi))
 172       IF (IFL.EQ.2) THEN
 173 *        track not parallel to the phi1 plane ?
 174          UN=X(4)*S1-X(5)*C1
 175          IF (UN.NE.0.) THEN
 176             S=(X(2)*C1-X(1)*S1)/UN
 177             IF (S.GE.0.) THEN
 178                ZI=X(3)+S*X(6)
 179                IF (ABS(ZI).LE.P(3)) THEN
 180                   XI=X(1)+S*X(4)
 181                   YI=X(2)+S*X(5)
 182                   R2=XI**2+YI**2
 183                   IF (P(1)**2.LE.R2.AND.R2.LE.P(2)**2) THEN
 184                      IF ((YI*CFIO-XI*SFIO).LE.0.) THEN
 185                         IF (S.LT.SNXT) SNXT=S
 186                      END IF
 187                   END IF
 188                END IF
 189             END IF
 190          END IF
 191 *        track not parallel to the phi2 plane ?
 192          UN=X(4)*S2-X(5)*C2
 193          IF (UN.NE.0.) THEN
 194             S=(X(2)*C2-X(1)*S2)/UN
 195             IF (S.GE.0.) THEN
 196                ZI=X(3)+S*X(6)
 197                IF (ABS(ZI).LE.P(3)) THEN
 198                   XI=X(1)+S*X(4)
 199                   YI=X(2)+S*X(5)
 200                   R2=XI**2+YI**2
 201                   IF (P(1)**2.LE.R2.AND.R2.LE.P(2)**2) THEN
 202                      IF ((YI*CFIO-XI*SFIO).GE.0.) THEN
 203                         IF (S.LT.SNXT) SNXT=S
 204                      END IF
 205                   END IF
 206                END IF
 207             END IF
 208          END IF
 209       END IF
 210       GO TO 999
 211 *
 212   101 SNXT=S
 213   999 END
 214