]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1995/10/24 10:20:55 cernlib | |
6 | * Geant | |
7 | * | |
8 | * | |
9 | #include "geant321/pilot.h" | |
10 | *CMZ : 3.21/02 29/03/94 15.41.30 by S.Giani | |
11 | *-- Author : | |
12 | SUBROUTINE GSAGTR(X,P,SAFETY,INSIDE) | |
13 | C. | |
14 | C. ****************************************************************** | |
15 | C. * * | |
16 | 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 * | |
22 | C. * SAFETY. * | |
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 * | |
32 | C. * the end faces. * | |
33 | C. * Called by : GSNGTR * | |
34 | C. * A.C.McPherson. 23rd April 1985. * | |
35 | C. * * | |
36 | C. * * | |
37 | C. ****************************************************************** | |
38 | C. | |
39 | DIMENSION X(3),P(30) | |
40 | C | |
41 | SAFETY=0.0 | |
42 | C | |
43 | C Check if point is inside. | |
44 | C | |
45 | IF(INSIDE.EQ.1) GO TO 999 | |
46 | C | |
47 | C Compute radius of cylinder. | |
48 | C | |
49 | RCYL=0.0 | |
50 | DO 10 I=1,4 | |
51 | I0=I*4+11 | |
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 | |
58 | 10 CONTINUE | |
59 | IF(RC2.GT.0.0) RCYL=SQRT(RC2) | |
60 | C | |
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. | |
74 | C | |
75 | IF(ABS(X(3)).GT.P(1)) SAFETY=ABS(X(3)-P(1)) | |
76 | TTH2=P(13)**2+P(14)**2 | |
77 | CTH2=1.0/(1.0+TTH2) | |
78 | DIR3=SQRT(CTH2) | |
79 | DIR1=P(13)*DIR3 | |
80 | DIR2=P(14)*DIR3 | |
81 | DX=DIR2*X(3)-DIR3*X(2) | |
82 | DY=DIR3*X(1)-DIR1*X(3) | |
83 | DZ=DIR1*X(2)-DIR2*X(1) | |
84 | D2=DX*DX+DY*DY+DZ*DZ | |
85 | IF(D2.LT.RCYL*RCYL) GO TO 999 | |
86 | C | |
87 | C Only Z component of vector is needed. | |
88 | C | |
89 | DDZ=DX*DIR2-DY*DIR1 | |
90 | D=SQRT(D2) | |
91 | Z=X(3)-DDZ*(D-RCYL)/D | |
92 | SAFETY=D-RCYL | |
93 | IF(ABS(Z).LT.P(1)) GO TO 999 | |
94 | SAFETY=SQRT(SAFETY*SAFETY+(ABS(Z)-P(1))**2/CTH2) | |
95 | 999 CONTINUE | |
96 | RETURN | |
97 | END |