5 * Revision 1.1.1.1 1995/10/24 10:20:57 cernlib
9 #include "geant321/pilot.h"
10 #if defined(CERNLIB_OLD)
11 *CMZ : 3.21/02 29/03/94 15.41.31 by S.Giani
13 SUBROUTINE GNEXT (X, SNEXT, SAFETY)
15 C. ******************************************************************
17 C. * SUBR. GNEXT (X, SNEXT, SAFETY) *
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 *
24 C. * Called by : User *
25 C. * Authors : S.Banerjee, R.Brun, F.Bruyant *
27 C. ******************************************************************
29 #include "geant321/gcbank.inc"
30 #include "geant321/gconsp.inc"
31 #include "geant321/gcshno.inc"
32 #include "geant321/gctmed.inc"
33 #include "geant321/gcvolu.inc"
34 #if defined(CERNLIB_USRJMP)
35 #include "geant321/gcjump.inc"
37 REAL X(6), X0(3), XC(6), XT(6)
42 DATA IDTYP / 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 2, 3, 1,
43 + 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 4, 3, 1, 1, 1,
46 C. ------------------------------------------------------------------
48 * *** Transform current point and direction into local reference system
50 IF (GRMAT(10,NLEVEL).EQ.0.) THEN
52 XC(I) = X(I) -GTRAN(I,NLEVEL)
56 * (later, code in line)
57 CALL GTRNSF (X, GTRAN(1,NLEVEL), GRMAT(1,NLEVEL), XC)
58 CALL GROT (X(4), GRMAT(1,NLEVEL), XC(4))
61 * *** Compute distance to boundaries
65 JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
67 IF (Q(JVO+3).EQ.0.) GO TO 300
69 IF (NIN.LT.0) GO TO 200
71 * *** Case with contents positioned
74 IF (ISEARC.GE.-1) GO TO 120
76 * ** Contents are ordered by (dynamic) GSORD, select neighbours
78 JSB = LQ(LQ(JVO-NIN-1))
83 INC = SIGN(1., XC(IAX+3))
85 CALL GFCOOR (XC, IAX, CX)
87 DR = XC(1)*XC(4) +XC(2)*XC(5)
88 IF (IAX.EQ.5) DR = DR +XC(3)*XC(6)
90 ELSE IF (IAX.EQ.6) THEN
91 INC = SIGN(1., XC(1)*XC(5)-XC(2)*XC(4))
93 INC = SIGN(1., XC(3)*(XC(1)*XC(4)+XC(2)*XC(5))
94 + -XC(6)*(XC(1)*XC(1)+XC(2)*XC(2)))
97 IDIV = LOCATF (Q(JSB+3), NSB, CX)
98 IF (IDIV.LT.0) IDIV = -IDIV
101 IF (INC.LT.0.AND.IAX.LE.3) THEN
102 SAFETY = Q(JSB+3) -CX
106 ELSE IF (IDIV.EQ.NSB) THEN
107 IF (INC.GT.0.AND.IAX.NE.7) THEN
108 SAFETY = CX -Q(JSB+2+NSB)
115 SAFETY = CX -Q(JSB+2+IDIV)
117 SAFETY = Q(JSB+3+IDIV) -CX
123 ELSE IF (IAX.EQ.6) THEN
124 IF (IDIV.EQ.0) IDIV = NSB
131 110 NCONT = IQ(JSC0+IDIV)
133 * ** Loop over (selected) contents
136 IF (IDIV.EQ.IDIVL) GO TO 400
138 * (following statement for IAX=6, when division NSB is empty)
139 IF (IDIV.GT.NSB) IDIV = 1
147 120 JNEAR = LQ(JVO-NIN-1)
148 IF (ISEARC.GT.0) THEN
149 #if !defined(CERNLIB_USRJMP)
150 CALL GUNEAR (ISEARC, 2, XC, JNEAR)
152 #if defined(CERNLIB_USRJMP)
153 CALL JUMPT4(JUNEAR,ISEARC, 2, XC, JNEAR)
155 IF (IQ(JNEAR+1).EQ.0) GO TO 300
159 IF (IQ(JNEAR+1).NE.0) THEN
165 130 IN = IQ(JNEAR+INEAR)
168 140 IN = IQ(JSCV+ICONT)
170 150 IF(IN.LE.0)GO TO 300
173 JVOT = LQ(JVOLUM-IVOT)
176 IF (NLEVEL.GE.NLDEV(NLEVEL)) THEN
177 * (case with JVOLUM structure locally developed)
178 JPAR = LQ(LQ(JVOLUM-LVOLUM(NLDEV(NLEVEL))))
179 DO 169 ILEV = NLDEV(NLEVEL), NLEVEL
180 IF (IQ(JPAR+1).EQ.0) THEN
181 IF (ILEV.EQ.NLEVEL) THEN
184 JPAR = LQ(JPAR-LINDEX(ILEV+1))
186 IF (JPAR.EQ.0) GO TO 170
187 ELSE IF (IQ(JPAR-3).GT.1) THEN
188 JPAR = LQ(JPAR-LINDEX(ILEV+1))
204 * * Compute distance to boundary of current content
206 180 IF (IROTT.EQ.0) THEN
208 XT(I) = XC(I) -Q(JIN+4+I)
212 * (later, code in line)
213 CALL GITRAN (XC, Q(JIN+5), IROTT, XT)
214 CALL GRMTD (XC(4), IROTT, XT(4))
221 CALL GNOBOX (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
222 ELSE IF (ISHT.EQ.2) THEN
223 CALL GNOTRA(XT,Q(JPAR+1),IACT,1,SNEXT,SNXT,SAFE)
224 ELSE IF (ISHT.EQ.3) THEN
225 CALL GNOTRA(XT,Q(JPAR+1),IACT,2,SNEXT,SNXT,SAFE)
227 CALL GNOTRP (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
229 ELSE IF (ISHT.LE.10) THEN
231 CALL GNOTUB(XT,Q(JPAR+1),IACT,1,SNEXT,SNXT,SAFE)
232 ELSE IF (ISHT.EQ.6) THEN
233 CALL GNOTUB(XT,Q(JPAR+1),IACT,2,SNEXT,SNXT,SAFE)
234 ELSE IF (ISHT.EQ.7) THEN
235 CALL GNOCON(XT,Q(JPAR+1),IACT,1,SNEXT,SNXT,SAFE)
236 ELSE IF (ISHT.EQ.8) THEN
237 CALL GNOCON(XT,Q(JPAR+1),IACT,2,SNEXT,SNXT,SAFE)
238 ELSE IF (ISHT.EQ.9) THEN
239 CALL GNOSPH (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
241 CALL GNOPAR (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
243 ELSE IF (ISHT.EQ.11) THEN
244 CALL GNOPGO (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
245 ELSE IF (ISHT.EQ.12) THEN
246 CALL GNOPCO (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
247 ELSE IF (ISHT.EQ.13) THEN
248 CALL GNOELT (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
249 ELSE IF (ISHT.EQ.14) THEN
250 CALL GNOHYP (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
251 ELSE IF (ISHT.EQ.28) THEN
252 CALL GSNGTR (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE,0)
253 ELSE IF (ISHT.EQ.NSCTUB) THEN
254 CALL GNOCTU (XT,Q(JPAR+1),IACT,SNEXT,SNXT,SAFE)
256 PRINT *, ' GNEXT : No code for shape ', ISHT
260 IF (SAFE.LT.SAFETY) SAFETY = SAFE
261 IF (SNXT.LT.SNEXT) THEN
263 IF (ISEARC.EQ.-2) THEN
264 IF (MOD(IQ(JSC0),2).NE.0) THEN
268 X0(I) = XC(I) + SNXT*XC(I+3)
271 IDIVB = LOCATF (Q(JSB+3), NSB, X0(IAX))
273 CALL GFCOOR (X0, IAX, CX)
274 IDIVB = LOCATF (Q(JSB+3), NSB, CX)
276 IF (IDIVB.LT.0) IDIVB = -IDIVB
283 ELSE IF (IDIVB.EQ.NSB) THEN
284 IF (IAX.NE.6) IDIVB = NSB - 1
290 IF (ISEARC.EQ.-2) THEN
291 IF (ICONT.EQ.NCONT) THEN
294 IF (IDIV.EQ.IDIVB) GO TO 300
295 IF (.NOT.BTEST(IQ(JVO),2)) THEN
301 * * Compute distance to boundary of current volume
303 JPAR = LQ(JGPAR-NLEVEL)
308 CALL GNBOX (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
309 ELSE IF (ISH.EQ.2) THEN
310 CALL GNTRAP (XC, Q(JPAR+1),IACT,1, SNEXT,SNXT,SAFE)
311 ELSE IF (ISH.EQ.3) THEN
312 CALL GNTRAP (XC, Q(JPAR+1),IACT,2, SNEXT,SNXT,SAFE)
314 CALL GNTRP (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
316 ELSE IF (ISH.LE.10) THEN
318 CALL GNTUBE (XC, Q(JPAR+1),IACT,1, SNEXT,SNXT,SAFE)
319 ELSE IF (ISH.EQ.6) THEN
320 CALL GNTUBE (XC, Q(JPAR+1),IACT,2, SNEXT,SNXT,SAFE)
321 ELSE IF (ISH.EQ.7) THEN
322 CALL GNCONE (XC, Q(JPAR+1),IACT,1, SNEXT,SNXT,SAFE)
323 ELSE IF (ISH.EQ.8) THEN
324 CALL GNCONE (XC, Q(JPAR+1),IACT,2, SNEXT,SNXT,SAFE)
325 ELSE IF (ISH.EQ.9) THEN
326 CALL GNSPHR (XC, Q(JPAR+1),IACT, SNEXT, SNXT, SAFE)
328 CALL GNPARA (XC, Q(JPAR+1),IACT, SNEXT, SNXT, SAFE)
330 ELSE IF (ISH.EQ.12) THEN
331 CALL GNPCON (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
332 ELSE IF (ISH.EQ.11) THEN
333 CALL GNPGON (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
334 ELSE IF (ISH.EQ.13) THEN
335 CALL GNELTU (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
336 ELSE IF (ISH.EQ.14) THEN
337 CALL GNHYPE (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
338 ELSE IF (ISH.EQ.28) THEN
339 CALL GSNGTR (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE,1)
340 ELSE IF (ISH.EQ.NSCTUB) THEN
341 CALL GNCTUB (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
343 PRINT *, ' GNEXT : No code for shape ', ISH
347 IF (SAFE.LT.SAFETY) SAFETY = SAFE
348 IF (SNXT.LT.SNEXT) SNEXT = SNXT
350 * * Check wether other pseudo-divisions have to be scanned
353 X0(I) = XC(I) +SNXT*XC(I+3)
356 IDIVL = LOCATF (Q(JSB+3), NSB, X0(IAX))
358 CALL GFCOOR (X0, IAX, CX)
359 IDIVL = LOCATF (Q(JSB+3), NSB, CX)
361 IF (IDIVL.LT.0) IDIVL = -IDIVL
368 ELSEIF (IDIVL.EQ.NSB)THEN
369 IF(IAX.NE.6)IDIVL=NSB-1
372 IF (IDIV.EQ.IDIVB) GO TO 400
374 193 IF ((IDIV-IDIVL)*INC.GE.0) GO TO 400
382 IF (INEAR.EQ.NNEAR) GO TO 300
387 * *** Case of volume incompletely divided
392 JVOT = LQ(JVOLUM-IVOT)
395 * ** Get the division parameters
397 IF (NLEVEL.LT.NLDEV(NLEVEL)) THEN
400 * (case with JVOLUM structure locally developed)
401 JPARM = LQ(LQ(JVOLUM-LVOLUM(NLDEV(NLEVEL))))
402 IF (NLEVEL.EQ.NLDEV(NLEVEL)) GO TO 215
403 DO 210 ILEV = NLDEV(NLEVEL), NLEVEL-1
404 IF (IQ(JPARM+1).EQ.0) THEN
405 JPARM = LQ(JPARM-LINDEX(ILEV+1))
406 IF (JPARM.EQ.0) GO TO 215
407 ELSE IF (IQ(JPARM-3).GT.1) THEN
408 JPARM = LQ(JPARM-LINDEX(ILEV+1))
412 IF (ILEV.EQ.NLEVEL-1) THEN
425 * ** Look at the first and the last divisions only
427 220 IDT = IDTYP(IAXIS, ISH)
430 IF (XC(IAXIS).LT.ORIG) THEN
435 ELSE IF (IDT.EQ.2) THEN
436 R = XC(1)**2 + XC(2)**2
437 IF (ISH.EQ.9) R = R + XC(3)**2
440 IF (ISH.EQ.5.OR.ISH.EQ.6.OR.ISH.EQ.9) THEN
447 PRINT *, ' GNEXT : Partially divided ',ISH,IAXIS
449 IF (NDIV.GT.1) IN2 = NDIV
451 ELSE IF (IDT.EQ.4) THEN
453 RXY = XC(1)**2 + XC(2)**2
455 IF (XC(3).NE.0.0) THEN
456 THET = RADDEG * ATAN (RXY/XC(3))
457 IF (THET.LT.0.0) THET = THET + 180.0
461 IF (THET.LE.ORIG) THEN
467 PRINT *, ' GNEXT : Partially divided ',ISH,IAXIS
469 IF (ISH.EQ.5.OR.ISH.EQ.7) THEN
471 IF (NDIV.GT.1) IN2 = NDIV
473 IF (XC(1).NE.0.0.OR.XC(2).NE.0.0) THEN
474 PHI = RADDEG * ATAN2 (XC(2), XC(1))
478 IF (ISH.EQ.6.OR.ISH.EQ.8) THEN
479 IF (PHI.LT.ORIG) THEN
486 IF (NDIV.GT.1) IN2 = NDIV
491 225 IF (IDT.EQ.1) THEN
495 X0(IAXIS) = ORIG + (IN - 0.5) * SDIV
496 IF (ISH.EQ.4.OR.(ISH.EQ.10.AND.IAXIS.NE.1)) THEN
497 CALL GCENT (IAXIS, X0)
500 XT(I) = XC(I) - X0(I)
503 ELSE IF (IDT.EQ.3) THEN
504 PH0 = DEGRAD * (ORIG + (IN - 0.5) * SDIV)
508 XT(I) = XC(I)*CPHR + XC(I+1)*SPHR
509 XT(I+1) = XC(I+1)*CPHR - XC(I)*SPHR
519 IF (IQ(JPARM-3).GT.1) THEN
532 CALL GNOBOX (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
533 ELSE IF (ISHT.EQ.2) THEN
534 CALL GNOTRA (XT, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
535 ELSE IF (ISHT.EQ.3) THEN
536 CALL GNOTRA (XT, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
538 CALL GNOTRP (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
540 ELSE IF (ISHT.LE.10) THEN
542 CALL GNOTUB (XT, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
543 ELSE IF (ISHT.EQ.6) THEN
544 CALL GNOTUB (XT, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
545 ELSE IF (ISHT.EQ.7) THEN
546 CALL GNOCON (XT, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
547 ELSE IF (ISHT.EQ.8) THEN
548 CALL GNOCON (XT, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
549 ELSE IF (ISHT.EQ.9) THEN
550 CALL GNOSPH (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
552 CALL GNOPAR (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
554 ELSE IF (ISHT.EQ.11) THEN
555 CALL GNOPGO (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
556 ELSE IF (ISHT.EQ.12) THEN
557 CALL GNOPCO (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
558 ELSE IF (ISHT.EQ.13) THEN
559 CALL GNOELT (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
560 ELSE IF (ISHT.EQ.14) THEN
561 CALL GNOHYP (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
562 ELSE IF (ISHT.EQ.28) THEN
563 CALL GSNGTR (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE,0)
564 ELSE IF (ISHT.EQ.NSCTUB) THEN
565 CALL GNOCTU (XT, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
567 PRINT *, ' GNEXT : No code for shape ', ISHT
571 IF (SAFE.LT.SAFETY) SAFETY = SAFE
572 IF (SNXT.LT.SNEXT) SNEXT = SNXT
581 * *** Calculate SNEXT and SAFETY with respect to the Mother
582 * *** SAFETY only for concave volumes if finite SNEXT
583 * *** has been found with respect to one of its contents
586 IF (SNEXT.LT.0.9*BIG) THEN
587 IF (.NOT.BTEST(IQ(JVO),2)) IACT = 0
589 JPAR = LQ(JGPAR-NLEVEL)
592 CALL GNBOX (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE )
593 ELSE IF (ISH.EQ.2) THEN
594 CALL GNTRAP (XC, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
595 ELSE IF (ISH.EQ.3) THEN
596 CALL GNTRAP (XC, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
598 CALL GNTRP (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
600 ELSE IF (ISH.LE.10) THEN
602 CALL GNTUBE (XC, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
603 ELSE IF (ISH.EQ.6) THEN
604 CALL GNTUBE (XC, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
605 ELSE IF (ISH.EQ.7) THEN
606 CALL GNCONE (XC, Q(JPAR+1), IACT, 1, SNEXT, SNXT, SAFE)
607 ELSE IF (ISH.EQ.8) THEN
608 CALL GNCONE (XC, Q(JPAR+1), IACT, 2, SNEXT, SNXT, SAFE)
609 ELSE IF (ISH.EQ.9) THEN
610 CALL GNSPHR (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
612 CALL GNPARA (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
614 ELSE IF (ISH.EQ.12) THEN
615 CALL GNPCON (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
616 ELSE IF (ISH.EQ.11) THEN
617 CALL GNPGON (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
618 ELSE IF (ISH.EQ.13) THEN
619 CALL GNELTU (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
620 ELSE IF (ISH.EQ.14) THEN
621 CALL GNHYPE (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
622 ELSE IF (ISH.EQ.28) THEN
623 CALL GSNGTR (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE,1)
624 ELSE IF (ISH.EQ.NSCTUB) THEN
625 CALL GNCTUB (XC, Q(JPAR+1), IACT, SNEXT, SNXT, SAFE)
627 PRINT *, ' GNEXT : No code for shape ', ISH
631 IF (SAFE.LT.SAFETY) SAFETY = SAFE
632 IF (SNXT.LT.SNEXT) SNEXT = SNXT
634 400 IF (GONLY(NLEVEL).EQ.0.) THEN
636 * *** Case of a 'NOT ONLY' volume -> step search
641 IF (ST.LE.0) GO TO 999
643 IF (ST.LE.EPSI3) THEN
654 XT(I) = X(I) +SN*X(I+3)
657 CALL GINVOL (XT, ISAME)
659 IF (ST.LT.EPSI2) GO TO 490
667 IF (ST.LT.EPSI2) THEN
676 IF (NN.GT.0) GO TO 420
679 490 IF (SN.LT.SNEXT) SNEXT = SN