]>
Commit | Line | Data |
---|---|---|
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) | |
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 |