]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/gtrak/gtnext.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Better printing for MAXSTEP
[u/mrichter/AliRoot.git] / GEANT321 / gtrak / gtnext.F
CommitLineData
fe4da5cc 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
14C.
15C. ******************************************************************
16C. * *
17C. * SUBR. GTNEXT *
18C. * *
19C. * Computes SAFETY and, only when new SAFETY is smaller than *
20C. * STEP, computes SNEXT. *
21C. * STEP has to be preset to BIG or to physical step size *
22C. * *
23C. * Called by : GTELEC, GTGAMA, GTHADR, GTMUON, GTNEUT, GTNINO *
24C. * *
25C. * Author : S.Giani (1993) *
26C. * *
27C. * This routine is now based on the new 'virtual divisions' *
28C. * algorithm to speed up the tracking. *
29C. * The tracking for MANY volumes is not anymore based on a step *
30C. * search: it is now based on a search through the list of *
31C. * 'possible overlapping volumes' built by GTMEDI. *
32C. * Boolean operations and divisions along arbitrary axis are *
33C. * now supported. *
34C. * *
35C. ******************************************************************
36C.
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)
51C.
52 PARAMETER (BIG1=0.9*BIG)
53C.
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
59C.
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/
63C.
64C. ------------------------------------------------------------------
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
74401 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
82C***** Code Expanded From Routine: GTRNSF
83C
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*
95C***** End of Code Expanded From Routine: GTRNSF
96C***** Code Expanded From Routine: GROT
97C
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*
105C***** 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*
311C***** 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*
330C***** End of Code Expanded From Routine: GITRAN
331C***** 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*
336C***** 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
956CCC debug
957 IF (ISWIT(9).EQ.123456789) THEN
958 PRINT *,' GTNEXT : SNEXT,SAFETY,INGOTO=',SNEXT,SAFETY,INGOTO
959 CALL GPCXYZ
960 ENDIF
961CCC
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