GEANT321/ggeom/gnpgo1.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 *********************************************************************
  13 ***** GNPGO1 ********************************************************
  14 *
  15 *     GNPGO1 ... 15-AUG-1991
  16 *     Version 1.1
  17 *     Rolf Nierhaus
  18 *
  19 *********************************************************************
  20 *
  21 *     Copyright   CERN,   Geneva  1991  -  Copyright  and  any  other
  22 *     appropriate legal protection of  these  computer  programs  and
  23 *     associated  documentation  reserved  in  all  countries  of the
  24 *     world.
  25 *
  26 *********************************************************************
  27 *
  28 *          Subroutine  GNPGO1 is called by GNPGON for the computation
  29 *     of SNXT, the distance from a point P  along  a  track  T  to  a
  30 *     boundary  surface  of a Geant volume V of shape PGON. The point
  31 *     P is inside the volume V.
  32 *
  33 *          V  is  generally  a composite volume consisting of several
  34 *     sections. The sections have  boundary  surfaces  orthogonal  to
  35 *     the   Z-axis.   Each  section  consists  generally  of  several
  36 *     sectors. Each sector is an  "elementary"  convex  volume.  This
  37 *     package  assumes it is either a hexahedron or a pentahedron. If
  38 *     it is a pentahedron, it has 6 vertices, of  which  two  are  on
  39 *     the  Z-axis.  All  sectors  of  the same section are congruent.
  40 *     Each section has the same number of sectors.
  41 *
  42 *          We  describe each surface by 6 parameters: the first three
  43 *     are   the   coordinates   of   a   point   on    the    surface
  44 *     XS(I),YS(I),ZS(I),  the  other  three are the components of the
  45 *     normal vector of the surface XN(I),YN(I),ZN(I). I is the  index
  46 *     of  the surface. We consider only one sector at a time, and the
  47 *     number of boundary  surfaces  is  never  larger  then  6.  Each
  48 *     surface  divides  the  space  into  two  regions:  the positive
  49 *     region and the negative region. We choose the direction of  the
  50 *     normal  vectors  of the boundary surfaces such that the bounded
  51 *     volume is within the positive region of each surface, that  is,
  52 *     the normal vector is pointing to the inside of the volume.
  53 *
  54 ***** Subroutine GNPGO1 *************************** 15-AUG-1991 *****
  55       SUBROUTINE GNPGO1(X,P,SNXT)
  56 #if !defined(CERNLIB_SINGLE)
  57       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  58 #endif
  59       REAL X(6),P(49),SNXT
  60 #if defined(CERNLIB_SINGLE)
  61       PARAMETER (F=0.01745329251994330,TPI=6.283185307179586)
  62 #endif
  63 #if !defined(CERNLIB_SINGLE)
  64       PARAMETER (F=0.01745329251994330D0,TPI=6.283185307179586D0)
  65 #endif
  66       PARAMETER (ONE=1,HALF=ONE/2,ZERO=0.)
  67       DIMENSION XS(6), YS(6), ZS(6), XN(6), YN(6), ZN(6)
  68       LOGICAL FLAG, FLG
  69       XP=X(1)
  70       YP=X(2)
  71       ZP=X(3)
  72       XD=X(4)
  73       YD=X(5)
  74       ZD=X(6)
  75       IMAX=P(4)-.5
  76 *     IMAX   ->   number of Z-sections
  77       JMAX=P(3)+.5
  78 *     JMAX   ->   number of Phi-sectors
  79       SNXT=0.
  80    10 CONTINUE
  81 *     Find current elementary volume
  82       IF (ZP.LE.P(5)) RETURN
  83 *     Current point (XP,YP,ZP) is below first section.
  84       DO 20 II=1,IMAX
  85          IF (ZP.LT.P(5+3*II)) THEN
  86             I=II
  87             GO TO 30
  88          ENDIF
  89    20 CONTINUE
  90       RETURN
  91 *     Current point (XP,YP,ZP) is above last section.
  92    30 CONTINUE
  93       IF (XP.EQ.0..AND.YP.EQ.0.) XP=1.E-20
  94       PHI=ATAN2(YP,XP)
  95       IF (PHI.LT.0.) PHI=PHI+TPI
  96       P1=F*P(1)
  97       PHI1=PHI-P1
  98       IF (PHI1.LT.0.) PHI1=PHI1+TPI
  99       IF (PHI1.GE.TPI) PHI1=PHI1-TPI
 100       IF (JMAX.EQ.1) THEN
 101           IF (ABS(PHI1-TPI).LT.1D-7) PHI1=0.
 102       ENDIF
 103       J=PHI1*P(3)/(F*P(2))+ONE
 104       IF (P(2).EQ.360.) THEN
 105          IF (J.LT.1) THEN
 106             J=J+JMAX
 107          ELSEIF (JMAX.LT.J) THEN
 108             J=J-JMAX
 109          END IF
 110       END IF
 111       IF (JMAX.LT.J.OR.J.LT.1) RETURN
 112 *     Current point is outside Phi-range.
 113 C*****  Code Expanded From Routine:  GNPGO2
 114 *     GNPGO2 finds the vector distance to the boundary surface
 115 *     of the current elementary volume.
 116 *     I is Z-section, J is Phi-sector.
 117 *     GNPGO2 calls GNPGO4 five or six times for the storage of
 118 *     the surface coefficients of its boundary surfaces.
 119 *
 120       INDEX = 2 + 3*I
 121       Z1 = P(INDEX)
 122       D1N = P(INDEX+1)
 123       D1X = P(INDEX+2)
 124       Z2 = P(INDEX+3)
 125       D2N = P(INDEX+4)
 126       D2X = P(INDEX+5)
 127       ZM = HALF*(Z1 + Z2)
 128       P11X = F*(P(1)+(J-1)*P(2)/JMAX)
 129       P2 = F*(P(1)+J*P(2)/JMAX)
 130       PP = HALF*(P11X + P2)
 131       COSP = COS(PP)
 132       SINP = SIN(PP)
 133       DMX = HALF*(D1X + D2X)
 134       DMN = HALF*(D1N + D2N)
 135       THX = ATAN((D2X - D1X)/(Z2 - Z1))
 136       COSTHX = COS(THX)
 137       SINTHX = SIN(THX)
 138       XNN = -COSP*COSTHX
 139       YNN = -SINP*COSTHX
 140 C*****  Code Expanded From Routine:  GNPGO4
 141 *     Store surface coefficients
 142       XS(5) = DMX*COSP
 143       YS(5) = DMX*SINP
 144       ZS(5) = ZM
 145       XN(5) = XNN
 146       YN(5) = YNN
 147       ZN(5) = SINTHX
 148 C*****  End of Code Expanded From Routine:  GNPGO4
 149 C*****  Code Expanded From Routine:  GNPGO9
 150 *          Logical   function   GNPGO9  returns  TRUE  if  the  point
 151 *     (XP,YP,ZP) is within the positive region of  the  surface  with
 152 *     index   I.   This   is  the  case  if  the  scalar  product  of
 153 *     (XP-XS,YP-YS,ZP-ZS) and (XN,YN,ZN) is positive (or zero).
 154       RESULT=(XP-XS(5))*XN(5)+(YP-YS(5))*YN(5)+(ZP-ZS(5))*ZN(5)
 155       FLG = 0. .LE. RESULT
 156       IF (.NOT.FLG) GO TO 50
 157       ISMAX = 5
 158       IF (DMN .NE. 0.) THEN
 159          ISMAX = 6
 160          THN = ATAN((D2N - D1N)/(Z2 - Z1))
 161          COSTHN = COS(THN)
 162          SINTHN = SIN(THN)
 163          XNN = COSP*COSTHN
 164          YNN = SINP*COSTHN
 165 C*****  Code Expanded From Routine:  GNPGO4
 166          XS(6) = DMN*COSP
 167          YS(6) = DMN*SINP
 168          ZS(6) = ZM
 169          XN(6) = XNN
 170          YN(6) = YNN
 171          ZN(6) = -SINTHN
 172 C*****  End of Code Expanded From Routine:  GNPGO4
 173 C*****  Code Expanded From Routine:  GNPGO9
 174          RESULT=(XP-XS(6))*XN(6)+(YP-YS(6))*YN(6)+(ZP-ZS(6))*ZN(6)
 175          FLG = 0. .LE. RESULT
 176          IF (.NOT.FLG) GO TO 50
 177       ENDIF
 178 C*****  Code Expanded From Routine:  GNPGO4
 179       XS(1) = ZERO
 180       YS(1) = ZERO
 181       ZS(1) = Z1
 182       XN(1) = ZERO
 183       YN(1) = ZERO
 184       ZN(1) = ONE
 185 C*****  End of Code Expanded From Routine:  GNPGO4
 186 C*****  Code Expanded From Routine:  GNPGO4
 187       XS(2) = ZERO
 188       YS(2) = ZERO
 189       ZS(2) = Z2
 190       XN(2) = ZERO
 191       YN(2) = ZERO
 192       ZN(2) = -ONE
 193 C*****  End of Code Expanded From Routine:  GNPGO4
 194 C*****  Code Expanded From Routine:  GNPGO4
 195       XS(3) = ZERO
 196       YS(3) = ZERO
 197       ZS(3) = ZM
 198       XN(3) = -SIN(P11X)
 199       YN(3) = COS(P11X)
 200       ZN(3) = ZERO
 201 C*****  End of Code Expanded From Routine:  GNPGO4
 202 C*****  Code Expanded From Routine:  GNPGO4
 203       XS(4) = ZERO
 204       YS(4) = ZERO
 205       ZS(4) = ZM
 206       XN(4) = SIN(P2)
 207       YN(4) = -COS(P2)
 208       ZN(4) = ZERO
 209 C*****  End of Code Expanded From Routine:  GNPGO4
 210 C*****  Code Expanded From Routine:  GNPGO5
 211 *     Vector distance to volume boundary
 212       SNXT1 = 1.E10
 213       DO 40  IS = 1, ISMAX
 214 C*****  Code Expanded From Routine:  GNPGO7
 215 *          To  find  the  distance  from  a  point (XP,YP,ZP) along a
 216 *     track  with  direction  cosines   (XD,YD,ZD)   to   a   surface
 217 *     (XS,YS,ZS)(XN,YN,ZN),  we  compute  first the scalar product of
 218 *     the  vector  (XS-XP,YS-YP,ZS-ZP)   with   the   normal   vector
 219 *     (XN,YN,ZN),  then  the scalar product of the vectors (XD,YD,ZD)
 220 *     and (XN,YN,ZN).  The  first  scalar  product  is  the  shortest
 221 *     distance  from  the  point  to  the  plane,  the  second scalar
 222 *     product is the cosine of the angle between the  track  and  the
 223 *     plane  normal.  The  quotient  is  the vector distance. If this
 224 *     vector distance is  positive  (or  zero)  we  set  the  logical
 225 *     variable  FLAG  TRUE.  GNPGO7  is  called with three parameters
 226 *     I,FLAG and DIST. I is the index of the  surface,  and  DIST  is
 227 *     the vector distance if FLAG is TRUE.
 228          SPPMSN = (XP - XS(IS))*XN(IS) + (YP - YS(IS))*YN(IS) + (ZP -
 229      +   ZS (IS))*ZN(IS)
 230          SPDN = XD*XN(IS) + YD*YN(IS) + ZD*ZN(IS)
 231          IF (SPDN .EQ. 0.) THEN
 232             DIST1 = 0.
 233          ELSE
 234             DIST1 = -(SPPMSN + .0001)/SPDN
 235          ENDIF
 236          FLAG = 0. .LT. DIST1
 237 C*****  End of Code Expanded From Routine:  GNPGO7
 238          IF (FLAG) SNXT1 = MIN(DIST1,SNXT1)
 239    40 CONTINUE
 240    50 CONTINUE
 241 C*****  End of Code Expanded From Routine:  GNPGO2
 242       IF (FLG) THEN
 243          SNXT=SNXT+SNXT1
 244          XP=XP+SNXT1*XD
 245          YP=YP+SNXT1*YD
 246          ZP=ZP+SNXT1*ZD
 247 *     The current point (XP,YP,ZP) is propagated along the track
 248 *     to the boundary of the current elementary volume.
 249          GO TO 10
 250       ENDIF
 251       END