--- /dev/null
+*
+* $Id$
+*
+* $Log$
+* Revision 1.3 1998/02/09 16:48:33 japost
+* Simone's fix for MANY volumes in gmedia.
+*
+* Revision 1.2 1996/09/30 14:25:07 ravndal
+* Windows NT related modifications
+*
+* Revision 1.1.1.1 1995/10/24 10:20:51 cernlib
+* Geant
+*
+*
+#include "geant321/pilot.h"
+#if !defined(CERNLIB_OLD)
+*CMZ : 3.21/02 29/03/94 15.24.17 by S.Giani
+*-- Author :
+ SUBROUTINE GMEDIA (X, NUMED)
+C.
+C. ******************************************************************
+C. * *
+C. * Finds in which volume/medium the point X is, and updates the *
+C. * common /GCVOLU/ and the structure JGPAR accordingly. *
+C. * *
+C. * NUMED returns the tracking medium number, or 0 if point is *
+C. * outside the experimental setup. *
+C. * *
+C. * Called by : GTREVE, GLTRAC, 'User' *
+C. * Authors : R.Brun, F.Bruyant, A.McPherson *
+C. * S.Giani. *
+C. * *
+C. * Modified by S.Giani (1993) to perform the search according *
+C. * to the new 'virtual divisions' algorithm and to build the *
+C. * stack of the 'possible overlapping volumes' in the case of *
+C. * MANY volumes. Any kind of boolean operation is now possible.*
+C. * Divisions along arbitrary axis are now possible. *
+C. * *
+C. ******************************************************************
+C.
+#include "geant321/gcflag.inc"
+#include "geant321/gckine.inc"
+#include "geant321/gcbank.inc"
+#include "geant321/gcvolu.inc"
+#include "geant321/gctrak.inc"
+#if defined(CERNLIB_USRJMP)
+#include "geant321/gcjump.inc"
+#endif
+#include "geant321/gcvdma.inc"
+#include "geant321/gchvir.inc"
+C.
+ DIMENSION X(*)
+ REAL XC(6)
+ LOGICAL BTEST
+ CHARACTER*4 NAME
+C.
+C. ------------------------------------------------------------------
+*
+ nvmany=0
+ nfmany=0
+ new2fl=0
+*
+ IF (NLEVEL.EQ.0) CALL GMEDIN
+*
+* SECTION I: The /GCVOLU/ table contains the initial guess for a path
+* in the geometry tree on which X may be found. Look along this
+* path until X is found inside. This is the starting position.
+* If this is an ONLY volume with no daughters, we are done;
+* otherwise reset search record variables, proceed to section II.
+*
+* The information contained in INFROM has to be invalidated
+* because it has no meaning for the subsequent tracking. INFR
+* is a local variable used to optimise the search in the
+* geometry tree.
+*
+ INFROM = 0
+*
+* *** Check if point is in current volume
+*
+ INFR = 0
+ JVIN = 0
+C***** Code Expanded From Routine: GTRNSF
+C
+ 100 IF (GRMAT(10,NLEVEL) .EQ. 0.) THEN
+ XC(1) = X(1) - GTRAN(1,NLEVEL)
+ XC(2) = X(2) - GTRAN(2,NLEVEL)
+ XC(3) = X(3) - GTRAN(3,NLEVEL)
+*
+ ELSE
+ XL1 = X(1) - GTRAN(1,NLEVEL)
+ XL2 = X(2) - GTRAN(2,NLEVEL)
+ XL3 = X(3) - GTRAN(3,NLEVEL)
+ XC(1) = XL1*GRMAT(1,NLEVEL) + XL2*GRMAT(2,NLEVEL) + XL3*GRMAT(3
+ + ,NLEVEL)
+ XC(2) = XL1*GRMAT(4,NLEVEL) + XL2*GRMAT(5,NLEVEL) + XL3*GRMAT(6
+ + ,NLEVEL)
+ XC(3) = XL1*GRMAT(7,NLEVEL) + XL2*GRMAT(8,NLEVEL) + XL3*GRMAT(9
+ + ,NLEVEL)
+
+ ENDIF
+ xc(4)=0.
+ xc(5)=0.
+ xc(6)=0.
+C***** End of Code Expanded From Routine: GTRNSF
+*
+ JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
+ JPAR = LQ(JGPAR-NLEVEL)
+ CALL GINME (XC, Q(JVO+2), Q(JPAR+1), IYES)
+ IF (IYES.EQ.0) THEN
+*
+* ** Point not in current volume, go up the tree
+*
+ IF (NLEVEL.GT.1) THEN
+ NLEVEL = NLEVEL -1
+ JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
+ NIN = Q(JVO+3)
+ IF(NIN.GT.0) THEN
+*
+* Do not set INFR whne going up the tree. GMEDIA can be called
+* by the user and it should not assume that the previous
+* position has something to do with the current search. INFR
+* is otherwise useful when searching in a 'MANY' volume
+* configuration. This statement is commented for the above reason.
+*
+* INFR =LINDEX(NLEVEL+1)
+ ELSE
+ INFR =0
+ ENDIF
+ GO TO 100
+ ELSE
+*
+* * Point is outside setup
+*
+ NUMED = 0
+ GO TO 999
+ ENDIF
+ ENDIF
+*
+* ** Point is in current volume
+*
+ IF(INFR .GT.0) THEN
+ JIN=LQ(JVO-INFR )
+ IQ(JIN) = IBSET(IQ(JIN),4)
+ JVIN = JIN
+ ENDIF
+* To avoid starting from the protuding part of a MANY volume
+ IF(GONLY(NLEVEL).EQ.0.) THEN
+ NLEVEL = NLEVEL -1
+ GO TO 100
+ ENDIF
+ NLMIN = NLEVEL
+ NLMANY = 0
+*
+* SECTION II: X is found inside current node at NLEVEL in /GCVOLU/.
+* Search all contents recursively for any containing X.
+* Take the first one found, if any, and continue at that
+* level, incrementing NLEVEL and extending /GCVOLU/ tables.
+* This is continued until a level is reached where X is not
+* found in any of the contents, or there are no contents.
+* Note: Since Section II is re-entered from Section III, a blocking word
+* is used to mark those contents already checked. Upon exit from Section
+* II, these blocking words are cleared at NLEVEL, but may remain set in
+* levels between NLEVEL-1 and NLMIN, if any. They must be cleared at exit.
+*
+* ** Check contents, if any
+*
+ 200 JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
+ NIN = Q(JVO+3)
+ if(raytra.eq.1..and.imyse.eq.1)then
+ CALL UHTOC(NAMES(NLEVEL),4,NAME,4)
+ CALL GFIND(NAME,'SEEN',ISSEEN)
+ if(isseen.eq.-2.or.isseen.eq.-1)goto 300
+ endif
+*
+* * Case with no contents
+*
+ IF (NIN.EQ.0) THEN
+ GO TO 300
+*
+* * Case with contents defined by division
+*
+ ELSEIF (NIN.LT.0) THEN
+ CALL GMEDIV (JVO, IN, XC, 1)
+ IF (IN.GT.0) THEN
+ INFR = 0
+ GO TO 200
+ ENDIF
+*
+* * Case with contents positioned
+*
+ ELSE
+ if(nin.gt.1)then
+ clmoth=q(jvirt+4*(LVOLUM(NLEVEL)-1)+3)
+ chmoth=q(jvirt+4*(LVOLUM(NLEVEL)-1)+4)
+ ndivto=q(jvirt+4*(LVOLUM(NLEVEL)-1)+2)
+ iaxis =q(jvirt+4*(LVOLUM(NLEVEL)-1)+1)
+ if(iaxis.le.3)then
+ ivdiv=((xc(iaxis)-clmoth)*ndivto/(chmoth-clmoth))+1
+ if(ivdiv.lt.1)then
+ ivdiv=1
+ elseif(ivdiv.gt.ndivto)then
+ ivdiv=ndivto
+ endif
+ else
+ call gfcoor(xc,iaxis,cx)
+ if(iaxis.eq.6)then
+ if((cx-clmoth).lt.-1.)then
+ cx=cx+360.
+ elseif((cx-chmoth).gt.1.)then
+ cx=cx-360.
+ endif
+ if(cx.gt.chmoth)then
+ cx=chmoth
+ elseif(cx.lt.clmoth)then
+ cx=clmoth
+ endif
+ endif
+ ivdiv=((cx-clmoth)*ndivto/(chmoth-clmoth))+1
+ if(ivdiv.lt.1)then
+ ivdiv=1
+ elseif(ivdiv.gt.ndivto)then
+ ivdiv=ndivto
+ endif
+ endif
+ jvdiv=lq(jvirt-LVOLUM(NLEVEL))
+ iofset=iq(jvdiv+ivdiv)
+ ncont=iq(jvdiv+iofset+1)
+ jcont=jvdiv+iofset+1
+ if(ncont.eq.0)goto 260
+ else
+ JCONT = LQ(JVO-NIN-1)+1
+ NCONT = 1
+ endif
+*
+* For each selected content in turn, check if point is inside
+*
+ DO 259 ICONT=1,NCONT
+ if(nin.eq.1)then
+ in=1
+ else
+ IN = IQ(JCONT+ICONT)
+ endif
+ IF(IN.EQ.0) THEN
+*
+* If the value IQ(JCONT+ICONT)=0 then we are back in the mother.
+* So jump to 260, the search is finished. Clean-up should be done
+* only up to ICONT-1, so we set:
+*
+ NCONT=ICONT-1
+ GOTO 260
+ ELSE
+ JIN = LQ(JVO-IN)
+ IF (.NOT.BTEST(IQ(JIN),4)) THEN
+ CALL GMEPOS (JVO, IN, XC, 1)
+ IF (IN.GT.0) THEN
+ new2fl=0
+ IF (GONLY(NLEVEL).NE.0.) THEN
+ NLMANY = 0
+ nvmany = 0
+ nfmany = 0
+ ENDIF
+ INFR = 0
+ GO TO 200
+ ELSE
+ IQ(JIN) = IBSET(IQ(JIN),4)
+ ENDIF
+ ENDIF
+ ENDIF
+ 259 CONTINUE
+*
+ 260 IF(NCONT.EQ.NIN) THEN
+ DO 268 IN=1,NIN
+ JIN = LQ(JVO-IN)
+ IQ(JIN) = IBCLR(IQ(JIN),4)
+ 268 CONTINUE
+ ELSE
+ DO 269 ICONT=1,NCONT
+ if(nin.eq.1)then
+ in=1
+ else
+ IN = IQ(JCONT+ICONT)
+ endif
+ JIN = LQ(JVO-IN)
+ IQ(JIN) = IBCLR(IQ(JIN),4)
+ 269 CONTINUE
+ IF(INFR .GT.0) THEN
+ JIN = LQ(JVO-INFR )
+ IQ(JIN) = IBCLR(IQ(JIN),4)
+ ENDIF
+ ENDIF
+*
+ ENDIF
+*
+* SECTION III: X is found at current node (NLEVEL in /GCVOLU) but not in
+* any of its contents, if any. If this is a MANY volume,
+* save it as a candidate best-choice, and continue the search
+* by backing up the tree one node and proceed to Section II.
+* If this is an ONLY volume, proceed to Section IV.
+*
+* *** Point is in current volume/medium, and not in any content
+*
+ 300 IF (GONLY(NLEVEL).EQ.0.) THEN
+*
+* ** Lowest level is 'NOT ONLY'
+*
+ IF (NLEVEL.GT.NLMANY) THEN
+ CALL GSCVOL
+ NLMANY = NLEVEL
+ nfmany=nvmany+1
+ ENDIF
+ if(new2fl.eq.0)then
+ nvmany=nvmany+1
+ manyle(nvmany)=nlevel
+ do 401 i = 1,nlevel
+ manyna(nvmany,i)=names(i)
+ manynu(nvmany,i)=number(i)
+ 401 continue
+ endif
+*
+* * Go up the tree up to a volume with positioned contents
+*
+ new2fl=-1
+ 310 INFR = LINDEX(NLEVEL)
+ NLEVEL = NLEVEL -1
+ JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
+ NIN = Q(JVO+3)
+ IF (NIN.LT.0) GO TO 310
+*
+C***** Code Expanded From Routine: GTRNSF
+C
+ IF (GRMAT(10,NLEVEL) .EQ. 0.) THEN
+ XC(1) = X(1) - GTRAN(1,NLEVEL)
+ XC(2) = X(2) - GTRAN(2,NLEVEL)
+ XC(3) = X(3) - GTRAN(3,NLEVEL)
+*
+ ELSE
+ XL1 = X(1) - GTRAN(1,NLEVEL)
+ XL2 = X(2) - GTRAN(2,NLEVEL)
+ XL3 = X(3) - GTRAN(3,NLEVEL)
+ XC(1) = XL1*GRMAT(1,NLEVEL) + XL2*GRMAT(2,NLEVEL) + XL3*
+ + GRMAT(3,NLEVEL)
+ XC(2) = XL1*GRMAT(4,NLEVEL) + XL2*GRMAT(5,NLEVEL) + XL3*
+ + GRMAT(6,NLEVEL)
+ XC(3) = XL1*GRMAT(7,NLEVEL) + XL2*GRMAT(8,NLEVEL) + XL3*
+ + GRMAT(9,NLEVEL)
+
+ ENDIF
+C***** End of Code Expanded From Routine: GTRNSF
+*
+ JIN = LQ(JVO-INFR )
+ IQ(JIN) = IBSET(IQ(JIN),4)
+ NLMIN = MIN(NLEVEL,NLMIN)
+ GO TO 200
+ ENDIF
+*
+* SECTION IV: This is the end of the search. The current node (NLEVEL
+* in /GCVOLU/) is the lowest ONLY volume in which X is found.
+* If X was also found in any of its contents, they are MANY
+* volumes: the best-choice is the one among them at the greatest
+* level in the tree, and it is stored. Otherwise the current
+* volume is the solution. Before exit, all of the blocking
+* words leftover in the tree must be reset to zero.
+* Note: A valid structure is assumed, in which no ONLY volumes overlap.
+* If this rule is violated, or if a daughter is not entirely contained
+* within the mother volume, the results are unpredictable.
+*
+ DO 419 NL=NLMIN,NLEVEL-1
+ JVO = LQ(JVOLUM-LVOLUM(NL))
+ NIN = Q(JVO+3)
+ DO 418 IN=1,NIN
+ JIN = LQ(JVO-IN)
+ IQ(JIN) = IBCLR(IQ(JIN),4)
+ 418 CONTINUE
+ 419 CONTINUE
+*
+ if(nlmany.eq.0)then
+ nvmany=0
+ nfmany=0
+ endif
+ IF (NLMANY.GT.0) CALL GFCVOL
+ JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
+ IF(JVIN.NE.0) IQ(JVIN) = IBCLR(IQ(JVIN),4)
+ NUMED = Q(JVO+4)
+* END GMEDIA
+ 999 IF(JGSTAT.NE.0) CALL GFSTAT(2)
+ END
+#endif