GEANT321/ggeom/gnsphr.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.3  1996/03/28 08:50:23  cernlib
   6 * In the call to function MAX  use NULL istead of 0. to get the proper
   7 * type of argument.
   8 *
   9 * Revision 1.2  1996/02/27 10:12:06  ravndal
  10 * Precision problem (neg. Sqrt) solved
  11 *
  12 * Revision 1.1.1.1  1995/10/24 10:20:54  cernlib
  13 * Geant
  14 *
  15 *
  16 #include "geant321/pilot.h"
  17 *CMZ :  3.21/02 29/03/94  15.41.30  by  S.Giani
  18 *-- Author :
  19       SUBROUTINE GNSPHR (X, PAR, IACT, SNEXT, SNXT, SAFE)
  20 C.
  21 C.    ******************************************************************
  22 C.    *                                                                *
  23 C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'SPHE' VOLUME,  *
  24 C.    *       FROM INSIDE POINT X(1-3) ALONG DIRECTION X(4-6)          *
  25 C.    *                                                                *
  26 C.    *       PAR   (input)  : volume parameters                       *
  27 C.    *       IACT  (input)  : action flag                             *
  28 C.    *         = 0  Compute SAFE only                                 *
  29 C.    *         = 1  Compute SAFE, and SNXT only if SNEXT .GT.new SAFE *
  30 C.    *         = 2  Compute both SAFE and SNXT                        *
  31 C.    *         = 3  Compute SNXT only                                 *
  32 C.    *       SNEXT (input)  : see IACT = 1                            *
  33 C.    *       SNXT  (output) : distance to volume boundary             *
  34 C.    *       SAFE  (output) : shortest distance to any boundary       *
  35 C.    *                                                                *
  36 C.    *    ==>Called by : GNEXT, GTNEXT                                *
  37 C.    *         Author  A.McPherson,  P.Weidhaas  *********            *
  38 C.    *                                                                *
  39 C.    ******************************************************************
  40 C.
  41 #if !defined(CERNLIB_SINGLE)
  42       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  43 #endif
  44 #include "geant321/gconsp.inc"
  45       REAL X(6), PAR(6), SNXT, SNEXT, SAFE
  46       PARAMETER (ZERO=0,ONE=1,SMALL=ONE/BIG)
  47
  48 C---------------------------------------------------------------
  49
  50       SNXT = BIG
  51       R2   = X(1)*X(1) + X(2)*X(2) + X(3)*X(3)
  52       R    = SQRT (R2)
  53       RIN  = PAR(1)
  54       ROUT = PAR(2)
  55       SAF1 = R - RIN
  56       SAF2 = ROUT - R
  57
  58       IF (IACT .LT. 3) THEN
  59
  60 C       -------------------------------------------------
  61 C       |  Compute safety-distance 'SAFE' (P.Weidhaas)  |
  62 C       -------------------------------------------------
  63
  64          SAFE = MIN (SAF1, SAF2)
  65          IF( RIN .EQ. 0.0 ) THEN
  66             SAFE = ROUT-R
  67             SAF3 = SAFE
  68             SAF4 = SAFE
  69          ELSE
  70             SAF3 = R
  71             SAF4 = R
  72          ENDIF
  73          IF( R .GT. 1.E-5) THEN
  74             IF( PAR(4)-PAR(3).GE.180.) THEN
  75                SAF3 = BIG
  76             ELSE
  77                TH = ACOS(X(3)/SNGL(R))*RADDEG
  78                IF (TH.LT.0) TH = 180.+TH
  79                DTH1 = TH-PAR(3)
  80                DTH2 = PAR(4)-TH
  81                DTH = MIN(DTH1,DTH2)
  82                SAF3 = R*SIN(DTH*DEGRAD)
  83             ENDIF
  84             RXY2 = X(1)*X(1)+X(2)*X(2)
  85             IF( RXY2 .GT. 1.E-12.AND. PAR(6)-PAR(5) .LT. 360) THEN
  86                RXY = SQRT(RXY2)
  87                PHI=ATAN2(X(2),X(1))*RADDEG
  88                DPH = PHI-PAR(5)
  89                SG = SIGN(ONE,DPH)
  90                DPH = MOD( ABS(DPH), ONE*360 )
  91                IF(SG.LE.0.) DPH = 360.-DPH
  92                DPHT = PAR(6)-PAR(5)
  93                IF(DPHT .LT. 0.) DPHT = DPHT+360.
  94                DDPH = MIN(DPH,DPHT-DPH)
  95                IF( DDPH .GT. 90.) THEN
  96                   SAF4 = BIG
  97                ELSE
  98                   SAF4 = RXY*SIN(DDPH*DEGRAD)
  99                ENDIF
 100             ENDIF
 101          ENDIF
 102          SAFE = MIN(ONE*SAFE,SAF3,SAF4)
 103
 104          IF (IACT .EQ. 0) GO TO 999
 105          IF (IACT .EQ. 1) THEN
 106             IF (SNEXT .LT. SAFE) GO TO 999
 107          ENDIF
 108       ENDIF
 109
 110 C     ------------------------------------------------
 111 C     |  Compute vector-distance 'SNXT' (McPherson)  |
 112 C     ------------------------------------------------
 113
 114       IF(R.LT.1.0E-5) GO TO 70
 115
 116       C1=(X(1)*X(4)+X(2)*X(5)+X(3)*X(6)) / R
 117       RTMP = PAR(2)
 118       SGN=ONE
 119       AC=R*R * (ONE-C1*C1)
 120       IF(AC.GT.PAR(1)**2.OR.C1.GT.0.0) GO TO 10
 121       SGN=-ONE
 122       RTMP = PAR(1)
 123    10 CONTINUE
 124       RTMPAC=MAX(RTMP*RTMP-AC,ZERO)
 125       SNXT=SQRT(RTMPAC)*SGN-R*C1
 126       DPSGN=X(1)*X(5)-X(2)*X(4)
 127       PHI2=PAR(6)
 128       IF(DPSGN.LT.0.0) PHI2=PAR(5)
 129 C
 130       TSGN=SIN(PHI2*DEGRAD)
 131       TCSG=COS(PHI2*DEGRAD)
 132       DEN=X(4)*TSGN-X(5)*TCSG
 133       IF(DEN.EQ.0.0) GO TO 20
 134       SP=(X(2)*TCSG-X(1)*TSGN)/DEN
 135       IF(SP.LT.0.0) GO TO 20
 136       IF(ABS(TSGN).GT.1.E-6.AND.(X(2)+SP*X(5))*TSGN
 137      +       .LT.0.)GO TO 20
 138       IF(ABS(TCSG).GT.1.E-6.AND.(X(1)+SP*X(4))*TCSG
 139      +       .LT.0.)GO TO 20
 140 C
 141       IF(SP.LT.SNXT) SNXT=SP
 142 C
 143    20 CONTINUE
 144 C
 145       TH2=PAR(4)
 146       IBOUN=0
 147    30 CONTINUE
 148       IBOUN=IBOUN+1
 149 C
 150       STH=SIN(TH2*DEGRAD)
 151       CTH=COS(TH2*DEGRAD)
 152       STH2=STH*STH
 153       CTH2=CTH*CTH
 154       AA=STH2*X(6)**2-CTH2*(X(4)**2+X(5)**2)
 155       BB=STH2*X(6)*X(3)-CTH2*(X(4)*X(1)+
 156      +X(5)*X(2))
 157       CC=STH2*X(3)**2-CTH2*(X(1)**2+X(2)**2)
 158       SQ=BB*BB-AA*CC
 159       IF(SQ.LT.0.0) GO TO 60
 160       ST=BIG
 161       IF(ABS(BB).GE.SMALL) ST=-CC*0.5/BB
 162       ITRY=1
 163       IF(AA.EQ.0.0) GO TO 45
 164       T2=SQRT(SQ)/AA
 165       ITRY=0
 166       T1=-BB/AA
 167    40 CONTINUE
 168       ST=T1+T2
 169    45 CONTINUE
 170       IF((X(3)+ST*X(6))*CTH.LT.0.0) GO TO 50
 171       IF(ST.GT.0.0.AND.ST.LT.SNXT) SNXT=ST
 172    50 CONTINUE
 173       IF(ITRY.NE.0) GO TO 60
 174       ITRY=1
 175       T2=-T2
 176       GO TO 40
 177    60 CONTINUE
 178       IF(IBOUN.NE.1) GO TO 999
 179       TH2=PAR(3)
 180       GO TO 30
 181 C
 182    70 CONTINUE
 183 C              THIS BIT FOR X,Y,Z=0,0,0
 184 C
 185       SNXT=PAR(2)
 186 C              WE HAVE IGNORED THETA AND PHI BUT USERS SHOULDN'T
 187 C              USE THETA PHI SEGMENTATION AT R=0 ANYWAY.
 188 C
 189   999 CONTINUE
 190       END