GEANT321/ggeom/gvsafe.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:57  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.31  by  S.Giani
  11 *-- Author :
  12 *
  13       FUNCTION GVSAFE (XYZ, C1, C, NC)
  14 ************************************************************************
  15 *                                                                      *
  16 * GVSAFE Calculates the distance from a point                YB 870511 *
  17 *        to a surface, and the gradient in this point.                 *
  18 *        This is an approximation of the distance, which is smaller    *
  19 *        than the real distance.                                       *
  20 *                                                                      *
  21 * Input : XYZ     point coordinates                                    *
  22 *         C1      constant coefficient from the surface                *
  23 *         C(1)    number of non constant coefficients from             *
  24 *                       the surface                                    *
  25 *         C(2),C(3),... non constant coefficients from the surface     *
  26 *         NC      total number of coefficients from the surface        *
  27 *                                                                      *
  28 * Extra OUTPUT: COMMON/SLATE/ISLATE(2)  ISLATE(2)=ISIGN (from GVSIGN)  *
  29 *               SLATE(10)    - value of S(X)                           *
  30 *               SLATE(11:13) - gradient                                *
  31 *               SLATE(14)    - normalisation factor for distance       *
  32 *               SLATE(15)    - normalisation factor for gradient       *
  33 *                                                                      *
  34 ************************************************************************
  35       COMMON /SLATE/ SLATE (40)
  36       INTEGER        ISLATE(40)
  37       EQUIVALENCE    (SLATE,ISLATE)
  38       EQUIVALENCE    (JDUMM,SLATE(2))
  39       REAL           XYZ(3) , C(*)
  40       INTEGER        GVSIGN
  41 *-----------------------------------------------------------------------
  42       SLATE(14) = 1.
  43       SLATE(15) = 1.
  44       IF (NC.NE.2)              GO TO 200
  45 *
  46 *     case with simplified surface X=C0, Y=C0, Z=C0, X*2+Y*2=C0
  47 *     (happens only when initialisation is done)
  48 *
  49       IAX = C(2)
  50       SLATE(11) = 0.
  51       SLATE(12) = 0.
  52       SLATE(13) = 0.
  53       IF (IAX.LE.3)             THEN
  54         SLATE(10) = C1 + XYZ(IAX)
  55         SLATE(IAX+10) = 1.
  56         CC = SLATE(10)
  57       ELSE
  58         RXY2 = XYZ(1)**2 + XYZ(2)**2
  59         RXY = SQRT(RXY2)
  60         SLATE(10) = (RXY2 - C1**2)
  61         SLATE(11) = 2.*XYZ(1)
  62         SLATE(12) = 2.*XYZ(2)
  63         SLATE(15) = 2.*RXY
  64         SLATE(14) = (RXY-C1)
  65         CC = SLATE(10)/SLATE(14)
  66       ENDIF
  67       IF (SLATE(10))            101,102,103
  68   101 ISLATE(2) = -1
  69       GO TO 999
  70   102 ISLATE(2) = 0
  71       GO TO 999
  72   103 ISLATE(2) = +1
  73       GO TO 999
  74 *
  75 *     case with surfaces with 4, 7 or 10 coefficients (normal case)
  76 *
  77   200 JDUMM = GVSIGN (XYZ, C1, C, NC)
  78       CC    = SLATE(1)
  79       SLATE(10) = SLATE(1)
  80       SLATE(11) = C(2)
  81       SLATE(12) = C(3)
  82       SLATE(13) = C(4)
  83       IF (NC .EQ. 4)            GO TO 999
  84       IF (NC .EQ. 7)            THEN
  85         AA = 1.
  86       ELSE
  87 CCC     AA = C(5)**2+C(6)**2+C(7)**2CCC        AA = C(5)**2+C(6)**2+C(7)**2
  88 CCC  ++ 0.5*(C(8)**2+C(9)**2+C(10)**2)
  89 CCC        AA = SQRT(AA)
  90         AA = 2.0
  91       ENDIF
  92       CALL GVGRAD (XYZ, C, NC, SLATE(11))
  93       TT2  = SLATE(11)**2 + SLATE(12)**2 + SLATE(13)**2
  94       SLATE(14) = (TT2+4.*AA*ABS(CC))
  95       IF (SLATE(14)-TT2.LE.0.1*TT2)     THEN
  96         SLATE(15) = SQRT(SLATE(14))
  97         SLATE(14) = SLATE(15)
  98       ELSE
  99         SLATE(15) = SQRT(TT2)
 100         SLATE(14) = 0.5*(SLATE(15) + SQRT(SLATE(14)))
 101       ENDIF
 102       IF (ABS(SLATE(14)).LE.0.) SLATE(14) = 1.E-10
 103       CC = CC / SLATE(14)
 104
 105   999 GVSAFE = CC
 106
 107       END