GEANT321/cgpack/cghpla.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:19:44  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 *FCA :          05/01/99  09:39:51  by  Federico Carminati
  12 *               Converted the routine to double precision and added
  13 *               protections against divide by 0
  14 *-- Author :
  15       SUBROUTINE CGHPLA(IFACE,XYZ,ABCD)
  16 ************************************************************************
  17 *                                                                      *
  18 *     Name: CGHPLA                                                     *
  19 *     Author: E. Chernyaev                       Date:    08.08.88     *
  20 *                                                Revised:              *
  21 *                                                                      *
  22 *     Function: Compute face plane equation coefficients:              *
  23 *                    Ax + By + Cz + D = 0                              *
  24 *                                                                      *
  25 *     References: none                                                 *
  26 *                                                                      *
  27 *     Input: IFACE(*) - face                                           *
  28 *            XYZ(3,*) - node coordinates                               *
  29 *                                                                      *
  30 *     Output: ABCD(4) - plane equation coefficients                    *
  31 *                                                                      *
  32 *     Errors: none                                                     *
  33 *                                                                      *
  34 ************************************************************************
  35 #if defined(CERNLIB_DOUBLE)
  36       IMPLICIT        DOUBLE PRECISION (A-H,O-Z)
  37 #endif
  38       PARAMETER       (ONE=1,HALF=ONE/2,SMALL=1D-12)
  39       DIMENSION       GRAV(3)
  40       REAL            XYZ(3,*),ABCD(4)
  41 *SG
  42       INTEGER         IFACE(*)
  43 *SG
  44 *-
  45       NEDGE = IFACE(1)
  46       A      = 0
  47       B      = 0
  48       C      = 0
  49       GRAV(1)= 0
  50       GRAV(2)= 0
  51       GRAV(3)= 0
  52       JF     = 2
  53       DO 100 NE=1,NEDGE
  54         N1     = IFACE(JF)
  55         N2     = IFACE(JF+1)
  56         JF     = JF + 2
  57         A      = A + XYZ(2,N1)*XYZ(3,N2) - XYZ(2,N2)*XYZ(3,N1)
  58         B      = B + XYZ(3,N1)*XYZ(1,N2) - XYZ(3,N2)*XYZ(1,N1)
  59         C      = C + XYZ(1,N1)*XYZ(2,N2) - XYZ(1,N2)*XYZ(2,N1)
  60         GRAV(1)= GRAV(1) + XYZ(1,N1) + XYZ(1,N2)
  61         GRAV(2)= GRAV(2) + XYZ(2,N1) + XYZ(2,N2)
  62         GRAV(3)= GRAV(3) + XYZ(3,N1) + XYZ(3,N2)
  63   100   CONTINUE
  64       HNGRAV = HALF/NEDGE
  65       GRAV(1)= GRAV(1) * HNGRAV
  66       GRAV(2)= GRAV(2) * HNGRAV
  67       GRAV(3)= GRAV(3) * HNGRAV
  68       IF (ABS(A) .LT. SMALL) A=0
  69       IF (ABS(B) .LT. SMALL) B=0
  70       IF (ABS(C) .LT. SMALL) C=0
  71       AREAI  = ONE/MAX(SQRT(A*A + B*B + C*C),SMALL)
  72       ABCD(1)= A * AREAI
  73       ABCD(2)= B * AREAI
  74       ABCD(3)= C * AREAI
  75       ABCD(4)=-(A*GRAV(1) + B*GRAV(2) + C*GRAV(3)) * AREAI
  76       END