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