GEANT321/ggeom/gintco.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:51  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.29  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GINTCO (X, RLEFT, RRIGHT, DZ, TAU, TAUL)
  13
  14 C        ********************************************
  15 C        *  This subroutine finds the intersection  *
  16 C        *  of a given ray (described by array X)   *
  17 C        *  with a given cone (described by radii   *
  18 C        *  RLEFT and RRIGHT and half-distance DZ). *
  19 C        *  Output parameter is TAU, and inter-     *
  20 C        *  section point is X = XP + TAU*XD, etc.  *
  21 C        *                                          *
  22 C        *  Called by GNOCON                        *
  23 C        *  Programmed by:  Patrick Weidhaas        *
  24 C        *          CERN,   March 1988              *
  25 C        ********************************************
  26
  27 #include "geant321/gconsp.inc"
  28       DIMENSION X(6)
  29 #if !defined(CERNLIB_SINGLE)
  30       DOUBLE PRECISION XP,YP,ZP,XD,YD,ZD,S,T,U,V,W
  31       DOUBLE PRECISION DISCR,SQDISC
  32 #endif
  33 C----------------------------------------------------
  34
  35 C...... Point of origin of ray:
  36
  37       XP = X(1)
  38       YP = X(2)
  39       ZP = X(3)
  40
  41 C...... Direction cosines:
  42
  43       XD = X(4)
  44       YD = X(5)
  45       ZD = X(6)
  46
  47       TAU  = BIG
  48       TAUL = BIG
  49
  50       S = 0.5 * (RLEFT + RRIGHT)
  51       T = (RLEFT - RRIGHT) / DZ
  52
  53 C......  Cone equation is:    x**2 + y**2 - Az**2 + Bz + C = 0
  54
  55       A = 0.25 * T*T
  56       B = S * T
  57       C = -S*S
  58
  59 C......  To obtain "TAU", we must solve the quadratic equation
  60 C......  Ut**2 + Vt + W = 0 .
  61
  62       U = XD**2 + YD**2 - A*ZD**2
  63       V = 2.0 * (XP*XD + YP*YD - A*ZP*ZD) + B*ZD
  64       W = XP**2 + YP**2 - A*ZP**2 + B*ZP + C
  65
  66       DISCR = V*V - 4.0*U*W
  67       IF (DISCR .LE. 0.0) GO TO 999
  68       IF(U.EQ.0.)GO TO 999
  69       SQDISC = SQRT (DISCR)
  70       TAU1 = (-V + SQDISC) / (2.0*U)
  71       TAU2 = (-V - SQDISC) / (2.0*U)
  72
  73
  74 C......  Set TAU to the smallest positive root;
  75 C......  otherwise let TAU = BIG .
  76 C
  77 C......  If both roots are positive, set TAUL to
  78 C......  the larger one: it may be needed in the
  79 C......  case of a PHI-segmented cone.
  80
  81       IF (TAU1 .LT. 0.0) THEN
  82         IF (TAU2 .LT. 0.0) GO TO 999
  83         TAU = TAU2
  84       ELSE
  85         TAU = TAU1
  86         IF (TAU2 .GT. 0.0) THEN
  87           TAUL = TAU2
  88           IF (TAU2.LT.TAU1) THEN
  89             TAU = TAU2
  90             TAUL = TAU1
  91           ENDIF
  92         ENDIF
  93       ENDIF
  94
  95   999 CONTINUE
  96       END
  97