GEANT321/ggeom/gnosph.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:53  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 GNOSPH (X, PAR, IACT, SNEXT, SNXT, SAFE)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'SPHE' VOLUME,  *
  17 C.    *       FROM OUTSIDE POINT X(1-3) ALONG DIRECTION X(4-6)SPHERE   *
  18 C.    *                                                                *
  19 C.    *       PAR   (input)  : volume parameters                       *
  20 C.    *       IACT  (input)  : action flag                             *
  21 C.    *         = 0  Compute SAFE only                                 *
  22 C.    *         = 1  Compute SAFE, and SNXT only if SNEXT .GT.new SAFE *
  23 C.    *         = 2  Compute both SAFE and SNXT                        *
  24 C.    *         = 3  Compute SNXT only                                 *
  25 C.    *       SNEXT (input)  : see IACT = 1                            *
  26 C.    *       SNXT  (output) : distance to volume boundary             *
  27 C.    *       SAFE  (output) : shortest distance to any boundary       *
  28 C.    *                                                                *
  29 C.    *    ==>Called by : GNEXT, GTNEXT                                *
  30 C.    *         Author  A.McPherson,  P.Weidhaas  *********            *
  31 C.    *                                                                *
  32 C.    ******************************************************************
  33 C.
  34 #if !defined(CERNLIB_SINGLE)
  35       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  36 #endif
  37 #include "geant321/gconsp.inc"
  38       REAL X(6),PAR(6),SNEXT,SNXT,SAFE
  39       DIMENSION IS(4),SS(4)
  40 C.
  41 C.    ----------------------------------------------------------------
  42 C.
  43
  44       SNXT = BIG
  45       R2  = X(1)*X(1) + X(2)*X(2) + X(3)*X(3)
  46       R   = SQRT (R2)
  47
  48       IF (IACT .LT. 3) THEN
  49
  50 C       -------------------------------------------------
  51 C       |  Compute safety-distance 'SAFE' (P.Weidhaas)  |
  52 C       -------------------------------------------------
  53
  54         RIN  = PAR(1)
  55         ROUT = PAR(2)
  56         IF (R .LT. RIN) THEN
  57           SAFE  = RIN - R
  58         ELSEIF (R .GT. ROUT) THEN
  59           SAFE  = R - ROUT
  60         ENDIF
  61
  62
  63         IF (IACT .EQ. 0) GO TO 999
  64         IF (IACT .EQ. 1) THEN
  65           IF (SNEXT .LT. SAFE) GO TO 999
  66         ENDIF
  67       ENDIF
  68
  69
  70 C     ------------------------------------------------
  71 C     |  Compute vector-distance 'SNXT' (McPherson)  |
  72 C     ------------------------------------------------
  73
  74
  75       BA=X(1)*X(4)+X(2)*X(5)+X(3)*X(6)
  76       IF(R2.GE.PAR(2)*PAR(2).AND.BA.GE.0.0) GO TO 999
  77 C
  78       CA=R2-PAR(2)*PAR(2)
  79       DISC=BA*BA-CA
  80       IF(DISC.LE.0.0) GO TO 999
  81 C
  82       RDISC=SQRT(DISC)
  83       SMAX=-BA+RDISC
  84       SMIN=-BA-RDISC
  85 C
  86 C                   NOW DO RMIN
  87 C
  88       CA=R2-PAR(1)*PAR(1)
  89       DISC=BA*BA-CA
  90 C
  91       SMIN1=SMIN
  92       SMAX1=-1.0
  93       SMIN2=SMIN
  94       SMAX2=SMAX
  95 C
  96       IF(DISC.LE.0.0) GO TO 30
  97       RDISC=SQRT(DISC)
  98       SMIN2=-BA+RDISC
  99       SMAX1=-BA-RDISC
 100 C
 101    30 CONTINUE
 102 C
 103 C                   NOW DO THE PHI STUFF.
 104 C
 105       IP2=0
 106       SMNP1=0.0
 107       SMXP1=SMAX2
 108 C
 109       IF(PAR(6)-PAR(5).GE.360.0) GO TO  110
 110 C
 111       DPSGN=X(1)*X(5)-X(2)*X(4)
 112       PHL=PAR(5)/RADDEG
 113       PHH=PAR(6)/RADDEG
 114       ISMIN=0
 115       ISMAX=0
 116 C
 117       TSGN=SIN(PHL)
 118       TCSG=COS(PHL)
 119       DEN=X(4)*TSGN-X(5)*TCSG
 120       IF(DEN.EQ.0.0) GO TO 40
 121       SNL=(X(2)*TCSG-X(1)*TSGN)/DEN
 122       IF(ABS(TSGN).GT.1.E-6.AND.(X(2)+SNL*X(5))*TSGN.LT.0.) GO TO 40
 123       IF(ABS(TCSG).GT.1.E-6.AND.(X(1)+SNL*X(4))*TCSG.LT.0.) GO TO 40
 124 C
 125       ISMIN=1
 126       SMIN=SNL
 127       IF(DPSGN.GT.0.0) GO TO 40
 128       ISMIN=0
 129       ISMAX=1
 130       SMAX=SNL
 131 C
 132    40 CONTINUE
 133 C
 134       TSGN=SIN(PHH)
 135       TCSG=COS(PHH)
 136       DEN=X(4)*TSGN-X(5)*TCSG
 137       IF(DEN.EQ.0.0) GO TO 60
 138       SNH=(X(2)*TCSG-X(1)*TSGN)/DEN
 139       IF(ABS(TSGN).GT.1.E-6.AND.(X(2)+SNH*X(5))*TSGN.LT.0.) GO TO 60
 140       IF(ABS(TCSG).GT.1.E-6.AND.(X(1)+SNH*X(4))*TCSG.LT.0.) GO TO 60
 141       IF(DPSGN.LT.0.0) GO TO 50
 142       ISMAX=1
 143       SMAX=SNH
 144       GO TO 60
 145 C
 146    50 CONTINUE
 147       ISMIN=1
 148       SMIN=SNH
 149 C
 150    60 CONTINUE
 151 C
 152       IF(ISMIN.EQ.0.OR.ISMAX.EQ.0) GO TO 80
 153       IF(SMAX.LT.0.0.AND.SMAX.GT.SMIN) GO TO 999
 154       IF(SMIN.LT.0.0) SMIN=0.0
 155       IF(SMAX.LT.0.0) GO TO  100
 156       IF(SMAX.GT.SMIN) GO TO 70
 157 C
 158 C                         SMAX +VE AND LESS THAN SMIN - 2 ALLOWED
 159 C                         REGIONS.
 160 C
 161       IP2=1
 162       SMXP1=SMAX
 163       SMNP2=SMIN
 164       SMXP2=SMAX2
 165       GO TO  110
 166 C
 167    70 CONTINUE
 168 C
 169 C                         SMIN +VE AND SMAX GT SMIN: NORMAL SINGLE
 170 C                         REGION
 171 C
 172       SMNP1=SMIN
 173       SMXP1=SMAX
 174       GO TO  110
 175 C
 176    80 CONTINUE
 177       IF(ISMIN.EQ.1) GO TO  100
 178       IF(ISMAX.EQ.0) GO TO  90
 179 C
 180 C                    SMAX BUT NO SMIN
 181 C
 182       SMXP1=SMAX
 183       GO TO  110
 184 C
 185    90 CONTINUE
 186 C
 187 C                   NO SMIN OR SMAX: ALWAYS IN OR ALWAYS OUT.
 188 C
 189       DPH=PAR(5)-PAR(4)
 190       IF(DPH.LT.180.0.AND.DPH.GT.0.0) GO TO 999
 191       IF(DPH.LT.-180.0) GO TO 999
 192       GO TO  110
 193 C
 194   100 CONTINUE
 195 C
 196 C                  SMIN BUT NO SMAX.
 197 C
 198       SMNP1=SMIN
 199 C
 200   110 CONTINUE
 201 C
 202 C                  NOW DO THETA.
 203 C
 204       IT2=0
 205       SMNT1=0.0
 206       SMXT1=SMAX2
 207       IF(PAR(4)-PAR(3).GE.180.0) GO TO 360
 208 C
 209       TH=PAR(3)
 210       IT=1
 211       ITLN=0
 212       ITLX=0
 213       ITHN=0
 214       ITHX=0
 215 C
 216   120 CONTINUE
 217 C
 218       IF(TH.NE.90.0) GO TO 130
 219       IF(X(6).EQ.0.0) GO TO 220
 220 C
 221       ST=-X(3)/X(6)
 222       STST=-X(6)
 223       GO TO 180
 224 C
 225   130 CONTINUE
 226 C
 227       TT=TAN(TH/RADDEG)
 228       TT2=TT*TT
 229 C
 230       A=X(4)*X(4)+X(5)*X(5)-TT2*X(6)*X(6)
 231       B=X(1)*X(4)+X(2)*X(5)-TT2*X(3)*X(6)
 232       C=X(1)*X(1)+X(2)*X(2)-TT2*X(3)*X(3)
 233 C
 234       IF(A.NE.0.0) GO TO 140
 235       IF(B.EQ.0.0) GO TO 220
 236 C
 237       ST=-C*0.5/B
 238 C
 239       Z=X(3)+ST*X(6)
 240       IF(Z*TT.LT.0.0) GO TO 220
 241 C
 242       STST=(B+ST*A)/Z
 243       ITRY=2
 244 C
 245       GO TO 180
 246   140 CONTINUE
 247 C
 248       BA=B/A
 249       CA=C/A
 250       DISC=BA*BA-CA
 251       IF(DISC.LT.0.0) GO TO 220
 252 C
 253       RDISC=0.0
 254       IF(DISC.GT.0.0) RDISC=SQRT(DISC)
 255       ITRY=1
 256       ST=-BA-RDISC
 257 C
 258   150 CONTINUE
 259 C
 260       IF(ST.LT.0.0) GO TO 160
 261       Z=X(3)+ST*X(6)
 262       IF(Z.EQ.0.0.AND.ABS(A).LT.0.0) GO TO 170
 263       IF(RDISC.EQ.0.0) GO TO 160
 264       IF(Z*TT.LT.0.0) GO TO 160
 265 C
 266       STST=(B+ST*A)/Z
 267       GO TO 180
 268 C
 269   160 CONTINUE
 270 C
 271       IF(ITRY.EQ.2) GO TO 220
 272       ST=-BA+RDISC
 273       ITRY=2
 274       GO TO 150
 275 C
 276   170 CONTINUE
 277       STST=-X(6)
 278   180 CONTINUE
 279 C
 280       IF(IT.NE.1) GO TO 200
 281       IF(STST.GT.0.0) GO TO 190
 282       ITLX=1
 283       SMXTL=ST
 284       GO TO 160
 285 C
 286   190 CONTINUE
 287       ITLN=1
 288       SMNTL=ST
 289       GO TO 160
 290 C
 291   200 CONTINUE
 292       IF(STST.GT.0.0) GO TO 210
 293       ITHN=1
 294       SMNTH=ST
 295       GO TO 160
 296 C
 297   210 CONTINUE
 298       ITHX=1
 299       SMXTH=ST
 300       GO TO 160
 301 C
 302   220 CONTINUE
 303       IF(IT.EQ.2) GO TO 230
 304       IT=2
 305       TH=PAR(4)
 306       GO TO 120
 307 C
 308   230 CONTINUE
 309 C
 310 C              ORDER THE VARIOUS BOUNDARIES.
 311 C
 312       ICOUNT=0
 313       IF(ITLN.EQ.0.OR.SMNTL.LE.0.0) GO TO 240
 314       IS(1)=1
 315       SS(1)=SMNTL
 316       ICOUNT=1
 317 C
 318   240 CONTINUE
 319       IF(ITLX.EQ.0.OR.SMXTL.LE.0.0) GO TO 260
 320       IPL=ICOUNT+1
 321       IF(ICOUNT.EQ.0.OR.SMXTL.GT.SS(1)) GO TO 250
 322       IS(2)=IS(1)
 323       SS(2)=SS(1)
 324       IPL=1
 325   250 CONTINUE
 326       ICOUNT=ICOUNT+1
 327       IS(IPL)=2
 328       SS(IPL)=SMXTL
 329 C
 330   260 CONTINUE
 331       IST=3
 332       IF(ITHN.EQ.0.OR.SMNTH.LE.0.0) GO TO 320
 333       STEST=SMNTH
 334 C
 335   270 CONTINUE
 336       IPL=ICOUNT+1
 337       IF(ICOUNT.EQ.0) GO TO 310
 338       DO 280 IC=1,ICOUNT
 339       IC1=ICOUNT-IC+1
 340       IF(STEST.GT.SS(IC1)) GO TO 290
 341       IPL=IPL-1
 342   280 CONTINUE
 343 C
 344   290 CONTINUE
 345       IF(IPL.EQ.ICOUNT+1) GO TO 310
 346       IM=ICOUNT+1-IPL
 347       DO 300 I=1,IM
 348       I1=ICOUNT-I+1
 349       I2=I1+1
 350       IS(I2)=IS(I1)
 351       SS(I2)=SS(I1)
 352   300 CONTINUE
 353 C
 354   310 CONTINUE
 355       ICOUNT=ICOUNT+1
 356       IS(IPL)=IST
 357       SS(IPL)=STEST
 358 C
 359   320 CONTINUE
 360       IF(IST.EQ.4) GO TO 330
 361       IF(ITHX.EQ.0.OR.SMXTH.LE.0.0) GO TO 330
 362       IST=4
 363       STEST=SMXTH
 364       GO TO 270
 365 C
 366   330 CONTINUE
 367 C
 368 C               CHECK WHETHER 1ST IS MAX OR MIN.
 369 C
 370       IF(ICOUNT.EQ.0) GO TO 350
 371       IF(IS(1).EQ.2.OR.IS(1).EQ.4) GO TO 340
 372 C
 373 C               START WITH MIN.
 374 C
 375       SMNT1=SS(1)
 376       IF(ICOUNT.GE.2) SMXT1=SS(2)
 377       IF(ICOUNT.LE.2) GO TO 360
 378       IT2=1
 379       SMNT2=SS(3)
 380       SMXT2=SMAX2
 381       IF(ICOUNT.GE.4) SMXT2=SS(4)
 382       GO TO 360
 383 C
 384   340 CONTINUE
 385 C
 386 C              START WITH MAX SO 1ST MIN IS 0.0
 387 C
 388       SMNT1=0.0
 389       SMXT1=SS(1)
 390       IF(ICOUNT.LE.1) GO TO 360
 391       IT2=1
 392       SMNT2=SS(2)
 393       SMXT2=SMAX2
 394       IF(ICOUNT.GE.3) SMXT2=SS(3)
 395       GO TO 360
 396 C
 397   350 CONTINUE
 398 C
 399 C              NO INTERSECTIONS ALWAYS IN OR ALWAYS OUT.
 400 C
 401       R=X(1)*X(1)+X(2)*X(2)
 402       IF(R.GT.0.0) R=SQRT(R)
 403       TH=90.0
 404       IF(X(3).NE.0.0) TH=ATAN(R/X(3))*RADDEG
 405       IF(TH.LT.0.0) TH=180.0+TH
 406       IF(TH.LT.PAR(3).OR.TH.GT.PAR(4)) GO TO 999
 407 C
 408   360 CONTINUE
 409 C
 410 C              NOW FIND SMALLEST S ALOWED BY ALL.
 411 C
 412       IF(SMAX1.LE.SMIN1) GO TO 370
 413       SMAXR=SMAX1
 414       SMINR=SMIN1
 415       IRT=1
 416       GO TO 380
 417 C
 418   370 CONTINUE
 419       SMAXR=SMAX2
 420       SMINR=SMIN2
 421       IRT=2
 422 C
 423   380 CONTINUE
 424       IF(SMNP1.GT.SMAXR) GO TO 430
 425       IF(SMXP1.LT.SMINR) GO TO 390
 426       SMIN=SMINR
 427       SMAX=SMAXR
 428       IF(SMNP1.GT.SMIN) SMIN=SMNP1
 429       IF(SMXP1.LT.SMAX) SMAX=SMXP1
 430       IPT=1
 431       GO TO 400
 432 C
 433   390 CONTINUE
 434       IF(IP2.EQ.0) GO TO 430
 435       IF(SMNP2.GT.SMAXR) GO TO 430
 436       IF(SMXP2.LT.SMINR) GO TO 430
 437       SMIN=SMINR
 438       SMAX=SMAXR
 439       IF(SMNP2.GT.SMIN) SMIN=SMNP2
 440       IF(SMXP2.LT.SMAX) SMAX=SMXP2
 441       IPT=2
 442 C
 443   400 CONTINUE
 444 C
 445       IF(SMNT1.GT.SMAX) GO TO 420
 446       IF(SMXT1.LT.SMIN) GO TO 410
 447       IF(SMNT1.GT.SMIN) SMIN=SMNT1
 448       GO TO 440
 449 C
 450   410 CONTINUE
 451       IF(IT2.EQ.0) GO TO 420
 452       IF(SMNT2.GT.SMAX) GO TO 420
 453       IF(SMXT2.LT.SMIN) GO TO 420
 454       IF(SMNT2.GT.SMIN) SMIN=SMNT2
 455       GO TO 440
 456 C
 457   420 CONTINUE
 458       IF(IPT.EQ.1) GO TO 390
 459   430 CONTINUE
 460       IF(IRT.EQ.1) GO TO 370
 461       GO TO 999
 462 C
 463   440 CONTINUE
 464       IF(SMIN.LE.0.)GO TO 999
 465       SNXT = SMIN
 466
 467   999 CONTINUE
 468       END