GEANT321/ggeom/gncone.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:51  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 GNCONE(X,P,IACT,IFL,SNEXT,SNXT,SAFE)
  13 C.    ******************************************************************
  14 C.    *                                                                *
  15 C.    *      Compute distance to intersection with boundary surface of *
  16 C     *      volume CONE or CONS, from point X(1),X(2),X(3) inside     *
  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 CONE, 2 for PHI segmented CONE     *
  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.    * 'CONE'    is a conical tube. It has 5 parameters :             *
  37 C.    *           the half length in z,                                *
  38 C.    *           the inside and outside radii at the low z limit,     *
  39 C.    *           and those at the high z limit.                       *
  40 C.    * 'CONS'    is a phi segment of a  conical tube.  It has 7       *
  41 C.    *           parameters, the same 5 as 'CONE' plus the phi limits.*
  42 C.    *           The segment starts at the first limit and  includes  *
  43 C.    *           increasing phi  value up  to the  second limit  or   *
  44 C.    *           that plus 360 degrees.                               *
  45 C.    *                                                                *
  46 C.    ******************************************************************
  47 #if !defined(CERNLIB_SINGLE)
  48       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  49       REAL SNEXT, SNXT, SAFE
  50       PARAMETER (F=0.01745329251994330D0)
  51 #endif
  52 #if defined(CERNLIB_SINGLE)
  53       PARAMETER (F=0.01745329251994330)
  54 #endif
  55       REAL X(6),P(7)
  56 *
  57 *     this part has to be moved outside the routine
  58       RO1=0.5*(P(4)+P(2))
  59       TG1=0.5*(P(4)-P(2))/P(1)
  60       CR1=1./SQRT(1.+TG1*TG1)
  61       RO2=0.5*(P(5)+P(3))
  62       TG2=0.5*(P(5)-P(3))/P(1)
  63       CR2=1./SQRT(1.+TG2*TG2)
  64       IF (IFL.EQ.2) THEN
  65          P6=P(6)*F
  66          P7=P(7)*F
  67          IF (P7.LT.P6) P7=P7+F*360.
  68          C1=COS(P6)
  69          S1=SIN(P6)
  70          C2=COS(P7)
  71          S2=SIN(P7)
  72          FIO=0.5*(P7+P6)
  73          CFIO=COS(FIO)
  74          SFIO=SIN(FIO)
  75       END IF
  76 *
  77       SNXT=1.E10
  78       R   =SQRT(X(1)**2+X(2)**2)
  79       RIN =TG1*X(3)+RO1
  80       ROUT=TG2*X(3)+RO2
  81 *
  82 *     Compute SAFE radius
  83       IF (IACT.LT.3) THEN
  84          SAF1=(R -RIN)*CR1
  85          SAF2=(ROUT-R)*CR2
  86          SAF3=P(1)-ABS(X(3))
  87          SAF4=1.E10
  88          IF (IFL.EQ.2) THEN
  89             IF ((X(2)*CFIO-X(1)*SFIO).LE.0.) THEN
  90                SAF4=ABS(X(1)*S1-X(2)*C1)
  91             ELSE
  92                SAF4=ABS(X(1)*S2-X(2)*C2)
  93             END IF
  94          END IF
  95          SAFE=MIN(SAF1,SAF2,SAF3,SAF4)
  96          IF (IACT.EQ.0) GO TO 999
  97          IF (IACT.EQ.1.AND.SNEXT.LE.SAFE) GO TO 999
  98       END IF
  99 *
 100 *     Intersection with z-plane
 101       IF (X(6).GT. 1.E-20) THEN
 102          SZ= (P(1)-X(3))/X(6)
 103       ELSEIF (X(6).LT.-1.E-20) THEN
 104          SZ=-(P(1)+X(3))/X(6)
 105       ELSE
 106          SZ= 1.E10
 107       END IF
 108 *
 109 *     Intersection with cones
 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 cone : x**2 + y**2 = (a*z+b)**2
 115 *
 116 *     (X(4)**2+X(5)**2-(a*x(6))**2)*t**2
 117 *     +2.*(X(1)*X(4)+X(2)*X(5)-a*x(6)*(a*x(3)+b))*t
 118 *     +X(1)**2+X(2)**2-(a*x(3)+b)**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 *
 124 *     Intersection with the inner cone
 125       SR1=1.E10
 126       IF (RO1.GT.0.) THEN
 127          U=T1-(TG1*X(6))**2
 128          V=T2- TG1*X(6)*RIN
 129          W=T3- RIN*RIN
 130 *        track not parallel to the cone ?
 131          IF (U.NE.0.) THEN
 132             B=V/U
 133             C=W/U
 134             D=B**2-C
 135             IF (D.GE.0.) THEN
 136                DS=SQRT(D)
 137                IF (DS.GE.ABS(B)) THEN
 138                   SR1= DS-B
 139                ELSEIF (B.LE.0.)  THEN
 140                   SR1=-DS-B
 141                END IF
 142             END IF
 143          ELSEIF (V.LT.0.) THEN
 144             SR1=-0.5*W/V
 145          END IF
 146       END IF
 147 *
 148 *     Intersection with the outer cone
 149       SR2=1.E10
 150       U=T1-(TG2*X(6))**2
 151       V=T2- TG2*X(6)*ROUT
 152       W=T3- ROUT*ROUT
 153 *     track not parallel to the cone ?
 154       IF (U.NE.0.) THEN
 155          B=V/U
 156          C=W/U
 157          D=B**2-C
 158          IF (D.GE.0.) THEN
 159             DS=SQRT(D)
 160             IF (DS.GE.ABS(B)) THEN
 161                SR2= DS-B
 162             ELSEIF (B.LE.0.)  THEN
 163                SR2=-DS-B
 164             END IF
 165          END IF
 166       ELSEIF (V.GT.0.) THEN
 167          SR2=-0.5*W/V
 168       END IF
 169 *
 170 *     Intersection with phi-planes
 171 *     x=r*cos(phi)=X(1)+t*X(4)
 172 *     y=r*sin(phi)=X(2)+t*X(5)
 173 *     z           =X(3)+t*X(6)
 174 *     t=(X(2)*cos(phi)-X(1)*sin(phi))/(X(4)*sin(phi)-X(5)*cos(phi))
 175 *
 176       SFI1=1.E10
 177       SFI2=1.E10
 178       IF (IFL.EQ.2) THEN
 179 *        track not parallel to the phi1 plane ?
 180          UN=X(4)*S1-X(5)*C1
 181          IF (UN.NE.0.) THEN
 182             S=(X(2)*C1-X(1)*S1)/UN
 183             IF (S.GE.0.) THEN
 184                XI=X(1)+S*X(4)
 185                YI=X(2)+S*X(5)
 186                IF ((YI*CFIO-XI*SFIO).LE.0.) SFI1=S
 187             END IF
 188          END IF
 189 *        track not parallel to the phi2 plane ?
 190          UN=X(4)*S2-X(5)*C2
 191          IF (UN.NE.0.) THEN
 192             S=(X(2)*C2-X(1)*S2)/UN
 193             IF (S.GE.0.) THEN
 194                XI=X(1)+S*X(4)
 195                YI=X(2)+S*X(5)
 196                IF ((YI*CFIO-XI*SFIO).GE.0.) SFI2=S
 197             END IF
 198          END IF
 199       END IF
 200 *
 201       SNXT=MIN(SR1,SR2,SZ,SFI1,SFI2)
 202   999 END
 203