This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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