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