GEANT321/ggeom/gsagtr.F

   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