]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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 |