5 * Revision 1.1.1.1 1995/10/24 10:20:55 cernlib
9 #include "geant321/pilot.h"
10 *CMZ : 3.21/02 29/03/94 15.41.30 by S.Giani
12 SUBROUTINE GSAGTR(X,P,SAFETY,INSIDE)
14 C. ******************************************************************
17 C. * SUBROUTINE GSAGTR(X,P,SAFETY,INSIDE) *
18 C. * Routine to cumpute the 'Safe' distance to nearest boundary *
19 C. * of the general twisted trapezoid. Point is given in X(1-3), *
20 C. * parameters of trapezoid are in P, if INSIDE = 1 point is *
21 C. * inside shape if 0 outside. 'SAFE' distance is returned in *
23 C. * I have not yet been able to come up with an exact *
24 C. * computation of this. From the outside I use an exscribed *
25 C. * cylinder with its axis as the line joining the centres of *
26 C. * the trapezia at the ends. As far as I can see there is no *
27 C. * straight forward way of finding a reliable conservative *
28 C. * estimate from the inside; so I set SAFETY to 0.0 for all *
29 C. * points inside the shape. The radius of the exscribed *
30 C. * cylinder is given by the longest of the eight distances *
31 C. * from the centre of the trapezium to each corner on each of *
33 C. * Called by : GSNGTR *
34 C. * A.C.McPherson. 23rd April 1985. *
37 C. ******************************************************************
43 C Check if point is inside.
45 IF(INSIDE.EQ.1) GO TO 999
47 C Compute radius of cylinder.
52 RC2=(P(I0)+(P(I0+2)-P(13))*P(1))**2+
53 + (P(I0+1)+(P(I0+3)-P(14))*P(1))**2
54 IF(RC2.GT.RCYL) RCYL=RC2
55 RC2=(P(I0)-(P(I0+2)-P(13))*P(1))**2+
56 + (P(I0+1)-(P(I0+3)-P(14))*P(1))**2
57 IF(RC2.GT.RCYL) RCYL=RC2
59 IF(RC2.GT.0.0) RCYL=SQRT(RC2)
61 C The direction cosines of the axis of the cylinder
62 C are computed and the distance from the point to the
63 C axis is calculated from the cross product of these
64 C direction cosines with the vector from the origin to
65 C the point. The cross product of this cross product
66 C with the direction cosines gives the vector from the
67 C axis to the point. Subtracting this times the ratio
68 C (D-RCYL)/D where D is the distance of the point from
69 C the axis and RCYL is the radius of the cylinder gives
70 C the point on the cylinder nearest to the point. If
71 C this is outside the z range of the shape then the
72 C distance along the cylinder surface to the z limit is
73 C added in quadrature.
75 IF(ABS(X(3)).GT.P(1)) SAFETY=ABS(X(3)-P(1))
76 TTH2=P(13)**2+P(14)**2
81 DX=DIR2*X(3)-DIR3*X(2)
82 DY=DIR3*X(1)-DIR1*X(3)
83 DZ=DIR1*X(2)-DIR2*X(1)
85 IF(D2.LT.RCYL*RCYL) GO TO 999
87 C Only Z component of vector is needed.
93 IF(ABS(Z).LT.P(1)) GO TO 999
94 SAFETY=SQRT(SAFETY*SAFETY+(ABS(Z)-P(1))**2/CTH2)