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)
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.
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)
52       IFLAG = 2
53       IF(PAR(2).EQ.360.0) IFLAG = 1
55       IF (IACT .LT. 3) THEN
57 C       -------------------------------------------------
58 C       |  Compute safety-distance 'SAFE' (P.Weidhaas)  |
59 C       -------------------------------------------------
62         SAFEZ  = 0.0
63 C
64 C......  Obtain axial distance "SAFEZ".
65 C
66         ZMIN   = PAR(4)
67         ZMAX   = PAR(3*NZ+1)
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
76 C----------------------------------------------------
77 C......  Prepare parameters for PHI-segmented cone:
78 C----------------------------------------------------
80       IF (IFLAG .EQ. 2) THEN
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)
91         SINPH1 = SIN(PHI1*DEGRAD)
92         COSPH1 = COS(PHI1*DEGRAD)
93         SINPH2 = SIN(PHI2*DEGRAD)
94         COSPH2 = COS(PHI2*DEGRAD)
97 C......  Set flag "IOPENP" if point (X(1),X(2),X(3)) lies in the
98 C......  PHI-opening.
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
110 C------------------   Start of loop over Z-sections  --------------
112         IPZ2=4
114         DO 150 IS=1,NSEC
116         IPZ1=IPZ2
117         IPZ2=IPZ1+3
119         SAFSEG = 0.0
120         SAFER  = 0.0
122         XT3=X(3) - 0.5 * (PAR(IPZ1)+PAR(IPZ2))
123         DZ   = 0.5 * (PAR(IPZ2)-PAR(IPZ1))
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
139         IF (PT2 .NE. PT4) GO TO 50
140         IF (PT3 .NE. PT5) GO TO 50
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**********************************************************
149         IF (R .LT. PT2) SAFER = PT2 - R
150         IF (R .GT. PT3) SAFER = R - PT3
152         IF (IFLAG .EQ. 2)  THEN
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********************************************************************
175 C......  Next eliminate those points whose angle PHI places them
176 C......  inside the given PHI-segment (IOPENP = 0).
178           IF (IOPENP .EQ. 0) GO TO 100
181 C......  Get coordinates of outer endpoints (at RMAX) of both PHI-segments.
183           XS1  = PT3 * COSPH1
184           YS1  = PT3 * SINPH1
185           XS2  = PT3 * COSPH2
186           YS2  = PT3 * SINPH2
188 C......  Get distances (squared) from the given point to each endpoint.
190           DISTS1 = (X(1) - XS1)**2  +  (X(2) - YS1)**2
191           DISTS2 = (X(1) - XS2)**2  +  (X(2) - YS2)**2
193 C......  Get distance to that PHI-segment whose endpoint
194 C......  is closest to the given point.
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
202         ENDIF
204         GO TO 100
207    50   CONTINUE
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*********************************************************
216         ZLENI = 0.5 / DZ
218         FACT1  = (PT4 - PT2) * ZLENI
219         FACT2  = (PT5 - PT3) * ZLENI
220         RIN   = PT2 + FACT1 * (DZ + XT3)
221         ROUT  = PT3 + FACT2 * (DZ + XT3)
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
232         IF (IFLAG .EQ. 2)  THEN
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********************************************************************
255 C......  Next eliminate those points whose angle PHI places them
256 C......  inside the given PHI-segment (IOPENP = 0).
258            IF (IOPENP .EQ. 0) GO TO 100
260 C......  Get coordinates of outer endpoints (at ROUT) of both PHI-segments.
262           IF (XT3 .LT. -DZ) THEN
263             ROUT = PT3
264           ELSEIF (XT3 .GT. DZ) THEN
265             ROUT = PT5
266           ENDIF
268           XS1  = ROUT * COSPH1
269           YS1  = ROUT * SINPH1
270           XS2  = ROUT * COSPH2
271           YS2  = ROUT * SINPH2
273 C......  Get distances (squared) from the given point to each endpoint.
275           DISTS1 = (X(1) - XS1)**2  +  (X(2) - YS1)**2
276           DISTS2 = (X(1) - XS2)**2  +  (X(2) - YS2)**2
278 C......  Obtain distance to that PHI-segment whose endpoint
279 C......  is closest to the given point.
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
287         ENDIF
290   100   CONTINUE
291         TSAFE  = MAX (SAFEZ, SAFER, SAFSEG)
293         IF (TSAFE .GT. 0.0) THEN
294           IF (TSAFE .LT. SAFE) SAFE  = TSAFE
295         ENDIF
297         IF (TSAFE .EQ. 0.0) THEN
298           IF (X(3) .LT. PAR(IPZ1)) GO TO 200
299         ENDIF
300   150   CONTINUE
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
310 C     ------------------------------------------------
311 C     |  Compute vector-distance 'SNXT' (McPherson)  |
312 C     ------------------------------------------------
314       TSNXT = BIG
316       DO 210 I=1, 6
317         XT(I) = X(I)
318   210 CONTINUE
320       IPZ2 = 4
322       DO 300  IS=1, NSEC
323         IPZ1  = IPZ2
324         IPZ2  = IPZ1 + 3
326         XT(3) = X(3) - 0.5 * (PAR(IPZ1) + PAR(IPZ2))
328         PT1 = 0.5 * (PAR(IPZ2) - PAR(IPZ1))
329         IF (PT1 .LE. 0.0) GO TO 300
331         IF (PAR(IPZ1+1) .NE. PAR(IPZ2+1)) GO TO 250
332         IF (PAR(IPZ1+2) .NE. PAR(IPZ2+2)) GO TO 250
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
343         CALL GNOTUB (XT, PT, 3, IFLAG, SNEXT, TSNXT, TSAFE)
344         GO TO 280
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)
356         CALL GNOCON (XT, PT, 3, IFLAG, SNEXT, TSNXT, TSAFE)
358   280 CONTINUE
359       IF (TSNXT .LT. SNXT)  SNXT = TSNXT
360   300 CONTINUE
362   999 CONTINUE
363       END