5 * Revision 1.1.1.1 1995/10/24 10:20:51 cernlib
9 #include "geant321/pilot.h"
10 *CMZ : 3.21/02 29/03/94 15.41.29 by S.Giani
12 SUBROUTINE GMEDIV (JVO, IN, XC, IFL)
14 C. ******************************************************************
16 C. * Updates the common /GCVOLU/ and the structure JGPAR *
17 C. * for contents defined by division. *
19 C. * For IFL nonzero, it also checks if the point XC is inside *
20 C. * the content. It returns IN = 0, if the point is outside. *
21 C. * Otherwise, it transforms XC in the local system. *
23 C. * For IFL zero, IN is returned 0, if IN > NDIV. *
25 C. * Input : JVO, IN, XC, IFL *
26 C. * Output : IN, XC *
28 C. * Called by : GDRAW, GMEDIA *
29 C. * Authors : S.Banerjee, R.Brun, F.Bruyant, A.McPherson *
31 C. ******************************************************************
33 #include "geant321/gcbank.inc"
34 #include "geant321/gconsp.inc"
35 #include "geant321/gcpoly.inc"
36 #include "geant321/gcvolu.inc"
37 #if !defined(CERNLIB_SINGLE)
38 DOUBLE PRECISION DPHIO, TPIDEG, ONE
43 PARAMETER (TPIDEG=360,ONE=1)
46 DATA IDTYP / 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 2, 3, 1,
47 + 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 4, 3, 1, 1, 1,
50 C. ------------------------------------------------------------------
56 JVOT = LQ(JVOLUM-IVOT)
57 IF (NLEVEL.LT.NLDEV(NLEVEL)) THEN
60 * (case with structure JVOLUM locally developed)
61 JPAR = LQ(LQ(JVOLUM-LVOLUM(NLDEV(NLEVEL))))
62 IF (NLEVEL.EQ.NLDEV(NLEVEL)) GO TO 20
63 DO 10 ILEV = NLDEV(NLEVEL), NLEVEL-1
64 IF (IQ(JPAR+1).EQ.0) THEN
65 JPAR = LQ(JPAR-LINDEX(ILEV+1))
66 IF (JPAR.EQ.0) GO TO 20
67 ELSE IF (IQ(JPAR-3).GT.1) THEN
68 JPAR = LQ(JPAR-LINDEX(ILEV+1))
72 IF (ILEV.EQ.NLEVEL-1) THEN
85 30 IDT = IDTYP(IAXIS,ISH)
89 * Division along X, Y or Z axis
94 XTT = XTT - Q(LQ(JGPAR-NLEVEL)+IAXIS+4) * XC(3)
96 YT = XC(2) - Q(LQ(JGPAR-NLEVEL)+6) * XC(3)
97 XTT = XTT - Q(LQ(JGPAR-NLEVEL)+4) * YT
101 IN = (XTT -ORIG)/SDIV +1
102 ELSE IF (IDT.EQ.2) THEN
104 * Division along R axis
106 R = XC(1)**2 + XC(2)**2
107 IF (ISH.EQ.9) R = R + XC(3)**2
109 IF (ISH.EQ.5.OR.ISH.EQ.6.OR.ISH.EQ.9) THEN
110 IN = (R - ORIG) / SDIV + 1
111 ELSE IF (ISH.EQ.7.OR.ISH.EQ.8) THEN
112 IPAR = LQ(JGPAR-NLEVEL)
113 DR = 0.5 * (Q(IPAR+4) - Q(IPAR+2)) / Q(IPAR+1)
114 RMN = 0.5 * (Q(IPAR+4) + Q(IPAR+2)) + DR * XC(3)
115 DR = 0.5 * (Q(IPAR+5) - Q(IPAR+3)) / Q(IPAR+1)
116 RMX = 0.5 * (Q(IPAR+5) + Q(IPAR+3)) + DR * XC(3)
117 STP = (RMX - RMN) / NDIV
118 IN = (R - RMN) / STP + 1
120 IPAR = LQ(JGPAR-NLEVEL)
127 IPT = IPT + 3 * IZSEC
131 IF ((XC(3)-Q(IPT+3*IZ))*(XC(3)-Q(IPT+3*IZ+3))
134 IPT = IPT + 3 * IZSEC
141 50 POR1 = (Q(IPT+3) - XC(3)) / (Q(IPT+3) - Q(IPT))
142 POR2 = (XC(3) - Q(IPT)) / (Q(IPT+3) - Q(IPT))
143 RMN = Q(IPT+1) * POR1 + Q(IPT+4) * POR2
144 RMX = Q(IPT+2) * POR1 + Q(IPT+5) * POR2
147 DPH = Q(IPAR+2) / NPDV
149 IF (XC(1).NE.0..OR.XC(2).NE.0.) THEN
150 PHI = RADDEG * ATAN2 (XC(2), XC(1))
156 PH0 = MOD( ABS(PH0), 360.0 )
157 IF(SG.LE.0.0) PH0 = 360.0-PH0
160 PH = DEGRAD * (Q(IPAR+1) + (IPSEC - 0.5) * DPH)
161 R = XC(1) * COS(PH) + XC(2) * SIN(PH)
163 STP = (RMX - RMN) / NDIV
164 IN = (R - RMN) / STP + 1
166 ELSE IF (IDT.EQ.3) THEN
168 * Division along Phi axis
170 IF (XC(1).NE.0..OR.XC(2).NE.0.) THEN
171 PHI = RADDEG * ATAN2 (XC(2), XC(1))
177 DPHIO = MOD( ABS(DPHIO), TPIDEG)
178 IF(SG.LE.0.0) DPHIO=TPIDEG-DPHIO
180 ELSE IF (IDT.EQ.4) THEN
182 * Division along Theta axis
184 IF (XC(3).NE.0.0) THEN
185 RXY = SQRT (XC(1)**2 + XC(2)**2)
186 THET = RADDEG * ATAN (RXY/XC(3))
187 IF (THET.LT.0.0) THET = THET + 180.0
191 IN = (THET - ORIG) / SDIV + 1
195 60 IF (IN.GT.NDIV) IN = 0
196 IF (IN.LE.0) GO TO 999
199 IF (IQ(JPAR-3).GT.1) THEN
211 * Volume found at deeper level
215 LVOLUM(NLEVEL) = IVOT
216 NAMES(NLEVEL) = IQ(JVOLUM+IVOT)
220 GONLY(NLEVEL) = GONLY(NL1)
221 IF (LQ(LQ(JVOLUM-IVOT)).EQ.0) THEN
222 NLDEV(NLEVEL) = NLDEV(NL1)
224 NLDEV(NLEVEL) = NLEVEL
231 X0(IAXIS) = ORIG + (IN - 0.5) * SDIV
232 IF (ISH.EQ.4.OR.(ISH.EQ.10.AND.IAXIS.NE.1)) THEN
233 CALL GCENT (IAXIS, X0)
236 XC(1) = XC(1) - X0(1)
237 XC(2) = XC(2) - X0(2)
238 XC(3) = XC(3) - X0(3)
240 C***** Code Expanded From Routine: GTRMUL
242 C. ------------------------------------------------------------------
244 IF (GRMAT(10,NL1) .EQ. 0.0) THEN
245 GTRAN(1,NLEVEL) = GTRAN(1,NL1) + X0(1)
246 GTRAN(2,NLEVEL) = GTRAN(2,NL1) + X0(2)
247 GTRAN(3,NLEVEL) = GTRAN(3,NL1) + X0(3)
249 GRMAT(I,NLEVEL) = GRMAT(I,NL1)
250 GRMAT(I+1,NLEVEL) = GRMAT(I+1,NL1)
254 DXTEM1 = X0(1)*GRMAT(1,NL1) + X0(2)*GRMAT(4,NL1) + X0(3)*
256 DXTEM2 = X0(1)*GRMAT(2,NL1) + X0(2)*GRMAT(5,NL1) + X0(3)*
258 DXTEM3 = X0(1)*GRMAT(3,NL1) + X0(2)*GRMAT(6,NL1) + X0(3)*
261 GRMAT(I,NLEVEL) = GRMAT(I,NL1)
262 GRMAT(I+1,NLEVEL) = GRMAT(I+1,NL1)
264 GTRAN(1,NLEVEL) = GTRAN(1,NL1) + DXTEM1
265 GTRAN(2,NLEVEL) = GTRAN(2,NL1) + DXTEM2
266 GTRAN(3,NLEVEL) = GTRAN(3,NL1) + DXTEM3
268 C***** End of Code Expanded From Routine: GTRMUL
270 ELSE IF (IDT.EQ.3.OR.IDT.EQ.4) THEN
272 PH0 = DEGRAD * (ORIG + (IN - 0.5) * SDIV)
280 GTRAN(1,NLEVEL) = GTRAN(1,NL1)
281 GRMAT(1,NLEVEL) = GRMAT(1,NL1)*CPHR + GRMAT(4,NL1)*SPHR
282 GRMAT(4,NLEVEL) = GRMAT(4,NL1)*CPHR - GRMAT(1,NL1)*SPHR
283 GRMAT(7,NLEVEL) = GRMAT(7,NL1)
284 GTRAN(2,NLEVEL) = GTRAN(2,NL1)
285 GRMAT(2,NLEVEL) = GRMAT(2,NL1)*CPHR + GRMAT(5,NL1)*SPHR
286 GRMAT(5,NLEVEL) = GRMAT(5,NL1)*CPHR - GRMAT(2,NL1)*SPHR
287 GRMAT(8,NLEVEL) = GRMAT(8,NL1)
288 GTRAN(3,NLEVEL) = GTRAN(3,NL1)
289 GRMAT(3,NLEVEL) = GRMAT(3,NL1)*CPHR + GRMAT(6,NL1)*SPHR
290 GRMAT(6,NLEVEL) = GRMAT(6,NL1)*CPHR - GRMAT(3,NL1)*SPHR
291 GRMAT(9,NLEVEL) = GRMAT(9,NL1)
293 XTT = XC(1) * CPHR + XC(2) * SPHR
294 XC(2) = XC(2) * CPHR - XC(1) * SPHR
297 IF (PH0.EQ.0.0.AND.GRMAT(10,NL1).EQ.0.0) THEN
298 GRMAT(10,NLEVEL) = 0.0
300 GRMAT(10,NLEVEL) = 1.0
302 IF (ISH.EQ.11) IPSEC = 1
305 GTRAN(1,NLEVEL) = GTRAN(1,NL1)
306 GTRAN(2,NLEVEL) = GTRAN(2,NL1)
307 GTRAN(3,NLEVEL) = GTRAN(3,NL1)
309 GRMAT(I,NLEVEL) = GRMAT(I,NL1)
310 GRMAT(I+1,NLEVEL) = GRMAT(I+1,NL1)
314 IQ(JGPAR+NLEVEL) = NPAR
315 LQ(JGPAR-NLEVEL) = JPAR