This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gnopco.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:20:52  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.29  by  S.Giani
11 *-- Author :
12       SUBROUTINE GNOPCO (X, PAR, IACT, SNEXT, SNXT, SAFE)
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'PCON' VOLUME,  *
17 C.    *       FROM OUTSIDE POINT X(1-3) ALONG DIRECTION X(4-6)         *
18 C.    *                                                                *
19 C.    *       PAR   (input)  : volume parameters                       *
20 C.    *       IACT  (input)  : action flag                             *
21 C.    *         = 0  Compute SAFE only                                 *
22 C.    *         = 1  Compute SAFE, and SNXT only if SNEXT .GT.new SAFE *
23 C.    *         = 2  Compute both SAFE and SNXT                        *
24 C.    *         = 3  Compute SNXT only                                 *
25 C.    *       SNEXT (input)  : see IACT = 1                            *
26 C.    *       SNXT  (output) : distance to volume boundary             *
27 C.    *       SAFE  (output) : shortest distance to any boundary       *
28 C.    *                                                                *
29 C.    *    ==>Called by : GNEXT, GTNEXT                                *
30 C.    *         Author  A.McPherson,  P.Weidhaas  *********            *
31 C.    *                                                                *
32 C.    ******************************************************************
33 C.
34 #include "geant321/gconsp.inc"
35       DIMENSION X(6), PAR(9), PT(7), XT(6)
36  
37       EQUIVALENCE (PT(1),PT1), (PT(2),PT2), (PT(3),PT3), (PT(4),PT4)
38       EQUIVALENCE (PT(5),PT5), (PT(6),PT6), (PT(7),PT7)
39 C.
40 C.   --------------------------------------------------------
41 C.
42  
43       SNXT = BIG
44       SAFE = BIG
45       R2 = X(1)*X(1) + X(2)*X(2)
46       R  = SQRT (R2)
47       NZ  = PAR(3)
48       NSEC= NZ - 1
49       PT6=PAR(1)
50       PT7=PAR(1)+PAR(2)
51  
52       IFLAG = 2
53       IF(PAR(2).EQ.360.0) IFLAG = 1
54  
55       IF (IACT .LT. 3) THEN
56  
57 C       -------------------------------------------------
58 C       |  Compute safety-distance 'SAFE' (P.Weidhaas)  |
59 C       -------------------------------------------------
60  
61  
62         SAFEZ  = 0.0
63 C
64 C......  Obtain axial distance "SAFEZ".
65 C
66         ZMIN   = PAR(4)
67         ZMAX   = PAR(3*NZ+1)
68  
69         IF (X(3) .LT. ZMIN) THEN
70           SAFEZ  = ZMIN - X(3)
71         ELSEIF (X(3) .GT. ZMAX) THEN
72           SAFEZ  = X(3) - ZMAX
73         ENDIF
74  
75  
76 C----------------------------------------------------
77 C......  Prepare parameters for PHI-segmented cone:
78 C----------------------------------------------------
79  
80       IF (IFLAG .EQ. 2) THEN
81  
82         PHI1 = MOD (PT6+360.0 , 360.0)
83         PHI2 = MOD (PT7+360.0 , 360.0)
84         IF ( X(1).NE.0.0 .OR. X(2).NE.0.0 ) THEN
85            PHI = ATAN2( X(2), X(1) ) * RADDEG
86         ELSE
87            PHI = 0.0
88         ENDIF
89         PHI  = MOD (PHI+360.0 , 360.0)
90  
91         SINPH1 = SIN(PHI1*DEGRAD)
92         COSPH1 = COS(PHI1*DEGRAD)
93         SINPH2 = SIN(PHI2*DEGRAD)
94         COSPH2 = COS(PHI2*DEGRAD)
95  
96  
97 C......  Set flag "IOPENP" if point (X(1),X(2),X(3)) lies in the
98 C......  PHI-opening.
99  
100         IOPENP = 0
101         IF (PHI2 .GT. PHI1) THEN
102           IF (PHI.GT.PHI2 .OR. PHI.LT.PHI1) IOPENP = 1
103         ELSE
104           IF (PHI.GT.PHI2 .AND. PHI.LT.PHI1) IOPENP = 1
105         ENDIF
106       ENDIF
107  
108  
109  
110 C------------------   Start of loop over Z-sections  --------------
111  
112         IPZ2=4
113  
114         DO 150 IS=1,NSEC
115  
116         IPZ1=IPZ2
117         IPZ2=IPZ1+3
118  
119         SAFSEG = 0.0
120         SAFER  = 0.0
121  
122         XT3=X(3) - 0.5 * (PAR(IPZ1)+PAR(IPZ2))
123         DZ   = 0.5 * (PAR(IPZ2)-PAR(IPZ1))
124  
125         PT2  = PAR(IPZ1+1)
126         PT3  = PAR(IPZ1+2)
127         PT4  = PAR(IPZ2+1)
128         PT5  = PAR(IPZ2+2)
129 C**** check DZ=0 segments
130         IF (DZ.LE.0.) THEN
131            IF ((R-PT2)*(R-PT4).LE.0. .OR. (R-PT3)*(R-PT5).LE.0.) THEN
132               SAFER = ABS(XT3)
133            ELSE
134               GO TO 150
135            ENDIF
136            GO TO 100
137         ENDIF
138  
139         IF (PT2 .NE. PT4) GO TO 50
140         IF (PT3 .NE. PT5) GO TO 50
141  
142 C**********************************************************
143 C
144 C......  The segment is a tube;  invoke the algorithm
145 C......  from "GNOTUB" inline to get "SAFER" and "SAFSEG".
146 C
147 C**********************************************************
148  
149         IF (R .LT. PT2) SAFER = PT2 - R
150         IF (R .GT. PT3) SAFER = R - PT3
151  
152         IF (IFLAG .EQ. 2)  THEN
153  
154 C********************************************************************
155 C......  Handle the case in which we have a PHI-segment of a tube.
156 C......  In addition to the radial distance (SAFER) and the
157 C......  axial distance (SAFEZ) we compute here the distance (SAFSEG)
158 C......  to the PHI-segment boundary that is closest to the point:
159 C
160 C......  SAFSEG is only calculated if PHI lies outside the interval
161 C......  [PHI1, PHI2]. Here PHI is the angle to the given point
162 C......  (thus we only consider SAFSEG if the point is outside the
163 C......  PHI-segment).
164 C
165 C......  Algorithm to find SAFSEG (same as in routine "GNTUBE"):
166 C
167 C......  For each PHI-boundary we find the distance from the given
168 C......  point to the outer (at RMAX) point of the segment boundary
169 C......  (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
170 C......  SAFSEG to be the distance to segment PHI1; else we set
171 C......  SAFSEG to be the distance to segment PHI2.
172 C********************************************************************
173  
174  
175 C......  Next eliminate those points whose angle PHI places them
176 C......  inside the given PHI-segment (IOPENP = 0).
177  
178           IF (IOPENP .EQ. 0) GO TO 100
179  
180  
181 C......  Get coordinates of outer endpoints (at RMAX) of both PHI-segments.
182  
183           XS1  = PT3 * COSPH1
184           YS1  = PT3 * SINPH1
185           XS2  = PT3 * COSPH2
186           YS2  = PT3 * SINPH2
187  
188 C......  Get distances (squared) from the given point to each endpoint.
189  
190           DISTS1 = (X(1) - XS1)**2  +  (X(2) - YS1)**2
191           DISTS2 = (X(1) - XS2)**2  +  (X(2) - YS2)**2
192  
193 C......  Get distance to that PHI-segment whose endpoint
194 C......  is closest to the given point.
195  
196           IF (DISTS1 .LE. DISTS2) THEN
197             SAFSEG = ABS(X(1) * SINPH1 - X(2) * COSPH1)
198           ELSE
199             SAFSEG = ABS(X(1) * SINPH2 - X(2) * COSPH2)
200           ENDIF
201  
202         ENDIF
203  
204         GO TO 100
205  
206  
207    50   CONTINUE
208  
209 C*********************************************************
210 C
211 C......  The segment is a cone; invoke the algorithm
212 C......  from "GNOCON" inline to get "SAFER" and "SAFSEG".
213 C
214 C*********************************************************
215  
216         ZLENI = 0.5 / DZ
217  
218         FACT1  = (PT4 - PT2) * ZLENI
219         FACT2  = (PT5 - PT3) * ZLENI
220         RIN   = PT2 + FACT1 * (DZ + XT3)
221         ROUT  = PT3 + FACT2 * (DZ + XT3)
222  
223         IF (R .LT. RIN) THEN
224           SAFER = (RIN - R) / SQRT(1.0 + FACT1 * FACT1)
225         ELSE
226           IF (R .GT. ROUT) THEN
227             SAFER = (R - ROUT) / SQRT(1.0 + FACT2 * FACT2)
228           ENDIF
229         ENDIF
230  
231  
232         IF (IFLAG .EQ. 2)  THEN
233  
234 C********************************************************************
235 C......  Handle the case in which we have a PHI-segment of a cone.
236 C......  In addition to the radial distance (SAFER) and the
237 C......  axial distance (SAFEZ) we compute here the distance (SAFSEG)
238 C......  to the PHI-segment boundary that is closest to the point:
239 C
240 C......  SAFSEG is only calculated if PHI lies outside the interval
241 C......  [PHI1, PHI2]. Here PHI is the angle to the given point
242 C......  (thus we only consider SAFSEG if the point is outside the
243 C......  PHI-segment).
244 C
245 C......  Algorithm to find SAFSEG (same as in routine "GNTUBE"):
246 C
247 C......  For each PHI-boundary we find the distance from the given
248 C......  point to the outer (at ROUT) point of the segment boundary
249 C......  (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
250 C......  SAFSEG to be the distance to segment PHI1; else we set
251 C......  SAFSEG to be the distance to segment PHI2.
252 C********************************************************************
253  
254  
255 C......  Next eliminate those points whose angle PHI places them
256 C......  inside the given PHI-segment (IOPENP = 0).
257  
258            IF (IOPENP .EQ. 0) GO TO 100
259  
260 C......  Get coordinates of outer endpoints (at ROUT) of both PHI-segments.
261  
262           IF (XT3 .LT. -DZ) THEN
263             ROUT = PT3
264           ELSEIF (XT3 .GT. DZ) THEN
265             ROUT = PT5
266           ENDIF
267  
268           XS1  = ROUT * COSPH1
269           YS1  = ROUT * SINPH1
270           XS2  = ROUT * COSPH2
271           YS2  = ROUT * SINPH2
272  
273 C......  Get distances (squared) from the given point to each endpoint.
274  
275           DISTS1 = (X(1) - XS1)**2  +  (X(2) - YS1)**2
276           DISTS2 = (X(1) - XS2)**2  +  (X(2) - YS2)**2
277  
278 C......  Obtain distance to that PHI-segment whose endpoint
279 C......  is closest to the given point.
280  
281           IF (DISTS1 .LE. DISTS2) THEN
282             SAFSEG = ABS(X(1) * SINPH1 - X(2) * COSPH1)
283           ELSE
284             SAFSEG = ABS(X(1) * SINPH2 - X(2) * COSPH2)
285           ENDIF
286  
287         ENDIF
288  
289  
290   100   CONTINUE
291         TSAFE  = MAX (SAFEZ, SAFER, SAFSEG)
292  
293         IF (TSAFE .GT. 0.0) THEN
294           IF (TSAFE .LT. SAFE) SAFE  = TSAFE
295         ENDIF
296  
297         IF (TSAFE .EQ. 0.0) THEN
298           IF (X(3) .LT. PAR(IPZ1)) GO TO 200
299         ENDIF
300   150   CONTINUE
301  
302   200   CONTINUE
303         IF (IACT .EQ. 0) GO TO 999
304         IF (IACT .EQ. 1) THEN
305           IF (SNEXT .LT. SAFE) GO TO 999
306         ENDIF
307       ENDIF
308  
309  
310 C     ------------------------------------------------
311 C     |  Compute vector-distance 'SNXT' (McPherson)  |
312 C     ------------------------------------------------
313  
314       TSNXT = BIG
315  
316       DO 210 I=1, 6
317         XT(I) = X(I)
318   210 CONTINUE
319  
320       IPZ2 = 4
321  
322       DO 300  IS=1, NSEC
323         IPZ1  = IPZ2
324         IPZ2  = IPZ1 + 3
325  
326         XT(3) = X(3) - 0.5 * (PAR(IPZ1) + PAR(IPZ2))
327  
328         PT1 = 0.5 * (PAR(IPZ2) - PAR(IPZ1))
329         IF (PT1 .LE. 0.0) GO TO 300
330  
331         IF (PAR(IPZ1+1) .NE. PAR(IPZ2+1)) GO TO 250
332         IF (PAR(IPZ1+2) .NE. PAR(IPZ2+2)) GO TO 250
333  
334 C
335 C......  This Z-section is a tube.
336 C
337         PT3 = PT1
338         PT1 = PAR(IPZ1+1)
339         PT2 = PAR(IPZ1+2)
340         PT4 = PT6
341         PT5 = PT7
342  
343         CALL GNOTUB (XT, PT, 3, IFLAG, SNEXT, TSNXT, TSAFE)
344         GO TO 280
345  
346  
347   250   CONTINUE
348 C
349 C......   This Z-section is a cone.
350 C
351         PT2 = PAR(IPZ1+1)
352         PT3 = PAR(IPZ1+2)
353         PT4 = PAR(IPZ2+1)
354         PT5 = PAR(IPZ2+2)
355  
356         CALL GNOCON (XT, PT, 3, IFLAG, SNEXT, TSNXT, TSAFE)
357  
358   280 CONTINUE
359       IF (TSNXT .LT. SNXT)  SNXT = TSNXT
360   300 CONTINUE
361  
362   999 CONTINUE
363       END