Moving to double precision and introducing an additional security check
[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gnpgon.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1999/05/18 15:55:17  fca
6 * AliRoot sources
7 *
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)
16 C.
17 C.    ******************************************************************
18 C.    *                                                                *
19 C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'PGON' VOLUME,  *
20 C.    *       FROM INSIDE POINT X(1-3) ALONG DIRECTION X(4-6)          *
21 C.    *                                                                *
22 C.    *       PAR   (input)  : volume parameters                       *
23 C.    *       IACT  (input)  : action flag                             *
24 C.    *         = 0  Compute SAFE only                                 *
25 C.    *         = 1  Compute SAFE, and SNXT only if SNEXT .GT.new SAFE *
26 C.    *         = 2  Compute both SAFE and SNXT                        *
27 C.    *         = 3  Compute SNXT only                                 *
28 C.    *       SNEXT (input)  : see IACT = 1                            *
29 C.    *       SNXT  (output) : distance to volume boundary             *
30 C.    *       SAFE  (output) : shortest distance to any boundary       *
31 C.    *                                                                *
32 C.    *    ==>Called by : GNEXT, GTNEXT                                *
33 C.    *         Author  A.McPherson,  P.Weidhaas  *********            *
34 C.    *                                                                *
35 C.    ******************************************************************
36 C.
37       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
38 #include "geant321/gconsp.inc"
39  
40       REAL X(6), PAR(50), SNEXT, SNXT, SAFE
41       DIMENSION TPAR(50)
42  
43 C-----------------------------------------------------------------
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  
68 C**************************************************************************
69 C
70 C......  Here we start the logic from "GNPCON" (which for reasons of
71 C......  efficiency and clarity has been implemented inline).
72 C
73 C**************************************************************************
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  
81 C
82 C......  First determine in which segment the particle is located.
83 C
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
93 C
94 C......  The particle is in the segment bounded by z-planes at
95 C......  Z1=PAR(IPL) and Z2=PAR(IPH), i.e.,  Z1 < X(3) < Z2.
96 C
97 C......  Set parameters for this segment and translate z-coordinate
98 C......  of point relative to center of this segment. this is done in
99 C......  preparation of invoking the algorithms used in "GNTUBE" and
100 C......  "GNCONE" (which for reasons of efficiency and clarity are
101 C......  implemented inline).
102 C
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  
124 C       -------------------------------------------------
125 C       |  Compute safety-distance 'SAFE' (P.Weidhaas)  |
126 C       -------------------------------------------------
127  
128         IFLAG = 2
129         IF (TPAR(2) .EQ. 360.0)  IFLAG = 1
130         SAFZ1 = DZ - ABS(XT3)
131         SAFEZ = MIN (SAFEZ,SAFZ1)
132 C
133 C......  Next determine whether the segment is a tube or a cone.
134 C
135         IF (PT2 .NE. PT4) GO TO 50
136         IF (PT3 .NE. PT5) GO TO 50
137  
138 C*********************************************************
139 C
140 C......  The segment is a tube; invoke the algorithm
141 C......  from routine "GNTUBE" inline to get "SAFER".
142 C
143 C*********************************************************
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  
157 C*********************************************************
158 C
159 C......  The segment is a cone; invoke the algorithm
160 C......  from routine "GNCONE" inline to get "SAFER".
161 C
162 C*********************************************************
163  
164  
165 C......  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  
175 C......  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  
187 C********************************************************************
188 C......  Here we handle the case of a phi-segment of a tube or cone.
189 C......  in addition to the radial distances (SAFR1, SAFR2) and the
190 C......  axial distances (SAFZ1, SAFZ2) we compute here the distance
191 C......  to the phi-segment boundary that is closest to the point:
192 C
193 C......  For each phi-boundary we find the distance from the given
194 C......  point to the outer (at R2) point of the segment boundary
195 C......  (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
196 C......  "SAFSEG" to be the distance to segment PHI1, else we set
197 C......  "SAFSEG" to be the distance to segment PHI2.
198 C*********************************************************************
199  
200  
201           COSPH1  = COS (PHI1)
202           SINPH1  = SIN (PHI1)
203           COSPH2  = COS (PHI2)
204           SINPH2  = SIN (PHI2)
205  
206 C......  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  
213 C......  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  
218 C......  Get distance to that phi-segment whose endpoint
219 C......  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  
232 C---------------------------------------------------------------------------
233 C
234 C......  Here we handle the case in which  SAFER < 0, i.e., the point is
235 C......  inside the polygon but outside the inscribed polycone. We must
236 C......  do an accurate calculation of "SAFER".
237 C
238 C---------------------------------------------------------------------------
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
271             ELSE
272               PRINT *, 'GNPGON: FACTC not calculated!!!'
273               FACTC=0
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  
292 C     ------------------------------------------------
293 C     |  Compute vector-distance 'SNXT' (Nierhaus)  |
294 C     ------------------------------------------------
295  
296  
297       CALL GNPGO1(X,PAR,SNXT)
298 C
299   999 CONTINUE
300       END
301