GEANT321/ggeom/gnext.F

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