CommitLineData
fe4da5cc 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 GNPGON (X, PAR, IACT, SNEXT, SNXT, SAFE)
13C.
14C. ******************************************************************
15C. * *
16C. * COMPUTE DISTANCE UP TO INTERSECTION WITH 'PGON' VOLUME, *
17C. * FROM INSIDE POINT X(1-3) ALONG DIRECTION X(4-6) *
18C. * *
19C. * PAR (input) : volume parameters *
20C. * IACT (input) : action flag *
21C. * = 0 Compute SAFE only *
22C. * = 1 Compute SAFE, and SNXT only if SNEXT .GT.new SAFE *
23C. * = 2 Compute both SAFE and SNXT *
24C. * = 3 Compute SNXT only *
25C. * SNEXT (input) : see IACT = 1 *
26C. * SNXT (output) : distance to volume boundary *
27C. * SAFE (output) : shortest distance to any boundary *
28C. * *
29C. * ==>Called by : GNEXT, GTNEXT *
30C. * Author A.McPherson, P.Weidhaas ********* *
31C. * *
32C. ******************************************************************
33C.
34#include "geant321/gconsp.inc"
35
36 DIMENSION X(6), PAR(50), TPAR(50)
37
38C-----------------------------------------------------------------
39
40 SNXT = BIG
41 SAFSEG = BIG
42 R2 = X(1)*X(1) + X(2)*X(2)
43 R = SQRT (R2)
44
45 PDIV = PAR(2) / PAR(3)
46 DELPHI = PDIV * DEGRAD
47 DPHI2 = 0.5 * DELPHI
48 CSDPH2 = COS (DPHI2)
49
50 TPAR(1) = PAR(1)
51 TPAR(2) = PAR(2)
52 TPAR(3) = PAR(3)
53 TPAR(4) = PAR(4)
54 NZ = PAR(4)
55
56 DO 5 I=1, NZ
57 I3 = 3*(I-1)
58 TPAR(5+I3) = PAR(5+I3)
59 TPAR(6+I3) = PAR(6+I3) / CSDPH2
60 TPAR(7+I3) = PAR(7+I3)
61 5 CONTINUE
62
63C**************************************************************************
64C
65C...... Here we start the logic from "GNPCON" (which for reasons of
66C...... efficiency and clarity has been implemented inline).
67C
68C**************************************************************************
69
70 ZMIN = TPAR(5)
71 ZMAX = TPAR(3*NZ+2)
72 SAFZ1 = X(3) - ZMIN
73 SAFZ2 = ZMAX - X(3)
74 SAFEZ = MIN (SAFZ1, SAFZ2)
75
76C
77C...... First determine in which segment the particle is located.
78C
79 DO 10 JPH=8, 3*NZ-1, 3
80 IF (X(3) .LT. TPAR(JPH)) THEN
81 IPH=JPH
82 GO TO 20
83 ENDIF
84 10 CONTINUE
85 IPH = 3*NZ+2
86
87 20 CONTINUE
88C
89C...... The particle is in the segment bounded by z-planes at
90C...... Z1=PAR(IPL) and Z2=PAR(IPH), i.e., Z1 < X(3) < Z2.
91C
92C...... Set parameters for this segment and translate z-coordinate
93C...... of point relative to center of this segment. this is done in
94C...... preparation of invoking the algorithms used in "GNTUBE" and
95C...... "GNCONE" (which for reasons of efficiency and clarity are
96C...... implemented inline).
97C
98 IPL = IPH - 3
99 DZ = 0.5 * (TPAR(IPH) - TPAR(IPL))
100 PT2 = TPAR(IPL+1)
101 PT3 = TPAR(IPL+2)
102 PT4 = TPAR(IPH+1)
103 PT5 = TPAR(IPH+2)
104 PT6 = TPAR(1)
105 PT7 = TPAR(1) + TPAR(2)
106 IF (PT7 .GT. 360.0) PT7 = PT7 - 360.0
107
108 XT3 = X(3) - 0.5 * (TPAR(IPL) + TPAR(IPH))
109
110 SAFZ2 = DZ + XT3
111 ZLENI = 0.5 / DZ
112
113 PHI1 = PT6 * DEGRAD
114 PHI2 = PT7 * DEGRAD
115 IF (PHI2.LE.PHI1) PHI2=PHI2+TWOPI
116
117 IF (IACT .LT. 3) THEN
118
119C -------------------------------------------------
120C | Compute safety-distance 'SAFE' (P.Weidhaas) |
121C -------------------------------------------------
122
123 IFLAG = 2
124 IF (TPAR(2) .EQ. 360.0) IFLAG = 1
125 SAFZ1 = DZ - ABS(XT3)
126 SAFEZ = MIN (SAFEZ,SAFZ1)
127C
128C...... Next determine whether the segment is a tube or a cone.
129C
130 IF (PT2 .NE. PT4) GO TO 50
131 IF (PT3 .NE. PT5) GO TO 50
132
133C*********************************************************
134C
135C...... The segment is a tube; invoke the algorithm
136C...... from routine "GNTUBE" inline to get "SAFER".
137C
138C*********************************************************
139
140 SAFR1 = BIG
141 IF(PT2.GT.0.) SAFR1 = R - PT2
142 SAFR2 = PT3 - R
143 SAFER = MIN (SAFR1, SAFR2)
144
145 IF (IFLAG .EQ. 2) GO TO 70
146
147 GO TO 100
148
149
150 50 CONTINUE
151
152C*********************************************************
153C
154C...... The segment is a cone; invoke the algorithm
155C...... from routine "GNCONE" inline to get "SAFER".
156C
157C*********************************************************
158
159
160C...... Compute radial distance to inner wall.
161
162 IF(PT2+PT4.GT.0.) THEN
163 FACT = (PT4 - PT2) * ZLENI
164 RAD1 = PT2 + FACT * SAFZ2
165 SAFR1 = (R - RAD1) / SQRT(1.0 + FACT*FACT)
166 ELSE
167 SAFR1 = BIG
168 ENDIF
169
170C...... Compute radial distance to outer wall.
171
172 FACT = (PT5 - PT3) * ZLENI
173 RAD2 = PT3 + FACT * SAFZ2
174 SAFR2 = (RAD2 - R) / SQRT(1.0 + FACT*FACT)
175
176 SAFER = MIN (SAFR1, SAFR2)
177
178 IF (IFLAG .EQ. 1) GO TO 100
179
180 70 CONTINUE
181
182C********************************************************************
183C...... Here we handle the case of a phi-segment of a tube or cone.
184C...... in addition to the radial distances (SAFR1, SAFR2) and the
185C...... axial distances (SAFZ1, SAFZ2) we compute here the distance
186C...... to the phi-segment boundary that is closest to the point:
187C
188C...... For each phi-boundary we find the distance from the given
189C...... point to the outer (at R2) point of the segment boundary
190C...... (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
191C...... "SAFSEG" to be the distance to segment PHI1, else we set
192C...... "SAFSEG" to be the distance to segment PHI2.
193C*********************************************************************
194
195
196 COSPH1 = COS (PHI1)
197 SINPH1 = SIN (PHI1)
198 COSPH2 = COS (PHI2)
199 SINPH2 = SIN (PHI2)
200
201C...... Get coordinates of outer endpoints (at R2) of both phi-segments.
202
203 XS1 = R * COSPH1
204 YS1 = R * SINPH1
205 XS2 = R * COSPH2
206 YS2 = R * SINPH2
207
208C...... Get distances (squared) from given point to each endpoint.
209
210 DISTS1 = (X(1) - XS1)**2 + (X(2) - YS1)**2
211 DISTS2 = (X(1) - XS2)**2 + (X(2) - YS2)**2
212
213C...... Get distance to that phi-segment whose endpoint
214C...... is closest to the given point.
215
216 IF (DISTS1 .LE. DISTS2) THEN
217 SAFSEG = ABS(SINPH1 * X(1) - COSPH1 * X(2))
218 ELSE
219 SAFSEG = ABS(SINPH2 * X(1) - COSPH2 * X(2))
220 ENDIF
221
222 100 CONTINUE
223
224
225 IF (SAFER .LE. 0.0) THEN
226
227C---------------------------------------------------------------------------
228C
229C...... Here we handle the case in which SAFER < 0, i.e., the point is
230C...... inside the polygon but outside the inscribed polycone. We must
231C...... do an accurate calculation of "SAFER".
232C
233C---------------------------------------------------------------------------
234
235 FACT = SAFZ2 * ZLENI
236 RAD1 = PT2 + FACT * (PT4 - PT2)
237 RAD2 = PT3 + FACT * (PT5 - PT3)
240
241 IF(X(1).EQ.0.)THEN
242 PHI = ATAN2(X(2), 1.E-8)
243 ELSE
244 PHI = ATAN2(X(2), X(1))
245 ENDIF
246 IF (PHI .LT. PHI1) PHI = PHI + TWOPI
247
248 DIST=0.
249 IF (PHI .GE. PHI1 .AND. PHI .LE. PHI2) THEN
250 PHIREL = PHI - PHI1
251 NSECTR = INT(PHIREL / DELPHI) + 1
252 PHICTR = PHI1 + (2.0*NSECTR - 1.0) * DPHI2
253 COSPHC = COS (PHICTR)
254 SINPHC = SIN (PHICTR)
255
256
257 IF (R2 .GE. RR2) THEN
258 SR2 = RAD2
259 DIST = ABS (COSPHC * X(1) + SINPHC * X(2) - SR2)
260 FACTC = (PT5 - PT3) * ZLENI
261 ELSEIF (R2 .LE. RR1) THEN
262 RIN = RAD1 * CSDPH2
263 SRIN = RIN
264 DIST = ABS (COSPHC * X(1) + SINPHC * X(2) - SRIN)
265 FACTC = (PT4 - PT2) * ZLENI
266 ENDIF
267
268 ENDIF
269
270 SAFER = DIST / SQRT(1.0 + FACTC*FACTC)
271
272 ENDIF
273
274 SAFE = MIN (SAFEZ, SAFER, SAFSEG)
275
276
277 IF (IACT .EQ. 0) GO TO 999
278 IF (IACT .EQ. 1) THEN
279 IF (SNEXT .LT. SAFE) GO TO 999
280 ENDIF
281 ENDIF
282
283
284C ------------------------------------------------
285C | Compute vector-distance 'SNXT' (Nierhaus) |
286C ------------------------------------------------
287
288
289 CALL GNPGO1(X,PAR,SNXT)
290C
291 999 CONTINUE
292 END
293