GEANT321/ggeom/gntrp.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:54  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 GNTRP (X, PAR, IACT, SNEXT, SNXT, SAFE)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'TRAP' 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.McPherson  *********                         *
  31 C.    *                                                                *
  32 C.    ******************************************************************
  33 C.
  34 #include "geant321/gconsp.inc"
  35       DIMENSION X(6), PAR(11)
  36 #if defined(CERNLIB_SINGLE)
  37       DIMENSION XT(6), P(11)
  38 #endif
  39 #if !defined(CERNLIB_SINGLE)
  40       DOUBLE PRECISION XT(6), P(11)
  41       DOUBLE PRECISION A,B,C,DSN,TOP1,TOP2,BOT1,BOT2
  42       DOUBLE PRECISION DX1,DX2,X1,X2,X3,DCX1,DCX2,DX,DY
  43       DOUBLE PRECISION SP1,SP2,DCX,DCY,TANG1,TANG2,A1,B1,A2,B2
  44       DOUBLE PRECISION XL1,XL2,XL3,XL4,CY,SY,SX,CX,C1,C2
  45       DOUBLE PRECISION TX,H0,DHDZ,CH
  46       DOUBLE PRECISION CHTAN,SHTAN
  47       DOUBLE PRECISION SN1,SN2,SN3,SN4,SN5,SN6,SN7,SN8,P15
  48 #endif
  49 C.
  50 C.               --------------------------------------
  51 C.
  52       SNXT = BIG
  53       SAFE = 0.0
  54       DO 1 I=1,6
  55    1  XT(I)=X(I)
  56       DO 2 I=1,11
  57    2  P(I)=PAR(I)
  58
  59       P15=0.5/P(1)
  60       DHDZ=(P(4)-P(8))*P15
  61
  62       IF (IACT .LT. 3) THEN
  63
  64
  65 C       -------------------------------------------------
  66 C       |  Compute safety-distance 'SAFE'  (McPherson)  |
  67 C       -------------------------------------------------
  68
  69
  70 C            CALCULATE RECTANGLE ON FACE AT Z=-DZ.
  71
  72         DX1=P(4)*P(7)
  73         X1=-DX1-P(5)
  74         X2=DX1-P(6)
  75         IF(X2.GT.X1) X1=X2
  76         X2=-DX1+P(5)
  77         X3=DX1+P(6)
  78         IF(X3.LT.X2) X2=X3
  79         DCX1=(X1+X2)*0.5
  80         DX1=(X2-X1)*0.5
  81 C
  82 C            CALCULATE RECTANGLE ON FACE AT Z=DZ.
  83 C
  84         DX2=P(8)*P(11)
  85         X1=-DX2-P(9)
  86         X2=DX2-P(10)
  87         IF(X2.GT.X1) X1=X2
  88         X2=-DX2+P(9)
  89         X3=DX2+P(10)
  90         IF(X3.LT.X2) X2=X3
  91         DCX2=(X1+X2)*0.5
  92         DX2=(X2-X1)*0.5
  93
  94 C
  95         TX=P(2)+(DCX2-DCX1)*P15
  96 C
  97 C            CALCULATE LOCAL RECTANGLE.
  98 C
  99         SP1=(P(1)-XT(3))*P15
 100         SP2=(P(1)+XT(3))*P15
 101 C
 102         DY=P(4)*SP1+P(8)*SP2
 103         DX=DX1*SP1+DX2*SP2
 104         DCX=(DCX1+DCX2)*0.5+TX*XT(3)
 105         DCY=P(3)*XT(3)
 106 C
 107 C            CHECK POINT IS INSIDE RECTANGLE.
 108 C
 109         IF(ABS(XT(1)-DCX).GT.DX) GO TO 10
 110 C
 111 C            CALCULATE ANGLE OF YZ PLANES.
 112 C
 113         TANG1=TX+(DX1-DX2)*P15
 114         TANG2=TX-(DX1-DX2)*P15
 115 C
 116 C            CALCULATE SAFETY FROM YZ PLANES.
 117 C
 118         SAF1=(XT(1)-DCX+DX)/SQRT(1.0+TANG1*TANG1)
 119         SAF2=(DX-XT(1)+DCX)/SQRT(1.0+TANG2*TANG2)
 120 C
 121 C            CALCULATE ANGLE OF XZ PLANES.
 122 C
 123         TANG1=P(3)+DHDZ
 124         TANG2=P(3)-DHDZ
 125 C
 126 C            CALCULATE SAFETY FROM XZ PLANES.
 127 C
 128         SAF3=(XT(2)-DCY+DY)/SQRT(1.0+TANG1*TANG1)
 129         SAF4=(DY-XT(2)+DCY)/SQRT(1.0+TANG2*TANG2)
 130 C
 131 C            CALCULATE SAFETY FROM XY PLANES.
 132 C
 133         SAF5=P(1)-ABS(XT(3))
 134 C
 135 C            OVERALL SAFETY.
 136 C
 137         SAFE = MIN(SAF1,SAF2,SAF3,SAF4,SAF5)
 138
 139
 140         IF (IACT .EQ. 0) GO TO 999
 141         IF (IACT .EQ. 1) THEN
 142           IF (SNEXT .LT. SAFE) GO TO 999
 143         ENDIF
 144       ENDIF
 145
 146
 147    10 CONTINUE
 148
 149 C     ------------------------------------------------
 150 C     |  Compute vector-distance 'SNXT' (McPherson)  |
 151 C     ------------------------------------------------
 152
 153
 154
 155 C           FIRST FIND S TO Z LIMITS.
 156
 157       SN1=BIG
 158       SN2=BIG
 159       SN3=BIG
 160       SN4=BIG
 161 C
 162       IF(XT(6).EQ.0.0) GOTO 15
 163       SN1=(P(1)-XT(3))/XT(6)
 164       SN2=-(P(1)+XT(3))/XT(6)
 165       IF(SN1.GE.0.) GOTO 15
 166       SNX=SN1
 167       SN1=SN2
 168       SN2=SNX
 169 15    CONTINUE
 170 C
 171 C           NOW Y LIMITS.
 172 C
 173       H0=(P(4)+P(8))*0.5
 174 C
 175       TOP1=H0-DHDZ*XT(3)
 176       TOP2=XT(2)-XT(3)*P(3)
 177       BOT1=XT(5)-XT(6)*P(3)
 178       BOT2=DHDZ*XT(6)
 179 C
 180       IF(BOT1+BOT2.NE.0.0) SN3=(TOP1-TOP2)/(BOT1+BOT2)
 181       IF(BOT2-BOT1.NE.0.0) SN4=(TOP1+TOP2)/(BOT2-BOT1)
 182 C
 183 C           NOW X LIMTS.
 184 C
 185       XL1=(P(5)+P(6)+P(9)+P(10))*0.25
 186       XL2=(P(5)-P(6)+P(9)-P(10))*0.25
 187       XL3=(P(5)+P(6)-P(9)-P(10))*0.25
 188       XL4=(P(5)-P(6)-P(9)+P(10))*0.25
 189 C
 190       CY=XT(2)-XT(3)*P(3)
 191       SY=XT(5)-XT(6)*P(3)
 192       CH=H0-DHDZ*XT(3)
 193 C
 194       A1=XL4*XT(6)*SY+XL3*XT(6)*XT(6)*DHDZ
 195       B1=XL4*(XT(6)*CY+XT(3)*SY)-XL3*XT(6)*(CH-XT(3)*DHDZ)-XL2*P(1)*SY-
 196      +XL1*XT(6)*DHDZ*P(1)
 197       C1=XL4*CY*XT(3)-XL3*XT(3)*CH-XL2*CY*P(1)+XL1*CH*P(1)
 198 C
 199       CHTAN=(P(7)*P(4)+P(11)*P(8))*0.5
 200       SHTAN=-(P(7)*P(4)-P(11)*P(8))*P15
 201       CHTAN=CHTAN+SHTAN*XT(3)
 202       SHTAN=SHTAN*XT(6)
 203       CX=XT(1)-XT(3)*P(2)
 204       SX=XT(4)-XT(6)*P(2)
 205 C
 206       A2=-P(1)*(DHDZ*XT(6)*SX+SY*SHTAN)
 207       B2=P(1)*(CH*SX-DHDZ*XT(6)*CX-CY*SHTAN-SY*CHTAN)
 208       C2=P(1)*(CH*CX-CY*CHTAN)
 209 C
 210       ISN56=0
 211       A=A1+A2
 212       B=B1+B2
 213       C=C1+C2
 214       IF(B*B-A*C*4.0.LT.0.0) GO TO 20
 215       IF(ABS(A).LT.1.0E-7) GO TO 19
 216       DSN=SQRT(B*B-A*C*4.0)
 217       SN5=(-B+DSN)*0.5/A
 218       SN6=(-B-DSN)*0.5/A
 219       ISN56=1
 220       GO TO 20
 221    19 CONTINUE
 222       IF(ABS(B).LT.1.0E-5) GO TO 20
 223       ISN56=1
 224       SN5=-C/B
 225       SN6=-C/B
 226    20 CONTINUE
 227 C
 228       ISN78=0
 229       A=A1-A2
 230       B=B1-B2
 231       C=C1-C2
 232       IF(B*B-A*C*4.0.LT.0.0) GO TO 30
 233       DSN=SQRT(B*B-A*C*4.0)
 234       IF(ABS(A).LT.1.0E-7) GO TO 25
 235       SN7=(-B+DSN)*0.5/A
 236       SN8=(-B-DSN)*0.5/A
 237       ISN78=1
 238       GO TO 30
 239    25 CONTINUE
 240       IF(ABS(B).LT.1.0E-5) GO TO 30
 241       ISN78=1
 242       SN7=-C/B
 243       SN8=-C/B
 244    30 CONTINUE
 245 C
 246       IF(SN2.GT.0.0.AND.SN2.LT.SN1) SN1=SN2
 247       IF(SN3.GT.0.0.AND.SN3.LT.SN1) SN1=SN3
 248       IF(SN4.GT.0.0.AND.SN4.LT.SN1) SN1=SN4
 249       IF(ISN56.EQ.0) GO TO 40
 250       IF(SN5.GT.0.0.AND.SN5.LT.SN1) SN1=SN5
 251       IF(SN6.GT.0.0.AND.SN6.LT.SN1) SN1=SN6
 252    40 CONTINUE
 253       IF(ISN78.EQ.0) GO TO 50
 254       IF(SN7.GT.0.0.AND.SN7.LT.SN1) SN1=SN7
 255       IF(SN8.GT.0.0.AND.SN8.LT.SN1) SN1=SN8
 256    50 CONTINUE
 257 C
 258       SNXT=SN1
 259 C
 260   999 CONTINUE
 261       END