]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/ggeom/gnext.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Minor corrections after big transformer changes
[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gnext.F
CommitLineData
fe4da5cc 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)
14C.
15C. ******************************************************************
16C. * *
17C. * SUBR. GNEXT (X, SNEXT, SAFETY) *
18C. * *
19C. * Computes SNEXT and SAFETY *
20C. * SNEXT (output) : distance to closest boundary *
21C. * from point X(1-3) along X(4-6) *
22C. * SAFETY (output) : shortest distance to any boundary *
23C. * *
24C. * Called by : User *
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. * Important notice: in case of MANY volumes, it might happen *
35C. * that at the end of GNEXT the volume where X was found to *
36C. * be is not anymore the current volume in GCVOLU (being the *
37C. * tree ready for the volume where the particles is entering).*
38C. * When this happens, the variable MYCOUN in the common block *
39C. * PHYCOU is greater than one; the user can take the proper *
40C. * actions (before calling GINVOL, for example) like restoring*
41C. * the old tree: this is possible because GMEDIA has saved *
42C. * the needed informations in MANYLE(NFMANY) (for nlevel), in *
43C. * MANYNA(NFMANY,i) (for names(i)) and in MANYNU(NFMANY,i) *
44C. * (for number(i)); this is sufficient to restore the old tree*
45C. * by calling GLVOLU. *
46C. * *
47C. ******************************************************************
48C.
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
65C.
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/
69C.
70C. ------------------------------------------------------------------
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