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