5 * Revision 1.1.1.1 1995/10/24 10:20:53 cernlib
9 #include "geant321/pilot.h"
10 *CMZ : 3.21/02 29/03/94 15.41.30 by S.Giani
12 SUBROUTINE GNPCON (X, PAR, IACT, SNEXT, SNXT, SAFE)
14 C. ******************************************************************
16 C. * COMPUTE DISTANCE UP TO INTERSECTION WITH 'PCON' VOLUME, *
17 C. * FROM INSIDE POINT X(1-3) ALONG DIRECTION X(4-6) *
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 *
29 C. * ==>Called by : GNEXT, GTNEXT *
30 C. * Author A.McPherson, P.Weidhaas ********* *
32 C. ******************************************************************
34 #include "geant321/gconsp.inc"
35 DIMENSION X(6), PAR(9), XT(6), PT(7)
37 EQUIVALENCE (PT(1), DZ), (PT(2), PT2), (PT(3), PT3)
38 EQUIVALENCE (PT(4), PT4), (PT(5), PT5), (PT(6), PT6)
39 EQUIVALENCE (PT(7), PT7)
41 EQUIVALENCE (XT(3), XT3)
43 C. ---------------------------------------------------
47 R2 = X(1)*X(1) + X(2)*X(2)
55 SAFEZ = MIN (SAFZ1, SAFZ2)
58 C...... First determine in which z-segment the particle is located.
60 DO 10 JPH=7, 3*(NZ-1)+1, 3
61 IF (X(3) .LE. PAR(JPH)) THEN
70 C...... The particle is in the segment bounded by z-planes at
71 C...... Z1=PAR(IPL) and Z2=PAR(IPH), i.e., Z1 < X(3) < Z2.
73 C...... Set parameters for this segment and translate z-coordinate
74 C...... of point relative to center of this segment.This is done in
75 C...... preparation of invoking the algorithms used in "GNTUBE" and
76 C...... "GNCONE" (which for reasons of efficiency and clarity are
77 C...... implemented inline).
80 DZ = 0.5 * (PAR(IPH) - PAR(IPL))
87 IF (PT7 .GT. 360.0) PT7 = PT7 - 360.0
89 XT3 = X(3) - 0.5 * (PAR(IPL) + PAR(IPH))
98 IF (PAR(2) .EQ. 360.0) IND = 1
100 IF (IACT .LT. 3) THEN
102 C -------------------------------------------------
103 C | Compute safety-distance 'SAFE' (P.Weidhaas) |
104 C -------------------------------------------------
105 SAFZ1 = DZ - ABS(XT3)
106 SAFEZ = MIN (SAFEZ,SAFZ1)
109 C...... Determine whether the segment is a tube or a cone.
112 IF (PT2 .NE. PT4) GO TO 50
113 IF (PT3 .NE. PT5) GO TO 50
115 C*********************************************************
117 C...... The segment is a tube: invoke the algorithm
118 C...... from routine "GNTUBE" inline to get "SAFER".
120 C*********************************************************
124 SAFER = MIN (SAFR1, SAFR2)
126 IF (IND .EQ. 2) GO TO 70
133 C*********************************************************
135 C...... The segment is a cone: invoke the algorithm
136 C...... from routine "GNCONE" inline to get "SAFER".
138 C*********************************************************
143 C...... Compute radial distance to inner wall.
145 FACT = (PT4 - PT2) * ZLENI
146 R1 = PT2 + FACT * SAFZ2
147 SAFR1 = (R - R1) / SQRT(1.0 + FACT*FACT)
149 C...... Compute radial distance to outer wall.
151 FACT = (PT5 - PT3) * ZLENI
152 R2 = PT3 + FACT * SAFZ2
153 SAFR2 = (R2 - R) / SQRT(1.0 + FACT*FACT)
154 SAFER = MIN (SAFR1, SAFR2)
156 IF (IND .EQ. 1) GO TO 100
160 C********************************************************************
161 C...... Here we handle the case of a PHI-segment of a tube or cone.
162 C...... In addition to the radial distances (SAFR1, SAFR2) and the
163 C...... axial distances (SAFZ1, SAFZ2) we compute here the distance
164 C...... to the PHI-segment boundary that is closest to the point.
166 C...... For each PHI-boundary we find the distance from the given
167 C...... point to the outer (at R2) point of the segment boundary
168 C...... (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
169 C...... "SAFSEG" to be the distance to segment PHI1, else we set
170 C...... "SAFSEG" to be the distance to segment PHI2.
171 C*********************************************************************
175 IF (PHI2 .LT. PHI1) PHI2 = PHI2 + TWOPI
182 C...... Get coordinates of outer endpoints (at R2) of both PHI-segments.
189 C...... Get distances (squared) from given point to each endpoint.
191 DISTS1 = (X(1) - XS1)**2 + (X(2) - YS1)**2
192 DISTS2 = (X(1) - XS2)**2 + (X(2) - YS2)**2
194 C...... Get distance to that PHI-segment whose endpoint
195 C...... is closest to the given point.
197 IF (DISTS1 .LE. DISTS2) THEN
198 SAFSEG = ABS(X(1) * SINPH1 - X(2) * COSPH1)
200 SAFSEG = ABS(X(1) * SINPH2 - X(2) * COSPH2)
206 SAFE = MIN (SAFEZ, SAFER, SAFSEG)
208 IF (IACT .EQ. 0) GO TO 999
209 IF (IACT .EQ. 1) THEN
210 IF (SNEXT .LT. SAFE) GO TO 999
215 C ------------------------------------------------
216 C | Compute vector-distance 'SNXT' (McPherson) |
217 C ------------------------------------------------
219 * Avoid rounding effects induced by translation ********************
221 IF (ABS(XT(3)).GE.DZ) XT(3) = (1.-0.000001)*XT(3)
223 IF (PT2 .NE. PT4) GO TO 200
224 IF (PT3 .NE. PT5) GO TO 200
233 CALL GNTUBE (XT, PT, 3, IND, SNEXT, SNXT, TSAFE)
237 CALL GNCONE (XT, PT, 3, IND, SNEXT, SNXT, TSAFE)