GEANT321/ggeom/gnocon.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 GNOCON(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) 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 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       PARAMETER (F=0.01745329251994330D0)
  50 #endif
  51 #if defined(CERNLIB_SINGLE)
  52       PARAMETER (F=0.01745329251994330)
  53 #endif
  54       REAL X(6),P(7),SNEXT,SNXT,SAFE
  55       PARAMETER (ONE=1,HALF=ONE/2,ZERO=0)
  56 *
  57 *     this part has to be moved outside the routine
  58       RO1=HALF*(P(4)+P(2))
  59       TG1=HALF*(P(4)-P(2))/P(1)
  60       CR1=ONE/SQRT(ONE+TG1*TG1)
  61       ZV1=1.E10
  62       IF (P(2).NE.P(4)) ZV1=-RO1/TG1
  63       RO2=HALF*(P(5)+P(3))
  64       TG2=HALF*(P(5)-P(3))/P(1)
  65       CR2=ONE/SQRT(ONE+TG2*TG2)
  66       ZV2=1.E10
  67       IF (P(3).NE.P(5)) ZV2=-RO2/TG2
  68       IF (IFL.EQ.2) THEN
  69          P6=P(6)*F
  70          P7=P(7)*F
  71          IF (P7.LT.P6) P7=P7+F*360
  72          C1=COS(P6)
  73          S1=SIN(P6)
  74          C2=COS(P7)
  75          S2=SIN(P7)
  76          FIO=HALF*(P7+P6)
  77          CFIO=COS(FIO)
  78          SFIO=SIN(FIO)
  79          DFI=HALF*(P7-P6)
  80          CDFI=COS(DFI)
  81 *        SDFI=SIN(DFI)
  82       END IF
  83 *
  84       SNXT=1.E10
  85       R   =SQRT(X(1)**2+X(2)**2)
  86       RIN =ABS(TG1*X(3)+RO1)
  87       ROUT=ABS(TG2*X(3)+RO2)
  88 *
  89 *     Compute SAFE radius
  90       IF (IACT.LT.3) THEN
  91          SAF1=(RIN -R)*CR1
  92          SAF2=(R-ROUT)*CR2
  93          SAF3=ABS(X(3))-P(1)
  94          SAF4=0.
  95          IF (IFL.EQ.2.AND.R.GT.0.) THEN
  96             CPSI=(X(1)*CFIO+X(2)*SFIO)/R
  97             IF (CPSI.LT.CDFI) THEN
  98                IF ((X(2)*CFIO-X(1)*SFIO).LE.0.) THEN
  99                   SAF4=ABS(X(1)*S1-X(2)*C1)
 100                ELSE
 101                   SAF4=ABS(X(1)*S2-X(2)*C2)
 102                END IF
 103             END IF
 104          END IF
 105          SAFE=MAX(SAF1,SAF2,SAF3,SAF4)
 106          IF (IACT.EQ.0) GO TO 999
 107          IF (IACT.EQ.1.AND.SNEXT.LE.SAFE) GO TO 999
 108       END IF
 109 *
 110 *     Intersection with z-plane
 111 *     only points outside the z range need to be considered
 112       IF (ABS(X(3)).GE.P(1)) THEN
 113          IF (X(3)*X(6).LT.0.) THEN
 114             S=(ABS(X(3))-P(1))/ABS(X(6))
 115             XI=X(1)+S*X(4)
 116             YI=X(2)+S*X(5)
 117             RIQ=XI**2+YI**2
 118             IF (X(3).LE.0.) THEN
 119                R1Q=P(2)**2
 120                R2Q=P(3)**2
 121             ELSE
 122                R1Q=P(4)**2
 123                R2Q=P(5)**2
 124             END IF
 125             IF (R1Q.LE.RIQ.AND.RIQ.LE.R2Q) THEN
 126                IF (IFL.EQ.1.OR.RIQ.LE.0.)  GO TO 101
 127                CPSI=(XI*CFIO+YI*SFIO)/SQRT(RIQ)
 128                IF (CPSI.GE.CDFI) GO TO 101
 129             END IF
 130          END IF
 131       END IF
 132 *
 133 *     Intersection with cones
 134 *     Intersection point (x,y,z)
 135 *     (x,y,z) is on track: x=X(1)+t*X(4)
 136 *                          y=X(2)+t*X(5)
 137 *                          z=X(3)+t*X(6)
 138 *     (x,y,z) is on cone : x**2 + y**2 = (a*z+b)**2
 139 *
 140 *     (X(4)**2+X(5)**2-(a*x(6))**2)*t**2
 141 *     +2.*(X(1)*X(4)+X(2)*X(5)-a*x(6)*(a*x(3)+b))*t
 142 *     +X(1)**2+X(2)**2-(a*x(3)+b)**2=0
 143 *
 144       T1=X(4)**2+X(5)**2
 145       T2=X(1)*X(4)+X(2)*X(5)
 146       T3=X(1)**2+X(2)**2
 147 *
 148 *     Intersection with the outer cone
 149 *     only points outside the outer cone need to be considered
 150       SR2=1.E10
 151       IF ((ZV2*X(3).GT.ZV2*ZV2).OR.(R.GT.ROUT)) THEN
 152          U=T1-(TG2*X(6))**2
 153          V=T2- TG2*X(6)*(TG2*X(3)+RO2)
 154          W=T3- ROUT*ROUT
 155 *        track not parallel to the cone ?
 156          IF (U.NE.0.) THEN
 157             B=V/U
 158             C=W/U
 159             D=B**2-C
 160             IF (D.GE.0.) THEN
 161 *              compute the smallest root first
 162                IR=-1
 163    11          S=-B+IR*SQRT(D)
 164                IF (S.GE.0.) THEN
 165                   ZI=X(3)+S*X(6)
 166                   IF (ABS(ZI).LE.P(1)) THEN
 167                      IF (IFL.EQ.1.OR.ZI.EQ.ZV2) THEN
 168                         SR2=S
 169                      ELSE
 170                         XI=X(1)+S*X(4)
 171                         YI=X(2)+S*X(5)
 172                         RI=TG2*ZI+RO2
 173                         CPSI=(XI*CFIO+YI*SFIO)/RI
 174                         IF (CPSI.GE.CDFI) SR2=S
 175                      END IF
 176                   END IF
 177                END IF
 178                IF (IR.EQ.-1) THEN
 179                   IF (SR2.EQ.S) GO TO 101
 180 *                 smallest root not ok. Try the biggest one
 181                   IR=+1
 182                   GO TO 11
 183                END IF
 184             END IF
 185          ELSEIF (V.NE.0.) THEN
 186             S=-HALF*W/V
 187             IF (S.GE.0.) THEN
 188                ZI=X(3)+S*X(6)
 189                IF (ABS(ZI).LE.P(1)) THEN
 190                   IF (IFL.EQ.1.OR.ZI.EQ.ZV2)  GO TO 101
 191                   XI=X(1)+S*X(4)
 192                   YI=X(2)+S*X(5)
 193                   RI=TG2*ZI+RO2
 194                   CPSI=(XI*CFIO+YI*SFIO)/RI
 195                   IF (CPSI.GE.CDFI) GO TO 101
 196                END IF
 197             END IF
 198          END IF
 199       END IF
 200 *
 201 *     Intersection with the inner cone
 202       SR1=1.E10
 203       IF (RO1.GT.0.) THEN
 204          U=T1-(TG1*X(6))**2
 205          V=T2- TG1*X(6)*(TG1*X(3)+RO1)
 206          W=T3- RIN*RIN
 207 *        track not parallel to the cone ?
 208          IF (U.NE.0.) THEN
 209             B=V/U
 210             C=W/U
 211             D=B**2-C
 212             IF (D.GE.0.) THEN
 213 *              compute the smallest root first
 214                IR=-1
 215    21          S=-B+IR*SQRT(D)
 216                IF (S.GE.0.) THEN
 217                   ZI=X(3)+S*X(6)
 218                   IF (ABS(ZI).LE.P(1)) THEN
 219                      IF (IFL.EQ.1.OR.ZI.EQ.ZV1) THEN
 220                         SR1=S
 221                      ELSE
 222                         XI=X(1)+S*X(4)
 223                         YI=X(2)+S*X(5)
 224                         RI=TG1*ZI+RO1
 225                         CPSI=(XI*CFIO+YI*SFIO)/RI
 226                         IF (CPSI.GE.CDFI) SR1=S
 227                      END IF
 228                   END IF
 229                END IF
 230                IF (IR.EQ.-1.AND.SR1.GT.9.E9) THEN
 231 *                 smallest root not ok. Try the biggest one
 232                   IR=+1
 233                   GO TO 21
 234                END IF
 235             END IF
 236          ELSEIF (V.NE.0.) THEN
 237             S=-HALF*W/V
 238             IF (S.GE.0.) THEN
 239                ZI=X(3)+S*X(6)
 240                IF (ABS(ZI).LE.P(1)) THEN
 241                   IF (IFL.EQ.1.OR.ZI.EQ.ZV1) THEN
 242                      SR1=S
 243                   ELSE
 244                      XI=X(1)+S*X(4)
 245                      YI=X(2)+S*X(5)
 246                      RI=TG1*ZI+RO1
 247                      CPSI=(XI*CFIO+YI*SFIO)/RI
 248                      IF (CPSI.GE.CDFI) SR1=S
 249                   END IF
 250                END IF
 251             END IF
 252          END IF
 253       END IF
 254 *
 255       SNXT=MIN(SR1,SR2)
 256 *
 257 *     Intersection with phi-planes
 258 *     x=r*cos(phi)=X(1)+t*X(4)
 259 *     y=r*sin(phi)=X(2)+t*X(5)
 260 *     z           =X(3)+t*X(6)
 261 *     t=(X(2)*cos(phi)-X(1)*sin(phi))/(X(4)*sin(phi)-X(5)*cos(phi))
 262 *
 263       IF (IFL.EQ.2) THEN
 264 *        track not parallel to the phi1 plane ?
 265          UN=X(4)*S1-X(5)*C1
 266          IF (UN.NE.0.) THEN
 267             S=(X(2)*C1-X(1)*S1)/UN
 268             IF (S.GE.0.) THEN
 269                ZI=X(3)+S*X(6)
 270                IF (ABS(ZI).LE.P(1)) THEN
 271                   XI=X(1)+S*X(4)
 272                   YI=X(2)+S*X(5)
 273                   RIQ=XI**2+YI**2
 274                   R1Q=(TG1*ZI+RO1)**2
 275                   R2Q=(TG2*ZI+RO2)**2
 276                   IF (R1Q.LE.RIQ.AND.RIQ.LE.R2Q) THEN
 277                      IF ((YI*CFIO-XI*SFIO).LE.0.) THEN
 278                         IF (S.LT.SNXT) SNXT=S
 279                      END IF
 280                   END IF
 281
 282               END IF
 283             END IF
 284          END IF
 285 *        track not parallel to the phi2 plane ?
 286          UN=X(4)*S2-X(5)*C2
 287          IF(UN.NE.0.) THEN
 288             S=(X(2)*C2-X(1)*S2)/UN
 289             IF (S.GE.0.) THEN
 290                ZI=X(3)+S*X(6)
 291                IF (ABS(ZI).LE.P(1)) THEN
 292                   XI=X(1)+S*X(4)
 293                   YI=X(2)+S*X(5)
 294                   RIQ=XI**2+YI**2
 295                   R1Q=(TG1*ZI+RO1)**2
 296                   R2Q=(TG2*ZI+RO2)**2
 297                   IF (R1Q.LE.RIQ.AND.RIQ.LE.R2Q) THEN
 298                      IF ((YI*CFIO-XI*SFIO).GE.0.) THEN
 299                         IF (S.LT.SNXT) SNXT=S
 300                      END IF
 301                   END IF
 302                END IF
 303             END IF
 304          END IF
 305       END IF
 306       GO TO 999
 307 *
 308   101 SNXT=S
 309   999 END