GEANT321/ggeom/gsngtr.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.2  1996/02/27 09:54:45  ravndal
   6 * Double precision for some locals, thanks to G.Poulard
   7 *
   8 * Revision 1.1.1.1  1995/10/24 10:20:56  cernlib
   9 * Geant
  10 *
  11 *
  12 #include "geant321/pilot.h"
  13 *CMZ :  3.21/02 29/03/94  15.41.30  by  S.Giani
  14 *-- Author :
  15       SUBROUTINE GSNGTR(X,P,IACT,SNEXT,SNXT,SAFE,INSIDE)
  16 C.
  17 C.    ******************************************************************
  18 C.    *                                                                *
  19 C.    *                                                                *
  20 C.    *    Routine to determine the shortest distance from the point   *
  21 C.    *    X(1-3) to the boundary of the shape of type GTRA defined    *
  22 C.    *    by the parameters P along the vector X(4-6). If INSIDE is   *
  23 C.    *    1 then the point is inside the shape and the distance is    *
  24 C.    *    returned as SNEXT. If INSIDE is 0 then the point is         *
  25 C.    *    outside the shape and if the line hits the shape then       *
  26 C.    *    if the new distance is less than the                        *
  27 C.    *    old value of SNEXT the new distance is returned as SNEXT.   *
  28 C.    *                                                                *
  29 C.    *          Called by : GNEXT, GTNEXT                             *
  30 C.    *          A.C.McPherson   22nd April 1985.                      *
  31 C.    *                                                                *
  32 C.    *                                                                *
  33 C.    ******************************************************************
  34 C.
  35 #include "geant321/gconsp.inc"
  36 #include "geant321/gcunit.inc"
  37 C
  38 #if !defined(CERNLIB_SINGLE)
  39       DOUBLE PRECISION X0,Y0,DXDZ,DYDZ,A,B,C,DISC,X1,X2,X3,SN,CP,SMALL
  40 #endif
  41       PARAMETER (SMALL=1E-10)
  42 C
  43       DIMENSION X(6),P(30),SN(2,5),IOUT(5),X0(4),Y0(4),DXDZ(4),DYDZ(4)
  44 C.
  45 C.                ---------------------------------------------
  46 C.
  47 C
  48 C               Compute Safety distance
  49 C
  50       DOUBLE PRECISION SL,SL1,SM,SM1
  51       IF(IACT.LT.3) CALL GSAGTR(X,P,SAFE,INSIDE)
  52       SNXT=BIG
  53       IF (IACT .EQ. 0) GO TO 999
  54       IF (IACT .EQ. 1) THEN
  55         IF (SNEXT .LT. SAFE) GO TO 999
  56       ENDIF
  57 C
  58 C               First compute the distance along the line to the
  59 C               boundaries.
  60 C
  61 C               The distance to the planes defined by z=+/-P(1).
  62 C
  63       IF(X(6).EQ.0.0) THEN
  64           SN(1,1)=BIG
  65           SN(2,1)=BIG
  66           GOTO 10
  67       ENDIF
  68       SN(1,1)=(-P(1)-X(3))/X(6)
  69       SN(2,1)=(P(1)-X(3))/X(6)
  70       IF(X(6).GT.0.0) GO TO 10
  71       ST=SN(2,1)
  72       SN(2,1)=SN(1,1)
  73       SN(1,1)=ST
  74    10 CONTINUE
  75 C
  76 C               The distance to the remaining four surfaces.
  77 C
  78       DO 20 I=1,4
  79       X0(I)=P(I*4+11)
  80       Y0(I)=P(I*4+12)
  81       DXDZ(I)=P(I*4+13)
  82       DYDZ(I)=P(I*4+14)
  83    20 CONTINUE
  84 C
  85       DO 65 I=1,4
  86       J=I+1
  87       IF(J.EQ.5) J=1
  88 C
  89       A=(X(4)-DXDZ(I)*X(6))*(DYDZ(J)-DYDZ(I))*X(6) -
  90      +(X(5)-DYDZ(I)*X(6))*(DXDZ(J)-DXDZ(I))*X(6)
  91 C
  92       B=(X(4)-DXDZ(I)*X(6))*(Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))*X(3)) +
  93      +(X(1)-X0(I)-DXDZ(I)*X(3))*(DYDZ(J)-DYDZ(I))*X(6) -
  94      +(X(5)-DYDZ(I)*X(6))*(X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))*X(3)) -
  95      +(X(2)-Y0(I)-DYDZ(I)*X(3))*(DXDZ(J)-DXDZ(I))*X(6)
  96 C
  97       C=(X(1)-X0(I)-DXDZ(I)*X(3))*(Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))*X(3))
  98      + - (X(2)-Y0(I)-DYDZ(I)*X(3))*(X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))*X(3))
  99 C
 100       IOUT(I+1)=0
 101       IF(C.GT.0.0) IOUT(I+1)=1
 102 C
 103 C             The solutions are in the normal form:
 104 C             s = (-B+/-SQRT(B*B-4.0*A*C))*0.5/A
 105 C
 106       SN(1,I+1)=BIG
 107       SN(2,I+1)=BIG
 108       IF(ABS(A).GT.1.0E-10) GO TO 30
 109 C
 110 C             A = 0 only one solution.
 111 C
 112       IF(ABS(B).LT.1.0E-10) GO TO 60
 113 C
 114       SN(1,I+1)=-C/B
 115       GO TO 60
 116 C
 117    30 CONTINUE
 118       IF(ABS(C).GT.1.0E-10) GO TO 40
 119       SN(1,I+1)=0.0
 120       SN(2,I+1)=0.0
 121       IF(ABS(B).LT.1.0E-10) GO TO 60
 122       SN(1,I+1)=-C/B
 123       IF(C.EQ.0.0) SN(1,I+1)=SIGN(SMALL,B)
 124       SN(2,I+1)=-B/A
 125       GO TO 50
 126 C
 127    40 CONTINUE
 128       DISC=B*B-A*C*4.0
 129       IF(DISC.LT.0.0) GO TO 60
 130       IF(DISC.GT.0.0) DISC=SQRT(DISC)
 131       SN(1,I+1)=(-B-DISC)*0.5/A
 132       SN(2,I+1)=(-B+DISC)*0.5/A
 133 C
 134    50 CONTINUE
 135       IF(SN(2,I+1).GT.SN(1,I+1)) GO TO 60
 136       ST=SN(2,I+1)
 137       SN(2,I+1)=SN(1,I+1)
 138       SN(1,I+1)=ST
 139 C
 140    60 CONTINUE
 141 C
 142       DO 65 K=1,2
 143       IF(ABS(SN(K,I+1)).GT.1.0E+05.OR.ABS(SN(K,I+1)).LT.1.0E-05)
 144      +GO TO 65
 145 C
 146       X1=X(1)+SN(K,I+1)*X(4)
 147       X2=X(2)+SN(K,I+1)*X(5)
 148       X3=X(3)+SN(K,I+1)*X(6)
 149       CP=(X1-X0(I)-DXDZ(I)*X3)*(Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))*X3)
 150      + - (X2-Y0(I)-DYDZ(I)*X3)*(X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))*X3)
 151       CP=CP/SQRT((X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))*X3)**2+
 152      +   (Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))*X3)**2)
 153 C
 154       IF(ABS(CP).LT.0.0001) GO TO 65
 155       IF(ABS(CP/SN(K,I+1)).LT.1.0E-06) GO TO 65
 156       WRITE(CHMAIL,1020) I,K,SN(K,I+1)
 157       CALL GMAIL(0,0)
 158       WRITE(CHMAIL,1021) X
 159       CALL GMAIL(0,0)
 160       WRITE(CHMAIL,1022) X1,X2,X3,CP
 161       CALL GMAIL(0,0)
 162       WRITE(CHMAIL,1023) A,B,C,DISC
 163       CALL GMAIL(0,0)
 164       WRITE(CHMAIL,1024) INSIDE
 165       CALL GMAIL(0,0)
 166       WRITE(CHMAIL,1025) X0
 167       CALL GMAIL(0,0)
 168       WRITE(CHMAIL,1026) Y0
 169       CALL GMAIL(0,0)
 170       WRITE(CHMAIL,1027) DXDZ
 171       CALL GMAIL(0,0)
 172       WRITE(CHMAIL,1028) DYDZ
 173       CALL GMAIL(0,0)
 174  1020 FORMAT('0 GSNGTR ERROR - I,K =',2I2,' SN =',E13.5)
 175  1021 FORMAT(' X =',6F15.6)
 176  1022 FORMAT(' X1,X2,X3 =',3F15.6,' CP =',E15.6)
 177  1023 FORMAT(' A =',E15.6,' B =',E15.6,' C =',E15.6,' DISC =',E15.6)
 178  1024 FORMAT(' INSIDE =',I3)
 179  1025 FORMAT('   X0 =',4E15.6)
 180  1026 FORMAT('   Y0 =',4E15.6)
 181  1027 FORMAT(' DXDZ =',4E15.6)
 182  1028 FORMAT(' DYDZ =',4E15.6)
 183 C
 184    65 CONTINUE
 185 C
 186 C
 187 C             Have computed the two distances for the z planes and
 188 C             the four surfaces. Combine them accordingly as to
 189 C             whether the point is inside or outside the shape.
 190 C
 191       IF(INSIDE.EQ.0) GO TO 80
 192 C
 193 C             Point is inside shape.
 194 C
 195       DO 70 I=1,5
 196       DO 70 J=1,2
 197       IF(SN(J,I).GT.0.0.AND.SN(J,I).LT.SNXT) SNXT=SN(J,I)
 198    70 CONTINUE
 199       GO TO 999
 200 C
 201    80 CONTINUE
 202 C
 203 C             Point is outside shape.
 204 C
 205       IOUT(1)=0
 206       IF(ABS(X(3)).GT.P(1)) IOUT(1)=1
 207 C
 208 C             For each of five sets of SN and IOUT, IOUT(I) equal to 1
 209 C             indicates that the point is outside the shape by the Ith
 210 C             test, SN(1,I) is the distance to the first change in the
 211 C             test and SN(2,I) is the distance to the second change.
 212 C             The remaining logic just attempts to find a distance when
 213 C             the line is inside by all five tests, bearing in mind that
 214 C             for some tests the line can start inside, leave and return
 215 C             inside.
 216 C
 217       SL=-1.0
 218       SM=BIG
 219       SM1=BIG
 220       DO 100 I=1,5
 221       IF(IOUT(I).EQ.0) GO TO 90
 222       IF(SN(2,I).LT.0.0) GO TO 999
 223       IF(SN(1,I).LT.0.0.AND.SN(2,I).GT.SL) SL=SN(2,I)
 224       IF(SN(1,I).GT.SL) SL=SN(1,I)
 225       IF(SN(1,I).GE.0.0.AND.SN(2,I).LT.SM) SM=SN(2,I)
 226       GO TO 100
 227    90 CONTINUE
 228       IF(SN(1,I).LT.0.0.AND.SN(2,I).GE.0.0.AND.SN(2,I).LT.SM) SM=SN(2,I)
 229       IF(SN(1,I).LT.0.0.OR.SN(1,I).GT.SM1) GO TO 100
 230       IF(SN(1,I).GE.SN(2,I)) GO TO 100
 231       SM1=SN(1,I)
 232       SL1=SN(2,I)
 233   100 CONTINUE
 234 C
 235 C             SL is the largest of the five distances to the first
 236 C             time the line is inside. SM is the smallest to the
 237 C             last time the point is inside. SM1 is the smallest
 238 C             distance to when the line is temporarily outside
 239 C             one of the tests.
 240 C
 241       IF(SM.LE.SL) GO TO 999
 242       IF(SM1.GT.SL) GO TO 130
 243 C
 244   110 CONTINUE
 245 C
 246 C             In this loop SL is updated by the return after SM1
 247 C             if SM1 is less than SL.
 248 C
 249       SL=SL1
 250       IF(SM.LE.SL) GO TO 999
 251       SM1=SM
 252 C
 253       DO 120 I=1,5
 254       IF(IOUT(I).EQ.1) GO TO 120
 255       IF(SN(2,I).LE.SL.OR.SN(1,I).GT.SM1) GO TO 120
 256       IF(SN(1,I).GE.SN(2,I)) GO TO 120
 257       SM1=SN(1,I)
 258       SL1=SN(2,I)
 259   120 CONTINUE
 260 C
 261       IF(SM1.GT.SL) GO TO 130
 262 C
 263       GO TO 110
 264   130 CONTINUE
 265 C
 266       IF(SL.LT.SNXT) SNXT=SL
 267 C
 268   999 CONTINUE
 269       END