GEANT321/gtrak/gtnext.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:44  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 #if !defined(CERNLIB_OLD)
  11 *CMZ :  3.21/04 21/03/95  16.13.08  by  S.Giani
  12 *-- Author :
  13       SUBROUTINE GTNEXT
  14 C.
  15 C.    ******************************************************************
  16 C.    *                                                                *
  17 C.    *   SUBR. GTNEXT                                                 *
  18 C.    *                                                                *
  19 C.    *   Computes SAFETY and, only when new SAFETY is smaller than    *
  20 C.    *    STEP, computes SNEXT.                                       *
  21 C.    *   STEP has to be preset to BIG or to physical step size        *
  22 C.    *                                                                *
  23 C.    *   Called by : GTELEC, GTGAMA, GTHADR, GTMUON, GTNEUT, GTNINO   *
  24 C.    *                                                                *
  25 C.    *   Author   : S.Giani (1993)                                    *
  26 C.    *                                                                *
  27 C.    *   This routine is now based on the new 'virtual divisions'     *
  28 C.    *    algorithm to speed up the tracking.                         *
  29 C.    *   The tracking for MANY volumes is not anymore based on a step *
  30 C.    *    search: it is now based on a search through the list of     *
  31 C.    *    'possible overlapping volumes' built by GTMEDI.             *
  32 C.    *    Boolean operations and divisions along arbitrary axis are   *
  33 C.    *     now supported.                                             *
  34 C.    *                                                                *
  35 C.    ******************************************************************
  36 C.
  37 #include "geant321/gcbank.inc"
  38 #include "geant321/gcflag.inc"
  39 #include "geant321/gconsp.inc"
  40 #include "geant321/gcstak.inc"
  41 #include "geant321/gctmed.inc"
  42 #include "geant321/gctrak.inc"
  43 #include "geant321/gcvolu.inc"
  44 #include "geant321/gcshno.inc"
  45 #if defined(CERNLIB_USRJMP)
  46 #include "geant321/gcjump.inc"
  47 #endif
  48 #include "geant321/gchvir.inc"
  49 #include "geant321/gcvdma.inc"
  50       DIMENSION NUMTMP(15),NAMTMP(15)
  51 C.
  52       PARAMETER (BIG1=0.9*BIG)
  53 C.
  54       CHARACTER*4 NAME
  55       dimension iarrin(500),cxm(3),xxm(6)
  56       REAL      X0(6), XC(6), XT(6)
  57       INTEGER   IDTYP(3,12)
  58       LOGICAL   BTEST
  59 C.
  60       DATA  IDTYP / 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 2, 3, 1,
  61      +              2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 4, 3, 1, 1, 1,
  62      +              2, 3, 1, 2, 3, 1/
  63 C.
  64 C.    ------------------------------------------------------------------
  65 *
  66 * * *** Transform current point and direction into local reference system
  67 *
  68       mycoun=0
  69       myinfr=0
  70       newfl=0
  71       manyfl=0
  72       tsafet=big
  73       tsnext=big
  74 401   IF (GRMAT(10,NLEVEL).EQ.0.) THEN
  75          XC(1) = VECT(1) - GTRAN(1,NLEVEL)
  76          XC(2) = VECT(2) - GTRAN(2,NLEVEL)
  77          XC(3) = VECT(3) - GTRAN(3,NLEVEL)
  78          XC(4) = VECT(4)
  79          XC(5) = VECT(5)
  80          XC(6) = VECT(6)
  81       ELSE
  82 C*****  Code Expanded From Routine:  GTRNSF
  83 C
  84 *
  85          XL1 = VECT(1) - GTRAN(1,NLEVEL)
  86          XL2 = VECT(2) - GTRAN(2,NLEVEL)
  87          XL3 = VECT(3) - GTRAN(3,NLEVEL)
  88          XC(1) = XL1*GRMAT(1,NLEVEL) + XL2*GRMAT(2,NLEVEL) + XL3*
  89      1      GRMAT(3,NLEVEL)
  90          XC(2) = XL1*GRMAT(4,NLEVEL) + XL2*GRMAT(5,NLEVEL) + XL3*
  91      1      GRMAT(6,NLEVEL)
  92          XC(3) = XL1*GRMAT(7,NLEVEL) + XL2*GRMAT(8,NLEVEL) + XL3*
  93      1      GRMAT(9,NLEVEL)
  94 *
  95 C*****  End of Code Expanded From Routine:  GTRNSF
  96 C*****  Code Expanded From Routine:  GROT
  97 C
  98          XC(4) = VECT(4)*GRMAT(1,NLEVEL) + VECT(5)*GRMAT(2,NLEVEL) +
  99      1      VECT(6)*GRMAT(3,NLEVEL)
 100          XC(5) = VECT(4)*GRMAT(4,NLEVEL) + VECT(5)*GRMAT(5,NLEVEL) +
 101      1      VECT(6)*GRMAT(6,NLEVEL)
 102          XC(6) = VECT(4)*GRMAT(7,NLEVEL) + VECT(5)*GRMAT(8,NLEVEL) +
 103      1      VECT(6)*GRMAT(9,NLEVEL)
 104 *
 105 C*****  End of Code Expanded From Routine:  GROT
 106       ENDIF
 107 *
 108 * *** Compute distance to boundaries
 109 *
 110       SNEXT  = STEP
 111       SAFETY = BIG
 112       INGOTO = 0
 113       JVO    = LQ(JVOLUM-LVOLUM(NLEVEL))
 114       ISH    = Q(JVO+2)
 115       IF (Q(JVO+3).EQ.0.) GO TO 300
 116       if(raytra.eq.1..and.imyse.eq.1)then
 117             CALL UHTOC(NAMES(NLEVEL),4,NAME,4)
 118             CALL GFIND(NAME,'SEEN',ISSEEN)
 119             if(isseen.eq.-2.or.isseen.eq.-1)goto 300
 120       endif
 121       NIN = Q(JVO+3)
 122       IF (NIN.LT.0) GO TO 200
 123 *
 124 * *** Case with contents positioned
 125 *
 126       sneold=SNEXT
 127       nnn=0
 128       nflag=0
 129       mmm=0
 130       snxtot=0.
 131  111  if(nin.gt.1)then
 132         if(nnn.gt.0)goto 112
 133         clmoth=q(jvirt+4*(LVOLUM(NLEVEL)-1)+3)
 134         chmoth=q(jvirt+4*(LVOLUM(NLEVEL)-1)+4)
 135         ndivto=q(jvirt+4*(LVOLUM(NLEVEL)-1)+2)
 136         iaxis =q(jvirt+4*(LVOLUM(NLEVEL)-1)+1)
 137         if(iaxis.eq.4)then
 138          do 1 i=1,6
 139           xxm(i)=xc(i)
 140  1       continue
 141         endif
 142         divthi=(chmoth-clmoth)/ndivto
 143         if(iaxis.le.3)then
 144           cx=xc(iaxis)
 145           if(xc(iaxis+3).ge.0.)then
 146             inc=1
 147           else
 148             inc=-1
 149           endif
 150           xvdiv=((cx-clmoth)*ndivto/(chmoth-clmoth))+1
 151           ivdiv=xvdiv
 152           if((xvdiv-ivdiv).lt.0.0001.and.inc.eq.-1)ivdiv=ivdiv-1
 153           if(ivdiv.lt.1)then
 154             ivdiv=1
 155           elseif(ivdiv.gt.ndivto)then
 156             ivdiv=ndivto
 157           endif
 158         else
 159           call gfcoor(xc,iaxis,cx)
 160           if(iaxis.eq.4)then
 161             dr= xc(1)*xc(4)+xc(2)*xc(5)
 162 *            if(dr.eq.0.)print *,'dr.eq.0.'
 163             if(dr.ge.0.)then
 164               inc=1
 165             else
 166               inc=-1
 167             endif
 168           elseif(iaxis.eq.6)then
 169             if((cx-clmoth).lt.-1.)then
 170               cx=cx+360.
 171             elseif((cx-chmoth).gt.1.)then
 172               cx=cx-360.
 173             endif
 174             if(cx.gt.chmoth)then
 175               cx=chmoth
 176             elseif(cx.lt.clmoth)then
 177               cx=clmoth
 178             endif
 179             dfi=xc(1)*xc(5)-xc(2)*xc(4)
 180             if(dfi.ge.0)then
 181               inc=1
 182             else
 183               inc=-1
 184             endif
 185           endif
 186           xvdiv=((cx-clmoth)*ndivto/(chmoth-clmoth))+1
 187           ivdiv=xvdiv
 188           if((xvdiv-ivdiv).lt.0.0001.and.inc.eq.-1)ivdiv=ivdiv-1
 189           if(ivdiv.lt.1)then
 190             ivdiv=1
 191           elseif(ivdiv.gt.ndivto)then
 192             ivdiv=ndivto
 193           endif
 194         endif
 195         jvdiv=lq(jvirt-LVOLUM(NLEVEL))
 196  112    iofset=iq(jvdiv+ivdiv)
 197         jcont2=jvdiv+iofset+1
 198         ncont=iq(jcont2)
 199         if(ncont.eq.0)then
 200           idmi=iq(jcont2+1)
 201           idma=iq(jcont2+2)
 202           llflag=0
 203         elseif(ncont.eq.1)then
 204           idmi=iq(jcont2+2)
 205           idma=iq(jcont2+3)
 206           in=iq(jcont2+1)
 207         else
 208           idmi=iq(jcont2+ncont+1)
 209           idma=iq(jcont2+ncont+2)
 210           iii=1
 211           in=iq(jcont2+iii)
 212         endif
 213         if(nnn.eq.0)then
 214          cxold=cx
 215          if(inc.gt.0)then
 216           cmin=clmoth+(idmi-1)*(chmoth-clmoth)/ndivto
 217           if(iaxis.ne.6)then
 218            safety=min(safety,(cxold-cmin))
 219           else
 220            safefi=min(90.,(cxold-cmin))
 221            saferr=sqrt(xc(1)**2+xc(2)**2)
 222            safe22=saferr*sin(safefi)
 223            safety=min(safety,safe22)
 224           endif
 225          else
 226           cmax=(clmoth+(idma-1)*(chmoth-clmoth)/ndivto)+divthi
 227           if(iaxis.ne.6)then
 228            safety=min(safety,(cmax-cxold))
 229           else
 230            safefi=min(90.,(cmax-cxold))
 231            saferr=sqrt(xc(1)**2+xc(2)**2)
 232            safe22=saferr*sin(safefi)
 233            safety=min(safety,safe22)
 234           endif
 235          endif
 236         endif
 237         if(ncont.eq.0)goto 181
 238       elseif(nin.eq.1)then
 239         in=1
 240       endif
 241 *
 242   150 if(nin.gt.1.and.ncont.gt.1)then
 243         in=iq(jcont2+iii)
 244       endif
 245       if(nin.gt.0)then
 246 *        if(infrom.gt.0.and.myinfr.eq.0.and.newfl.eq.0)then
 247 *          if(in.eq.infrom)goto 171
 248 *        endif
 249         jin=lq(jvo-in)
 250         if(.NOT.BTEST(iq(jin),4))then
 251         else
 252           goto 171
 253         endif
 254       endif
 255       if(nin.gt.1)then
 256         llflag=0
 257         if(mmm.le.500)then
 258          do 151 ll=1,mmm
 259           if(iarrin(ll).eq.in)then
 260             llflag=1
 261             goto 171
 262           endif
 263  151     continue
 264         endif
 265         if(llflag.eq.0)then
 266           mmm=mmm+1
 267           if(mmm.le.500)then
 268             iarrin(mmm)=in
 269           endif
 270         endif
 271       endif
 272       IF (IN.LT.0)  GO TO 300
 273       JIN   = LQ(JVO-IN)
 274       IVOT  = Q(JIN+2)
 275       JVOT  = LQ(JVOLUM-IVOT)
 276       IROTT = Q(JIN+4)
 277 *
 278       IF (BTEST(IQ(JVOT),1)) THEN
 279 *       (case with JVOLUM structure locally developed)
 280          JPAR = LQ(LQ(JVOLUM-LVOLUM(NLDEV(NLEVEL))))
 281          DO 169 ILEV = NLDEV(NLEVEL), NLEVEL
 282             IF (IQ(JPAR+1).EQ.0) THEN
 283                IF (ILEV.EQ.NLEVEL) THEN
 284                   JPAR = LQ(JPAR-IN)
 285                ELSE
 286                   JPAR = LQ(JPAR-LINDEX(ILEV+1))
 287                ENDIF
 288                IF (JPAR.EQ.0) GO TO 170
 289             ELSE IF (IQ(JPAR-3).GT.1) THEN
 290                JPAR = LQ(JPAR-LINDEX(ILEV+1))
 291             ELSE
 292                JPAR = LQ(JPAR-1)
 293             ENDIF
 294   169    CONTINUE
 295          JPAR = JPAR + 5
 296          NPAR = IQ(JPAR)
 297          GO TO 179
 298       ENDIF
 299 *     (normal case)
 300   170 NPAR = Q(JVOT+5)
 301       IF (NPAR.EQ.0) THEN
 302          JPAR = JIN +9
 303          NPAR = Q(JPAR)
 304       ELSE
 305          JPAR = JVOT +6
 306       ENDIF
 307  179  if((nin.eq.1).or.(nin.gt.1.and.llflag.eq.0))then
 308 *
 309 *   * Compute distance to boundary of current content
 310 *
 311 C*****  Code Expanded From Routine:  GITRAN
 312   180 IF (IROTT .EQ. 0) THEN
 313          XT(1) = XC(1) - Q(5+JIN)
 314          XT(2) = XC(2) - Q(6+JIN)
 315          XT(3) = XC(3) - Q(7+JIN)
 316 *
 317          XT(4) = XC(4)
 318          XT(5) = XC(5)
 319          XT(6) = XC(6)
 320 *
 321       ELSE
 322          XL1 = XC(1) - Q(5+JIN)
 323          XL2 = XC(2) - Q(6+JIN)
 324          XL3 = XC(3) - Q(7+JIN)
 325          JR = LQ(JROTM-IROTT)
 326          XT(1) = XL1*Q(JR+1) + XL2*Q(JR+2) + XL3*Q(JR+3)
 327          XT(2) = XL1*Q(JR+4) + XL2*Q(JR+5) + XL3*Q(JR+6)
 328          XT(3) = XL1*Q(JR+7) + XL2*Q(JR+8) + XL3*Q(JR+9)
 329 *
 330 C*****  End of Code Expanded From Routine:  GITRAN
 331 C*****  Code Expanded From Routine:  GRMTD
 332          XT(4)=XC(4)*Q(JR+1)+XC(5)*Q(JR+2)+XC(6)*Q(JR+3)
 333          XT(5)=XC(4)*Q(JR+4)+XC(5)*Q(JR+5)+XC(6)*Q(JR+6)
 334          XT(6)=XC(4)*Q(JR+7)+XC(5)*Q(JR+8)+XC(6)*Q(JR+9)
 335 *
 336 C*****  End of Code Expanded From Routine:  GRMTD
 337       ENDIF
 338 *
 339       IACT = 1
 340       ISHT = Q(JVOT+2)
 341       IF (ISHT.LT.5) THEN
 342          IF (ISHT.EQ.1) THEN
 343             CALL GNOBOX (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
 344          ELSE IF (ISHT.EQ.2) THEN
 345             CALL GNOTRA(XT,Q(JPAR+1),IACT,1,SNEXT,SNXT,SAFE)
 346          ELSE IF (ISHT.EQ.3) THEN
 347             CALL GNOTRA(XT,Q(JPAR+1),IACT,2,SNEXT,SNXT,SAFE)
 348          ELSE
 349             CALL GNOTRP (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
 350          ENDIF
 351       ELSE IF (ISHT.LE.10) THEN
 352          IF (ISHT.EQ.5) THEN
 353             CALL GNOTUB(XT,Q(JPAR+1),IACT,1,SNEXT,SNXT,SAFE)
 354          ELSE IF (ISHT.EQ.6) THEN
 355             CALL GNOTUB(XT,Q(JPAR+1),IACT,2,SNEXT,SNXT,SAFE)
 356          ELSE IF (ISHT.EQ.7) THEN
 357             CALL GNOCON(XT,Q(JPAR+1),IACT,1,SNEXT,SNXT,SAFE)
 358          ELSE IF (ISHT.EQ.8) THEN
 359             CALL GNOCON(XT,Q(JPAR+1),IACT,2,SNEXT,SNXT,SAFE)
 360          ELSE IF (ISHT.EQ.9) THEN
 361             CALL GNOSPH (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
 362          ELSE
 363             CALL GNOPAR (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
 364          ENDIF
 365       ELSE IF (ISHT.EQ.11) THEN
 366          CALL GNOPGO (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
 367       ELSE IF (ISHT.EQ.12) THEN
 368          CALL GNOPCO (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
 369       ELSE IF (ISHT.EQ.13) THEN
 370          CALL GNOELT (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
 371       ELSE IF (ISHT.EQ.14) THEN
 372          CALL GNOHYP (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
 373       ELSE IF (ISHT.EQ.28) THEN
 374          CALL GSNGTR (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE,0)
 375       ELSE IF (ISHT.EQ.NSCTUB) THEN
 376          CALL GNOCTU (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
 377       ELSE
 378          PRINT *, ' GTNEXT : No code for shape ', ISHT
 379          STOP
 380       ENDIF
 381 *
 382       safe=max(safe,0.)
 383       if(snxt.le.-prec)snxt=big1
 384       snxt=max(snxt,0.)
 385       IF (SAFE.LT.SAFETY) SAFETY = SAFE
 386       IF (SNXT.LE.MIN(SNEXT,BIG1)) THEN
 387          INGOTO = IN
 388          SNEXT  = SNXT
 389          IGNEXT = 1
 390          LQ(JGPAR-NLEVEL-1) = JPAR
 391          IQ(JGPAR+NLEVEL+1) = NPAR
 392       ENDIF
 393       endif
 394  171  if(nin.eq.1)then
 395         goto 300
 396       elseif(nin.ge.1.and.ncont.gt.1)then
 397            iii=iii+1
 398            if(iii.le.ncont)goto 150
 399       endif
 400 *
 401 *   *         Compute distance to boundary of current volume
 402 *
 403  181  if(nnn.eq.0)then
 404                JPAR = LQ(JGPAR-NLEVEL)
 405                IACT = 2
 406                ISH  = Q(JVO+2)
 407                IF (ISH.LT.5) THEN
 408                   IF (ISH.EQ.1) THEN
 409                      CALL GNBOX (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 410                   ELSE IF (ISH.EQ.2) THEN
 411                      CALL GNTRAP (XC, Q(JPAR+1),IACT,1, SNEXT,SNXT,SAFE)
 412                   ELSE IF (ISH.EQ.3) THEN
 413                      CALL GNTRAP (XC, Q(JPAR+1),IACT,2, SNEXT,SNXT,SAFE)
 414                   ELSE
 415                      CALL GNTRP (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 416                   ENDIF
 417                ELSE IF (ISH.LE.10) THEN
 418                   IF (ISH.EQ.5) THEN
 419                      CALL GNTUBE (XC, Q(JPAR+1),IACT,1, SNEXT,SNXT,SAFE)
 420                   ELSE IF (ISH.EQ.6) THEN
 421                      CALL GNTUBE (XC, Q(JPAR+1),IACT,2, SNEXT,SNXT,SAFE)
 422                   ELSE IF (ISH.EQ.7) THEN
 423                      CALL GNCONE (XC, Q(JPAR+1),IACT,1, SNEXT,SNXT,SAFE)
 424                   ELSE IF (ISH.EQ.8) THEN
 425                      CALL GNCONE (XC, Q(JPAR+1),IACT,2, SNEXT,SNXT,SAFE)
 426                   ELSE IF (ISH.EQ.9) THEN
 427                      CALL GNSPHR (XC, Q(JPAR+1),IACT, SNEXT, SNXT, SAFE)
 428                   ELSE
 429                      CALL GNPARA (XC, Q(JPAR+1),IACT, SNEXT, SNXT, SAFE)
 430                   ENDIF
 431                ELSE IF (ISH.EQ.12) THEN
 432                   CALL GNPCON (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 433                ELSE IF (ISH.EQ.11) THEN
 434                   CALL GNPGON (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 435                ELSE IF (ISH.EQ.13) THEN
 436                   CALL GNELTU (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 437                ELSE IF (ISH.EQ.14) THEN
 438                   CALL GNHYPE (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 439                ELSE IF (ISH.EQ.28) THEN
 440                   CALL GSNGTR (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE,1)
 441                ELSE IF (ISH.EQ.NSCTUB) THEN
 442                   CALL GNCTUB (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 443                ELSE
 444                   PRINT *, ' GTNEXT : No code for shape ', ISH
 445                   STOP
 446                ENDIF
 447 *
 448                safe=max(safe,0.)
 449                if(snxt.le.-prec)snxt=big1
 450                snxt=max(snxt,0.)
 451                IF (SAFE.LT.SAFETY) SAFETY = SAFE
 452                IF (SNXT.LE.MIN(SNEXT,BIG1)) THEN
 453                   SNEXT  = SNXT
 454                   IGNEXT = 1
 455                   INGOTO = 0
 456                ENDIF
 457       endif
 458       if(iaxis.eq.4)then
 459        if(idma.eq.ndivto.and.inc.gt.0)goto 400
 460         cxm(1)=clmoth+(idmi-1)*(chmoth-clmoth)/ndivto
 461         if(idmi.eq.idma)then
 462           cxm(2)=cxm(1)+divthi
 463         else
 464           cxm(2)=(clmoth+(idma-1)*(chmoth-clmoth)/ndivto)+divthi
 465         endif
 466         cxm(3)=20000.
 467         call gntube(xxm,cxm,3,1,SNEXT,snxnew,safe)
 468         if(snxnew.lt.0.)snxnew=big1
 469         snxnew=snxnew+.004
 470         snxtot=snxtot+snxnew
 471         if(snxtot.lt.SNEXT)then
 472           xxm(1)=xxm(1)+snxnew*xxm(4)
 473           xxm(2)=xxm(2)+snxnew*xxm(5)
 474           xxm(3)=xxm(3)+snxnew*xxm(6)
 475           call gfcoor(xxm,iaxis,cxnew)
 476           xevdiv=((cxnew-clmoth)*ndivto/(chmoth-clmoth))+1
 477           ivdiv=xevdiv
 478           dr= xxm(1)*xxm(4)+xxm(2)*xxm(5)
 479 *          if(dr.eq.0.)print *,'dr.eq.0.'
 480           if(dr.ge.0.)then
 481               inc=1
 482           else
 483               inc=-1
 484           endif
 485           if((xevdiv-ivdiv).lt.0.0001.and.inc.eq.-1)ivdiv=ivdiv-1
 486           if(ivdiv.lt.1)then
 487               ivdiv=1
 488           elseif(ivdiv.gt.ndivto)then
 489               ivdiv=ndivto
 490           endif
 491           nnn=nnn+1
 492           goto 111
 493         else
 494           if(inc.gt.0)then
 495            cmax=(clmoth+(idma-1)*(chmoth-clmoth)/ndivto)+divthi
 496            safety=min(safety,(cmax-cxold))
 497           else
 498            cmin=clmoth+(idmi-1)*(chmoth-clmoth)/ndivto
 499            safety=min(safety,(cxold-cmin))
 500           endif
 501           goto 400
 502         endif
 503       endif
 504           if(nnn.ne.0.and.SNEXT.eq.sneold)goto 199
 505                x0(1) = xc(1) + SNEXT*xc(4)
 506                x0(2) = xc(2) + SNEXT*xc(5)
 507                x0(3) = xc(3) + SNEXT*xc(6)
 508                x0(4) = xc(4)
 509                x0(5) = xc(5)
 510                x0(6) = xc(6)
 511           if(iaxis.le.3)then
 512             cx=x0(iaxis)
 513             xevdiv=((cx-clmoth)*ndivto/(chmoth-clmoth))+1
 514             ievdiv=xevdiv
 515             if((xevdiv-ievdiv).lt.0.0001.and.inc.eq.-1)ievdiv=ievdiv-1
 516             if(ievdiv.lt.1)then
 517               ievdiv=1
 518             elseif(ievdiv.gt.ndivto)then
 519               ievdiv=ndivto
 520             endif
 521           else
 522             call gfcoor(x0,iaxis,cx)
 523             if(iaxis.eq.6)then
 524              if((cx-clmoth).lt.-1.)then
 525               cx=cx+360.
 526              elseif((cx-chmoth).gt.1.)then
 527               cx=cx-360.
 528              endif
 529              if(cx.gt.chmoth)then
 530               cx=chmoth
 531              elseif(cx.lt.clmoth)then
 532               cx=clmoth
 533              endif
 534             endif
 535             xevdiv=((cx-clmoth)*ndivto/(chmoth-clmoth))+1
 536             ievdiv=xevdiv
 537             if((xevdiv-ievdiv).lt.0.0001.and.inc.eq.-1)ievdiv=ievdiv-1
 538             if(ievdiv.lt.1)then
 539               ievdiv=1
 540             elseif(ievdiv.gt.ndivto)then
 541               ievdiv=ndivto
 542             endif
 543           endif
 544  199      if(ievdiv.ge.idmi.and.ievdiv.le.idma)then
 545             if(inc.gt.0)then
 546              cmax=(clmoth+(idma-1)*(chmoth-clmoth)/ndivto)+divthi
 547              if(iaxis.ne.6)then
 548               safety=min(safety,(cmax-cxold))
 549              else
 550               safefi=min(90.,(cmax-cxold))
 551               safe22=saferr*sin(safefi)
 552               safety=min(safety,safe22)
 553              endif
 554             else
 555              cmin=clmoth+(idmi-1)*(chmoth-clmoth)/ndivto
 556              if(iaxis.ne.6)then
 557               safety=min(safety,(cxold-cmin))
 558              else
 559               safefi=min(90.,(cxold-cmin))
 560               safe22=saferr*sin(safefi)
 561               safety=min(safety,safe22)
 562              endif
 563             endif
 564             goto 400
 565           endif
 566           if(iaxis.eq.6.or.iaxis.le.3)then
 567            if(ievdiv.lt.idmi.and.inc.gt.0)then
 568             if(nnn.eq.0.and.iaxis.eq.6
 569      +      .and.(chmoth-clmoth).eq.360.)nflag=1
 570             if(nflag.eq.0)then
 571 *             print *,'ievdiv=',ievdiv,' ;idmi=',idmi,' inc.gt.0'
 572 *             print *,isht,'=isht; ',iaxis,'=iaxis; ',ish,'=ish;'
 573              if(iaxis.le.3)then
 574                cmax=(clmoth+(idma-1)*(chmoth-clmoth)/ndivto)+divthi
 575                safety=min(safety,abs(cmax-cxold))
 576              elseif(iaxis.eq.6)then
 577                cmax=(clmoth+(idma-1)*(chmoth-clmoth)/ndivto)+divthi
 578                safefi=min(90.,(cmax-cxold))
 579                safe22=saferr*sin(safefi)
 580                safety=min(safety,safe22)
 581              endif
 582              goto 400
 583             endif
 584            elseif(ievdiv.gt.idma.and.inc.lt.0)then
 585             if(nnn.eq.0.and.iaxis.eq.6
 586      +      .and.(chmoth-clmoth).eq.360.)nflag=1
 587             if(nflag.eq.0)then
 588 *             print *,'ievdiv=',ievdiv,' ;idma=',idma,' inc.lt.0'
 589 *             print *,isht,'=isht; ',iaxis,'=iaxis; ',ish,'=ish;'
 590              if(iaxis.le.3)then
 591                cmin=clmoth+(idmi-1)*(chmoth-clmoth)/ndivto
 592                safety=min(safety,abs(cxold-cmin))
 593              elseif(iaxis.eq.6)then
 594                cmin=clmoth+(idmi-1)*(chmoth-clmoth)/ndivto
 595                safefi=min(90.,(cxold-cmin))
 596                safe22=saferr*sin(safefi)
 597                safety=min(safety,safe22)
 598              endif
 599              goto 400
 600             endif
 601            endif
 602           endif
 603           nnn=nnn+1
 604           sneold=SNEXT
 605           if(inc.gt.0)then
 606             if(iaxis.eq.6)then
 607              if(idma.eq.ndivto.and.(chmoth-clmoth).eq.360.)then
 608                ivdiv=1
 609              else
 610                ivdiv=idma+1
 611              endif
 612             else
 613              ivdiv=idma+1
 614             endif
 615           else
 616             if(iaxis.eq.6)then
 617              if(idmi.eq.1.and.(chmoth-clmoth).eq.360.)then
 618                ivdiv=ndivto
 619              else
 620                ivdiv=idmi-1
 621              endif
 622             else
 623              ivdiv=idmi-1
 624             endif
 625           endif
 626           goto 111
 627 *
 628 * ***    Case of volume incompletely divided
 629 *
 630   200 JDIV   = LQ(JVO-1)
 631       IAXIS  = Q(JDIV+1)
 632       IVOT   = Q(JDIV+2)
 633       JVOT   = LQ(JVOLUM-IVOT)
 634       ISHT   = Q(JVOT+2)
 635 *
 636 *  ** Get the division parameters
 637 *
 638       IF (NLEVEL.LT.NLDEV(NLEVEL)) THEN
 639          JPARM = 0
 640       ELSE
 641 *        (case with JVOLUM structure locally developed)
 642          JPARM = LQ(LQ(JVOLUM-LVOLUM(NLDEV(NLEVEL))))
 643          IF (NLEVEL.EQ.NLDEV(NLEVEL)) GO TO 215
 644          DO 210 ILEV = NLDEV(NLEVEL), NLEVEL-1
 645             IF (IQ(JPARM+1).EQ.0) THEN
 646                JPARM = LQ(JPARM-LINDEX(ILEV+1))
 647                IF (JPARM.EQ.0) GO TO 215
 648             ELSE IF (IQ(JPARM-3).GT.1) THEN
 649                JPARM = LQ(JPARM-LINDEX(ILEV+1))
 650             ELSE
 651                JPARM = LQ(JPARM-1)
 652             ENDIF
 653             IF (ILEV.EQ.NLEVEL-1) THEN
 654                NDIV = IQ(JPARM+1)
 655                ORIG =  Q(JPARM+2)
 656                SDIV =  Q(JPARM+3)
 657             ENDIF
 658   210    CONTINUE
 659          GO TO 220
 660       ENDIF
 661 *     (normal case)
 662   215 NDIV = Q(JDIV+3)
 663       ORIG = Q(JDIV+4)
 664       SDIV = Q(JDIV+5)
 665 *
 666 *  ** Look at the first and the last divisions only
 667 *
 668   220 IDT  = IDTYP(IAXIS, ISH)
 669       IF (IDT.EQ.1) THEN
 670          IN2 = 0
 671          IF (XC(IAXIS).LT.ORIG) THEN
 672             IN  = 1
 673          ELSE
 674             IN  = NDIV
 675          ENDIF
 676       ELSE IF (IDT.EQ.2) THEN
 677          R   = XC(1)**2 + XC(2)**2
 678          IF (ISH.EQ.9) R = R + XC(3)**2
 679          R   = SQRT(R)
 680          IN2 = 0
 681          IF (ISH.EQ.5.OR.ISH.EQ.6.OR.ISH.EQ.9) THEN
 682             IF (R.LT.ORIG) THEN
 683                IN  = 1
 684             ELSE
 685                IN  = NDIV
 686             ENDIF
 687          ELSE
 688 **          PRINT *, ' GTNEXT : Partially divided ',ISH,IAXIS
 689             IN  = 1
 690             IF (NDIV.GT.1) IN2 = NDIV
 691          ENDIF
 692       ELSE IF (IDT.EQ.4) THEN
 693          IN2 = 0
 694          RXY = XC(1)**2 + XC(2)**2
 695          RXY = SQRT(RXY)
 696          IF (XC(3).NE.0.0) THEN
 697             THET = RADDEG * ATAN (RXY/XC(3))
 698             IF (THET.LT.0.0) THET = THET + 180.0
 699          ELSE
 700             THET = 90.
 701          ENDIF
 702          IF (THET.LE.ORIG) THEN
 703             IN  = 1
 704          ELSE
 705             IN  = NDIV
 706          ENDIF
 707       ELSE
 708          IN2 = 0
 709          IF (ISH.EQ.5.OR.ISH.EQ.7) THEN
 710             IN  = 1
 711             IF (NDIV.GT.1) IN2 = NDIV
 712          ELSE
 713             IF (XC(1).NE.0.0.OR.XC(2).NE.0.0) THEN
 714                PHI = RADDEG * ATAN2 (XC(2), XC(1))
 715             ELSE
 716                PHI = 0.0
 717             ENDIF
 718             IF (ISH.EQ.6.OR.ISH.EQ.8) THEN
 719                IF (PHI.LT.ORIG) THEN
 720                   IN  = 1
 721                ELSE
 722                   IN  = NDIV
 723                ENDIF
 724             ELSE
 725                IN  = 1
 726                IF (NDIV.GT.1) IN2 = NDIV
 727             ENDIF
 728          ENDIF
 729       ENDIF
 730 *
 731   225 IF (IDT.EQ.1) THEN
 732          X0(1) = 0.0
 733          X0(2) = 0.0
 734          X0(3) = 0.0
 735          X0(IAXIS) = ORIG + (IN - 0.5) * SDIV
 736          IF (ISH.EQ.4.OR.(ISH.EQ.10.AND.IAXIS.NE.1)) THEN
 737             CALL GCENT (IAXIS, X0)
 738          ENDIF
 739          XT(1) = XC(1) - X0(1)
 740          XT(2) = XC(2) - X0(2)
 741          XT(3) = XC(3) - X0(3)
 742          XT(4) = XC(4)
 743          XT(5) = XC(5)
 744          XT(6) = XC(6)
 745       ELSE IF (IDT.EQ.3) THEN
 746          PH0  = DEGRAD * (ORIG + (IN - 0.5) * SDIV)
 747          CPHR = COS(PH0)
 748          SPHR = SIN(PH0)
 749          XT(1) = XC(1)*CPHR + XC(2)*SPHR
 750          XT(2) = XC(2)*CPHR - XC(1)*SPHR
 751          XT(3) = XC(3)
 752          XT(4) = XC(4)*CPHR + XC(5)*SPHR
 753          XT(5) = XC(5)*CPHR - XC(4)*SPHR
 754          XT(6) = XC(6)
 755       ELSE
 756          DO 234 I = 1, 6, 2
 757             XT(I) = XC(I)
 758             XT(I+1) = XC(I+1)
 759   234    CONTINUE
 760       ENDIF
 761 *
 762       IF (JPARM.NE.0) THEN
 763          IF (IQ(JPARM-3).GT.1) THEN
 764             JPAR = LQ(JPARM-IN)
 765          ELSE
 766             JPAR = LQ(JPARM-1)
 767          ENDIF
 768          JPAR = JPAR + 5
 769       ELSE
 770          JPAR = JVOT + 6
 771       ENDIF
 772 *
 773       IACT = 1
 774       IF (ISHT.LT.5) THEN
 775          IF (ISHT.EQ.1) THEN
 776             CALL GNOBOX (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 777          ELSE IF (ISHT.EQ.2) THEN
 778             CALL GNOTRA (XT, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
 779          ELSE IF (ISHT.EQ.3) THEN
 780             CALL GNOTRA (XT, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
 781          ELSE
 782             CALL GNOTRP (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 783          ENDIF
 784       ELSE IF (ISHT.LE.10) THEN
 785          IF (ISHT.EQ.5) THEN
 786             CALL GNOTUB (XT, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
 787          ELSE IF (ISHT.EQ.6) THEN
 788             CALL GNOTUB (XT, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
 789          ELSE IF (ISHT.EQ.7) THEN
 790             CALL GNOCON (XT, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
 791          ELSE IF (ISHT.EQ.8) THEN
 792             CALL GNOCON (XT, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
 793          ELSE IF (ISHT.EQ.9) THEN
 794             CALL GNOSPH (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 795          ELSE
 796             CALL GNOPAR (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 797          ENDIF
 798       ELSE IF (ISHT.EQ.11) THEN
 799          CALL GNOPGO (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 800       ELSE IF (ISHT.EQ.12) THEN
 801          CALL GNOPCO (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 802       ELSE IF (ISHT.EQ.13) THEN
 803          CALL GNOELT (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 804       ELSE IF (ISHT.EQ.28) THEN
 805          CALL GSNGTR (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE,0)
 806       ELSE IF (ISHT.EQ.NSCTUB) THEN
 807          CALL GNOCTU (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 808       ELSE
 809          PRINT *, ' GTNEXT : No code for shape ', ISHT
 810          STOP
 811       ENDIF
 812 *
 813       safe=max(safe,0.)
 814       if(snxt.le.-prec)snxt=big1
 815       snxt=max(snxt,0.)
 816       IF (SAFE.LT.SAFETY) SAFETY = SAFE
 817       IF (SNXT.LE.MIN(SNEXT,BIG1)) THEN
 818          SNEXT  = SNXT
 819          IGNEXT = 1
 820          if(raytra.eq.1.)ingoto=-1
 821       ENDIF
 822 *
 823       IF (IN2.NE.0) THEN
 824          IF (IN2.NE.IN) THEN
 825             IN  = IN2
 826             GO TO 225
 827          ENDIF
 828       ENDIF
 829 *       (later, this section only for concave volumes if INGOTO >0
 830   300 IACT = 1
 831       IF (IGNEXT.NE.0) THEN
 832          IF (.NOT.BTEST(IQ(JVO),2)) IACT = 0
 833       ENDIF
 834       if(nin.eq.1.and.ignext.ne.0)then
 835         if(q(jin+8).eq.0.)iact=1
 836       endif
 837       JPAR = LQ(JGPAR-NLEVEL)
 838       IF (ISH.LT.5) THEN
 839          IF (ISH.EQ.1) THEN
 840             CALL GNBOX (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE )
 841          ELSE IF (ISH.EQ.2) THEN
 842             CALL GNTRAP (XC, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
 843          ELSE IF (ISH.EQ.3) THEN
 844             CALL GNTRAP (XC, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
 845          ELSE
 846             CALL GNTRP (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 847          ENDIF
 848       ELSE IF (ISH.LE.10) THEN
 849          IF (ISH.EQ.5) THEN
 850             CALL GNTUBE (XC, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
 851          ELSE IF (ISH.EQ.6) THEN
 852             CALL GNTUBE (XC, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
 853          ELSE IF (ISH.EQ.7) THEN
 854             CALL GNCONE (XC, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
 855          ELSE IF (ISH.EQ.8) THEN
 856             CALL GNCONE (XC, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
 857          ELSE IF (ISH.EQ.9) THEN
 858             CALL GNSPHR (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 859          ELSE
 860             CALL GNPARA (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 861          ENDIF
 862       ELSE IF (ISH.EQ.12) THEN
 863          CALL GNPCON (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 864       ELSE IF (ISH.EQ.11) THEN
 865          CALL GNPGON (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 866       ELSE IF (ISH.EQ.13) THEN
 867          CALL GNELTU (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 868       ELSE IF (ISH.EQ.14) THEN
 869          CALL GNHYPE (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 870       ELSE IF (ISH.EQ.28) THEN
 871          CALL GSNGTR (XC,Q(JPAR+1),  IACT, SNEXT, SNXT, SAFE,1)
 872       ELSE IF (ISH.EQ.NSCTUB) THEN
 873          CALL GNCTUB (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
 874       ELSE
 875          PRINT *, ' GTNEXT : No code for shape ', ISH
 876          STOP
 877       ENDIF
 878 *
 879       safe=max(safe,0.)
 880       if(snxt.le.-prec)snxt=big1
 881       snxt=max(snxt,0.)
 882       IF (SAFE.LT.SAFETY) SAFETY = SAFE
 883       IF (SNXT.LE.MIN(SNEXT,BIG1)) THEN
 884          SNEXT  = SNXT
 885          IGNEXT = 1
 886          INGOTO = 0
 887       ENDIF
 888 *
 889  400  if(iswit(9).eq.123456789.and.Q(JVO+3).gt.1.)then
 890         print *,'n. of checked objects = ',mmm
 891       endif
 892       if(myinfr.gt.0)then
 893         jin=lq(jvo-myinfr)
 894         iq(jin)=ibclr(iq(jin),4)
 895         myinfr=0
 896       endif
 897       if(gonly(nlevel).eq.0..or.nvmany.ne.0) THEN
 898          if(safety.lt.tsafet)tsafet=safety
 899          if(snext.lt.tsnext)then
 900           mycoun=mycoun+1
 901           tsnext=snext
 902           tignex=ignext
 903           tingot=ingoto
 904           call gscvol
 905           if(ingoto.gt.0)then
 906             iq(jgpar2+nlevel+1)=iq(jgpar+nlevel+1)
 907             lq(jgpar2-nlevel-1)=lq(jgpar-nlevel-1)
 908           endif
 909          endif
 910          if(gonly(nlevel).eq.0.)then
 911  404       continue
 912            if(gonly(nlevel-1).eq.0..or.newfl.eq.0)then
 913              if(gonly(nlevel-1).ne.0.)newfl=1
 914              nlevel=nlevel-1
 915              jvo=lq(jvolum-lvolum(nlevel))
 916              nin=q(jvo+3)
 917              if(nin.lt.0)goto 404
 918              myinfr=lindex(nlevel+1)
 919              jin=lq(jvo-myinfr)
 920              iq(jin)=ibset(iq(jin),4)
 921              ignext=0
 922              goto 401
 923            endif
 924          endif
 925  403   continue
 926        if(manyfl.lt.nvmany)then
 927          manyfl=manyfl+1
 928          if(manyfl.eq.nfmany)goto 403
 929          levtmp=manyle(manyfl)
 930          do 402 i=1,levtmp
 931           namtmp(i)=manyna(manyfl,i)
 932           numtmp(i)=manynu(manyfl,i)
 933  402     continue
 934          call glvolu(levtmp,namtmp,numtmp,ier)
 935          if(ier.ne.0)print *,'Fatal error in GLVOLU'
 936          ignext=0
 937          goto 401
 938        endif
 939        if(tsafet.le.safety)safety=tsafet
 940        if(tsnext.le.snext)then
 941          snext=tsnext
 942          ignext=tignex
 943          ingoto=tingot
 944          call gfcvol
 945          nlevin=nlevel
 946          if(ingoto.gt.0)then
 947           iq(jgpar+nlevel+1)=iq(jgpar2+nlevel+1)
 948           lq(jgpar-nlevel-1)=lq(jgpar2-nlevel-1)
 949          endif
 950        endif
 951       endif
 952 *
 953 * *** Attempt to rescue negative SNXT due to rounding errors
 954 *
 955   900 IF (SNXT.EQ.BIG1) THEN
 956 CCC debug
 957          IF (ISWIT(9).EQ.123456789) THEN
 958             PRINT *,' GTNEXT : SNEXT,SAFETY,INGOTO=',SNEXT,SAFETY,INGOTO
 959             CALL GPCXYZ
 960          ENDIF
 961 CCC
 962          SAFETY = 0.
 963          SNEXT  = 0.
 964          IGNEXT = 1
 965          INGOTO = 0
 966       ENDIF
 967       IF(JGSTAT.NE.0) CALL GFSTAT(3)
 968 *                                                             END GTNEXT
 969       END
 970 #endif