]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/ggeom/gvdphi.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Allow any Cherenkov-like particle to be transported
[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gvdphi.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1995/10/24 10:20:57 cernlib
6* Geant
7*
8*
9#include "geant321/pilot.h"
10*CMZ : 3.21/04 02/03/95 11.58.00 by S.Giani
11*-- Author :
12 SUBROUTINE GVDPHI(ISH,IROT,DX,PARS,CL,CH,IERR)
13C.
14C. **********************************************************
15C. * *
16C. * ROUTINE TO FIND THE PHI LIMITS OF THE OBJECT SHAPE *
17C. * ISH ROTATED BY MATRIX IROT AND DISPLACED BY VECTOR *
18C. * DX. IT HAS NPAR PARAMTERS IN THE ARRAY PARS. THE *
19C. * LOWER LIMIT IS RETURNED IN CL AND THE HIGHER IN CH. *
20C. * NOTE THE OBJECT IS CONTAINED IN THE RANGE OF *
21C. * INCREASING PHI FROM CL TO CH THOUGH CL AND CH ARE *
22C. * FORCED TO LIE IN THE RANGE 0.0 TO 360.0 SO THAT THE *
23C. * VALUE OF CL CAN BE HIGHER THAN THAT OF CH. *
24C. * *
25C. * ==>Called by : GVDLIM modified *
26C. * Author S.Giani ********* *
27C. * *
28C. **********************************************************
29C.
30#include "geant321/gcbank.inc"
31#include "geant321/gconsp.inc"
32#include "geant321/gcshno.inc"
33 DIMENSION DX(3),PARS(50),X(3),XT(3)
34C.
35C. -------------------------------------------
36C.
37 IERR=1
38C
39 DXS=DX(1)*DX(1)+DX(2)*DX(2)
40 IF(DXS.GT.0.0) DXS=SQRT(DXS)
41C
42 IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 40
43C
44C
45C CUBOIDS, TRAPEZOIDS, PARALLELEPIPEDS.
46C
47 IERR=0
48 CL=0.0
49 CH=360.0
50C
51C IF IN DOUBT SET IT TO FULL RANGE.
52C
53 IF(DXS.LE.0.0) GO TO 999
54C
55 PHC=90.
56 IF(DX(1).NE.0.)PHC=ATAN2(DX(2),DX(1))*RADDEG
57 IF(DX(1).EQ.0..AND.DX(2).LT.0.)PHC=-90.
58 IF(PHC.LT.0.0) PHC=PHC+360
59 PL=0.0
60 PH=0.0
61C
62 DO 30 IP=1,8
63C
64C THIS IS A LOOP OVER THE 8 CORNERS.
65C FIRST FIND THE LOCAL COORDINATES.
66C
67 IF(ISH.EQ.28) THEN
68C
69C General twisted trapezoid.
70C
71 IL=(IP+1)/2
72 I0=IL*4+11
73 IS=(IP-IL*2)*2+1
74 X(3)=PARS(1)*IS
75 X(1)=PARS(I0)+PARS(I0+2)*X(3)
76 X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
77 GO TO 20
78C
79 ENDIF
80C
81 IP3=ISH+2
82 IF(ISH.EQ.10) IP3=3
83 IF(ISH.EQ.4) IP3=1
84 X(3)=PARS(IP3)
85 IF(IP.LE.4) X(3)=-X(3)
86 IP2=3
87 IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
88 IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
89 IF(ISH.EQ.4) IP2=4
90 IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
91 X(2)=PARS(IP2)
92 IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
93 IP1=1
94 IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
95 IF(ISH.EQ.4) IP1=5
96 IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
97 IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
98 X(1)=PARS(IP1)
99 IF(MOD(IP,2).EQ.1) X(1)=-X(1)
100C
101 IF(ISH.NE.10) GO TO 10
102 X(1)=X(1)+X(2)*PARS(4)+X(3)*PARS(5)
103 X(2)=X(2)+X(3)*PARS(6)
104 10 CONTINUE
105C
106 IF(ISH.NE.4) GO TO 20
107 IP4=7
108 IF(X(3).GT.0.0) IP4=11
109 X(1)=X(1)+X(2)*PARS(IP4)+X(3)*PARS(2)
110 X(2)=X(2)+X(3)*PARS(3)
111 20 CONTINUE
112C
113C ROTATE.
114C
115 JROT=LQ(JROTM-IROT)
116 XT(1)=X(1)
117 XT(2)=X(2)
118 XT(3)=X(3)
119 IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
120C
121 XPT=DXS+(DX(1)*XT(1)+DX(2)*XT(2))/DXS
122 YPT=(DX(1)*XT(2)-DX(2)*XT(1))/DXS
123C
124 IF(YPT.EQ.0.0.AND.XPT.EQ.0.0) GO TO 999
125 P=ATAN2(YPT,XPT)
126 IF(P.GT.PI) P=P-PI*2.0
127 IF(P.LT.PL) PL=P
128 IF(P.GT.PH) PH=P
129C
130C
131 30 CONTINUE
132C
133C
134 IF(PH-PL.GT.PI) GO TO 999
135 CL=PHC+PL*RADDEG
136 CH=PHC+PH*RADDEG
137C
138*** SG = SIGN(1.0,CL)
139*** CL = MOD( ABS(CL),360.0 )
140*** IF(SG.LE.0.0) CL=360.-CL
141*** SG=SIGN(1.0,CH)
142*** CH = MOD( ABS(CH),360.0 )
143*** IF(SG.LE.0.0) CH=360.-CH
144C
145 GO TO 999
146C
147 40 CONTINUE
148 MYFLAG=0
149 IF(ISH.EQ.11)THEN
150 NZLAST=PARS(4)
151 IZLAST=2+3*NZLAST
152 CLZ=PARS(5)
153 CHZ=PARS(IZLAST)
154 DZ2=ABS(CHZ-CLZ)
155 DZ=DZ2*.5
156 TMPRAD=0.
157 TMPMIN=100000.
158 DO 145 I=7,IZLAST+2,3
159 IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
160 IF(PARS(I-1).LT.TMPMIN)TMPMIN=PARS(I-1)
161 145 CONTINUE
162 RLOW=TMPMIN
163 PHIMIN=PARS(1)
164 PHIMAX=PHIMIN+PARS(2)
165 AANG=ABS(PHIMAX-PHIMIN)
166 NANG=PARS(3)
167 AATMAX=NANG*360./AANG
168 LATMAX=AATMAX
169 ALA=AATMAX-LATMAX
170 IF(ALA.GT..5)LATMAX=LATMAX+1
171 AFINV=1./COS(PI/LATMAX)
172 FINV=ABS(AFINV)
173 RM=TMPRAD*FINV
174 IF(PARS(2).EQ.360)THEN
175 MYFLAG=5
176 ELSE
177 PHIMIN=PARS(1)*DEGRAD
178 PHIMAX=(PARS(1)+PARS(2))*DEGRAD
179 MYFLAG=6
180 ENDIF
181 ELSEIF(ISH.EQ.12)THEN
182 NZLAST=PARS(3)
183 IZLAST=1+3*NZLAST
184 CLZ=PARS(4)
185 CHZ=PARS(IZLAST)
186 DZ2=ABS(CHZ-CLZ)
187 DZ=DZ2*.5
188 TMPRAD=0.
189 TMPMIN=100000.
190 DO 146 I=6,IZLAST+2,3
191 IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
192 IF(PARS(I-1).LT.TMPMIN)TMPMIN=PARS(I-1)
193 146 CONTINUE
194 RM=TMPRAD
195 RLOW=TMPMIN
196 IF(PARS(2).EQ.360)THEN
197 MYFLAG=5
198 ELSE
199 PHIMIN=PARS(1)*DEGRAD
200 PHIMAX=(PARS(1)+PARS(2))*DEGRAD
201 MYFLAG=6
202 ENDIF
203 ELSEIF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14)THEN
204 GO TO 80
205 ENDIF
206C
207C TUBES AND CONES.
208C
209 IERR=0
210 CL=0.0
211 CH=360.0
212C
213C WHEN IN DOUBT SET TO FULL RANGE.
214C
215 IF(MYFLAG.EQ.0)RM=PARS(2)
216 IF(ISH.LE.6.OR.ISH.EQ.NSCTUB.OR.MYFLAG.NE.0) GO TO 50
217**
218 IF(ISH.EQ.13) THEN
219**
220** approxime to a cylinder whit radius
221** equal to the ellipse major axis
222**
223 IF(PARS(1).GT.RM) RM=PARS(1)
224 GOTO 50
225 ENDIF
226 IF(ISH.EQ.14) THEN
227 RM = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
228 GO TO 50
229 ENDIF
230 RM=PARS(3)
231 IF(PARS(5).GT.RM) RM=PARS(5)
232C
233 50 CONTINUE
234C
235 IF(DXS.GT.RM) GO TO 70
236 IF(ISH.EQ.5.OR.ISH.EQ.7.OR.ISH.EQ.14.OR.MYFLAG.EQ.5) GO TO 999
237 IF(ISH.EQ.13) GOTO 999
238* Here we treat the CONS and TUBS
239* This is the simple case, no rotation
240* Compute the position of the limits on
241* the X-Y plane.
242 IF(MYFLAG.EQ.0.AND.ISH.LE.6)THEN
243 PHIMIN=PARS(4)*DEGRAD
244 PHIMAX=PARS(5)*DEGRAD
245 RLOW=PARS(1)
246 ELSEIF(MYFLAG.EQ.0.AND.ISH.EQ.8)THEN
247 PHIMIN=PARS(6)*DEGRAD
248 PHIMAX=PARS(7)*DEGRAD
249 RLOW=MIN(PARS(2),PARS(4))
250 ENDIF
251 CL=PHIMIN*RADDEG
252 CH=PHIMAX*RADDEG
253 IF(DXS.NE.0.)THEN
254 DDX1 = DX(1)+RM*COS(PHIMIN)
255 DDY1 = DX(2)+RM*SIN(PHIMIN)
256 DDX2 = DX(1)+RM*COS(PHIMAX)
257 DDY2 = DX(2)+RM*SIN(PHIMAX)
258 CLA = ATAN2(DDY1,DDX1)*RADDEG
259 CHA = ATAN2(DDY2,DDX2)*RADDEG
260 DDX1 = DX(1)+RLOW*COS(PHIMIN)
261 DDY1 = DX(2)+RLOW*SIN(PHIMIN)
262 DDX2 = DX(1)+RLOW*COS(PHIMAX)
263 DDY2 = DX(2)+RLOW*SIN(PHIMAX)
264 CLB = ATAN2(DDY1,DDX1)*RADDEG
265 CHB = ATAN2(DDY2,DDX2)*RADDEG
266 CL=MIN(CLA,CLB,CHA,CHB)
267 CH=MAX(CLA,CLB,CHA,CHB)
268 IF((CH-CL).GT.(PHIMAX-PHIMIN)*RADDEG)THEN
269 IF(ISH.EQ.6.OR.ISH.EQ.8.OR.MYFLAG.EQ.6)THEN
270 IF(DXS.GT.RLOW)THEN
271 CL=0.
272 CH=360.
273 ENDIF
274 ELSE
275 CL=0.
276 CH=360.
277 ENDIF
278 ENDIF
279 ENDIF
280C
281 60 CONTINUE
282C
283 IF(IROT.EQ.0) GO TO 65
284 IF(CL.EQ.0..AND.CH.EQ.360.)GOTO 65
285 JROT=LQ(JROTM-IROT)
286 IF(Q(JROT+15).NE.0.0.AND.Q(JROT+15).NE.180.0)THEN
287 CL=0.
288 CH=360.
289 GO TO 999
290 ENDIF
291C
292 PHX=Q(JROT+12)
293 PHY=Q(JROT+14)
294 IF(PHY.LT.PHX) PHY=PHY+360.0
295 ISPH=1
296 IF(PHY-PHX.GT.180.0) ISPH=-1
297 IF(DXS.NE.0.)THEN
298 PHI1 = ISPH*PHIMIN+PHX*DEGRAD
299 PHI2 = ISPH*PHIMAX+PHX*DEGRAD
300 DDX1 = DX(1)+RM*COS(PHI1)
301 DDY1 = DX(2)+RM*SIN(PHI1)
302 DDX2 = DX(1)+RM*COS(PHI2)
303 DDY2 = DX(2)+RM*SIN(PHI2)
304 CLA = ATAN2(DDY1,DDX1)*RADDEG
305 CHA = ATAN2(DDY2,DDX2)*RADDEG
306 DDX1 = DX(1)+RLOW*COS(PHI1)
307 DDY1 = DX(2)+RLOW*SIN(PHI1)
308 DDX2 = DX(1)+RLOW*COS(PHI2)
309 DDY2 = DX(2)+RLOW*SIN(PHI2)
310 CLB = ATAN2(DDY1,DDX1)*RADDEG
311 CHB = ATAN2(DDY2,DDX2)*RADDEG
312 CL=MIN(CLA,CLB,CHA,CHB)
313 CH=MAX(CLA,CLB,CHA,CHB)
314 ELSE
315 CL=ISPH*CL+PHX
316 CH=ISPH*CH+PHX
317 ENDIF
318 IF(ISPH.EQ.1) GO TO 65
319 CHT=CH
320 CH=CL
321 CL=CHT
322C
323 65 CONTINUE
324C
325*** SG=SIGN(1.0,CL)
326*** CL = MOD( ABS(CL),360.0 )
327*** IF(SG.LE.0.0) CL=360.-CL
328*** SG=SIGN(1.0,CH)
329*** CH = MOD( ABS(CH),360.0 )
330*** IF(SG.LE.0.0) CH=360.-CH
331C
332 GO TO 999
333C
334 70 CONTINUE
335C
336C DISPLACEMENT GREATER THAN MAXIMUM RADIUS SO
337C ASSUME COMPLETE TUBE TO GENERATE 'WORST CASE'.
338C
339 IF(MYFLAG.EQ.0)DZ=PARS(3)
340 IF(ISH.EQ.NSCTUB) THEN
341 S1 = (1.0-PARS(8))*(1.0+PARS(8))
342 IF( S1 .GT. 0.0) S1 = SQRT(S1)
343 S2 = (1.0-PARS(11))*(1.0+PARS(11))
344 IF( S2 .GT. 0.0) S2 = SQRT(S2)
345 IF( S2 .GT. S1 ) S1 = S2
346 DZ = DZ+RM*S1
347 ELSEIF(ISH.GT.6.AND.ISH.NE.13.AND.ISH.NE.14.AND.MYFLAG.EQ.0)THEN
348 DZ=PARS(1)
349 ENDIF
350C
351 X(1)=0.0
352 X(2)=0.0
353 X(3)=1.0
354C
355C LOCAL Z AXIS.
356C
357 JROT=LQ(JROTM-IROT)
358 XT(1)=X(1)
359 XT(2)=X(2)
360 XT(3)=X(3)
361 IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
362C
363 COST=ABS(DX(1)*XT(1)+DX(2)*XT(2))
364 COS2=ABS(DX(1)*XT(2)-DX(2)*XT(1))
365 SINT=(DXS+COST)*(DXS-COST)
366 SIN2=(DXS+COS2)*(DXS-COS2)
367 IF(SINT.GT.0.0) SINT=SQRT(SINT)
368 IF(SIN2.GT.0.0) SIN2=SQRT(SIN2)
369C
370 XPT=DXS-(COST*DZ+SINT*RM)/DXS
371C
372 IF(XPT.LE.0.0) GO TO 999
373 YPT=(SIN2*RM+COS2*DZ)/DXS
374 DP=ATAN(YPT/XPT)
375C
376 P0=ATAN2(DX(2),DX(1))
377 CL=(P0-DP)*RADDEG
378 CH=(P0+DP)*RADDEG
379C
380*** SG=SIGN(1.0,CL)
381*** CL = MOD( ABS(CL),360.0 )
382*** IF(SG.LE.0.0) CL=360.-CL
383*** SG=SIGN(1.0,CH)
384*** CH = MOD( ABS(CH),360.0 )
385*** IF(SG.LE.0.0) CH=360.-CH
386C
387 GO TO 999
388C
389 80 CONTINUE
390 IF(ISH.GT.9.AND.MYFLAG.EQ.0) GO TO 999
391C
392C SPHERE.
393C
394 IERR=0
395 CL=0.0
396 CH=360.0
397C
398 IF(IROT.NE.0.OR.DXS.GT.0.0) GO TO 90
399C
400C UNROTATED AND CENTERED.
401C
402 CL=PARS(5)
403 CH=PARS(6)
404C
405*** SG=SIGN(1.0,CL)
406*** CL = MOD( ABS(CL),360.0 )
407*** IF(SG.LE.0.0) CL=360.-CL
408*** SG=SIGN(1.0,CH)
409*** CH = MOD( ABS(CH),360.0 )
410*** IF(SG.LE.0.0) CH=360.-CH
411C
412 GO TO 999
413C
414 90 CONTINUE
415C
416C ROTATED OR NOT CENTERED.
417C
418 IF(DXS.LT.PARS(2)) GO TO 999
419 P0=ATAN2(DX(2),DX(1))
420 DP=ASIN(PARS(2)/DXS)
421 CL=(P0-DP)*RADDEG
422 CH=(P0+DP)*RADDEG
423C
424*** SG=SIGN(1.0,CL)
425*** CL = MOD( ABS(CL),360.0 )
426*** IF(SG.LE.0.0) CL=360.-CL
427*** SG=SIGN(1.0,CH)
428*** CH = MOD( ABS(CH),360.0 )
429*** IF(SG.LE.0.0) CH=360.-CH
430C
431 999 CONTINUE
432 END