]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/ggeom/gvdthe.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gvdthe.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
d43b40e2 5* Revision 1.1.1.1 1999/05/18 15:55:17 fca
6* AliRoot sources
7*
fe4da5cc 8* Revision 1.1.1.1 1995/10/24 10:20:57 cernlib
9* Geant
10*
11*
12#include "geant321/pilot.h"
13*CMZ : 3.21/02 29/03/94 15.41.31 by S.Giani
14*-- Author :
15 SUBROUTINE GVDTHE(ISH,IROT,DX,PARS,CL,CH,IERR)
16C.
17C. ******************************************************************
18C. * *
19C. * ROUTINE TO COMPUTE THE THETA LIMITS FOR VOLUME OF SHAPE *
20C. * ISH ROTATED BY MATRIX IROT AND DISPLACED BY VECTOR DX. *
21C. * THE VOLUME HAS NPAR PARAMETERS IN THE ARRAY PARS. THE LOWER *
22C. * LIMIT IS RETURNED IN CL THE HIGHER IN CH. IF THE *
23C. * CALCULATION CANNOT BE MADE IERR IS SET TO 1 OTHERWISE IT *
24C. * IS SET TO 0. *
25C. * *
26C. * ==>Called by : GVDLIM *
27C. * Author S.Giani ********* *
28C. * *
29C. ******************************************************************
30C.
31#include "geant321/gcbank.inc"
32#include "geant321/gconsp.inc"
33#include "geant321/gcshno.inc"
d43b40e2 34 DIMENSION DX(3),PARS(100),X(3),XT(3),X1(3),X2(3),XT1(3),XT2(3)
fe4da5cc 35C.
36C. ----------------------------------------------
37C.
38 IERR=1
39C
40 DXS=DX(1)*DX(1)+DX(2)*DX(2)
41 IF(DXS.GT.0.0) DXS=SQRT(DXS)
42C
43 IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 130
44C
45C CUBES, TRAPEZOIDS AND PARALLELEPIPEDS.
46C
47 IERR=0
48 CL=0.0
49 CH=180.0
50C
51 IF(DXS.LE.0.0) GO TO 999
52C
53 TH=90.
54 IF(ABS(DX(3)).LT.1.0E-06)GO TO 5
55 TH=ATAN(DXS/DX(3))*RADDEG
56 IF(TH.LT.0.0) TH=TH+180.0
57 5 TL=TH
58C
59 DO 50 IP=1,8
60C
61C THIS IS A LOOP OVER THE 8 CORNERS.
62C FIRST FIND THE LOCAL COORDINATES.
63C
64 IF(ISH.EQ.28) THEN
65 goto 999
66C General twisted trapezoid.
67C
68* IL=(IP+1)/2
69* I0=IL*4+11
70* IS=(IP-IL*2)*2+1
71* X(3)=PARS(1)*IS
72* X(1)=PARS(I0)+PARS(I0+2)*X(3)
73* X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
74* GO TO 20
75C
76 ENDIF
77C
78 IP3=ISH+2
79 IF(ISH.EQ.10) IP3=3
80 IF(ISH.EQ.4) IP3=1
81 X(3)=PARS(IP3)
82 IF(IP.LE.4) X(3)=-X(3)
83 IP2=3
84 IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
85 IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
86 IF(ISH.EQ.4) IP2=4
87 IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
88 X(2)=PARS(IP2)
89 IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
90 IP1=1
91 IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
92 IF(ISH.EQ.4) IP1=5
93 IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
94 IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
95 X(1)=PARS(IP1)
96 IF(MOD(IP,2).EQ.1) X(1)=-X(1)
97C
98 IF(ISH.NE.10) GO TO 10
99 X(1)=X(1)+X(2)*PARS(4)+X(3)*PARS(5)
100 X(2)=X(2)+X(3)*PARS(6)
101 10 CONTINUE
102C
103 IF(ISH.NE.4) GO TO 20
104 IP4=7
105 IF(X(3).GT.0.0) IP4=11
106 X(1)=X(1)+X(2)*PARS(IP4)+X(3)*PARS(2)
107 X(2)=X(2)+X(3)*PARS(3)
108 20 CONTINUE
109C
110C ROTATE.
111C
112 JROT=LQ(JROTM-IROT)
113 XT(1)=X(1)
114 XT(2)=X(2)
115 XT(3)=X(3)
116 IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
117C
118 Z=DX(3)+XT(3)
119 R=(DX(1)+XT(1))**2+(DX(2)+XT(2))**2
120 RT=R+(DX(3)+XT(3))**2
121 IF(RT.LT.1.0E-10) GO TO 999
122C
123 IF(R.GT.0.0) GO TO 30
124 IF(Z.GT.0.0) T=0.0
125 IF(Z.LT.0.0) T=180.0
126 GO TO 40
127C
128 30 CONTINUE
129 T=90.0
130 IF(ABS(Z).LT.1.0E-6) GO TO 40
131C
132 R=SQRT(R)
133 T=ATAN(R/Z)*RADDEG
134 IF(T.EQ.0.0.AND.Z.LT.0.0) T=180.0
135 IF(T.LT.0.0) T=T+180.0
136C
137 40 CONTINUE
138 IF(T.GT.TH) TH=T
139 IF(T.LT.TL) TL=T
140C
141 50 CONTINUE
142C
143C THETA LIMITS SET FROM THE POINTS NOW DO THE EDGES.
144C
145 DO 120 IL=1,12
146C
147C FIND THE END POINT NUMBERS FOR EACH EDGE.
148C
149 IF(IL.GT.4) GO TO 60
150 IPP1=1
151 IF(IL.GT.2) IPP1=4
152 IPP2=2
153 IF(MOD(IL,2).EQ.1) IPP2=3
154C
155 GO TO 80
156 60 CONTINUE
157 IF(IL.LT.9) GO TO 70
158C
159 IPP1=5
160 IF(IL.GT.10) IPP1=8
161 IPP2=6
162 IF(MOD(IL,2).EQ.1) IPP2=7
163C
164 GO TO 80
165 70 CONTINUE
166C
167 IPP1=IL-4
168 IPP2=IL
169C
170 80 CONTINUE
171C
172C NOW GET THE POINTS AND ROTATE THEM.
173C
174 IF(ISH.EQ.28) THEN
175C
176C General twisted trapezoid.
177C
178 ILL=IPP1
179 IF(IPP1.EQ.3) ILL=4
180 IF(IPP1.EQ.4) ILL=3
181 I0=ILL*4+11
182 X1(3)=PARS(1)
183 IF(IPP1.LT.5) X1(3)=-X1(3)
184 X1(1)=PARS(I0)+PARS(I0+2)*X1(3)
185 X1(2)=PARS(I0+1)+PARS(I0+3)*X1(3)
186 ILL=IPP2
187 IF(IPP2.EQ.3) ILL=4
188 IF(IPP2.EQ.4) ILL=3
189 I0=ILL*4+11
190 X2(3)=PARS(1)
191 IF(IPP2.LT.5) X2(3)=-X2(3)
192 X2(1)=PARS(I0)+PARS(I0+2)*X2(3)
193 X2(2)=PARS(I0+1)+PARS(I0+3)*X2(3)
194C
195 GO TO 100
196C
197 ENDIF
198C
199 IP3=ISH+2
200 IF(ISH.EQ.10) IP3=3
201 IF(ISH.EQ.4) IP3=1
202 X1(3)=PARS(IP3)
203 IF(IPP1.LE.4) X1(3)=-X1(3)
204 IP2=3
205 IF(ISH.GT.2.AND.X1(3).GT.0.0) IP2=4
206 IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
207 IF(ISH.EQ.4) IP2=4
208 IF(ISH.EQ.4.AND.X1(3).GT.0.0) IP2=8
209 X1(2)=PARS(IP2)
210 IF(MOD(IPP1+3,4).LT.2) X1(2)=-X1(2)
211 IP1=1
212 IF(ISH.NE.1.AND.ISH.NE.10.AND.X1(3).GT.0.0) IP1=2
213 IF(ISH.EQ.4) IP1=5
214 IF(ISH.EQ.4.AND.X1(3).GT.0.0) IP1=IP1+4
215 IF(ISH.EQ.4.AND.X1(2).GT.0.0) IP1=IP1+1
216 X1(1)=PARS(IP1)
217 IF(MOD(IPP1,2).EQ.1) X1(1)=-X1(1)
218C
219 IP3=ISH+2
220 IF(ISH.EQ.10) IP3=3
221 IF(ISH.EQ.4) IP3=1
222 X2(3)=PARS(IP3)
223 IF(IPP2.LE.4) X2(3)=-X2(3)
224 IP2=3
225 IF(ISH.GT.2.AND.X2(3).GT.0.0) IP2=4
226 IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
227 IF(ISH.EQ.4) IP2=4
228 IF(ISH.EQ.4.AND.X2(3).GT.0.0) IP2=8
229 X2(2)=PARS(IP2)
230 IF(MOD(IPP2+3,4).LT.2) X2(2)=-X2(2)
231 IP1=1
232 IF(ISH.NE.1.AND.ISH.NE.10.AND.X2(3).GT.0.0) IP1=2
233 IF(ISH.EQ.4) IP1=5
234 IF(ISH.EQ.4.AND.X2(3).GT.0.0) IP1=IP1+4
235 IF(ISH.EQ.4.AND.X2(2).GT.0.0) IP1=IP1+1
236 X2(1)=PARS(IP1)
237 IF(MOD(IPP2,2).EQ.1) X2(1)=-X2(1)
238C
239 IF(ISH.NE.10) GO TO 90
240 X1(1)=X1(1)+X1(2)*PARS(4)+X1(3)*PARS(5)
241 X1(2)=X1(2)+X1(3)*PARS(6)
242 X2(1)=X2(1)+X2(2)*PARS(4)+X2(3)*PARS(5)
243 X2(2)=X2(2)+X2(3)*PARS(6)
244 90 CONTINUE
245C
246 IF(ISH.NE.4) GO TO 100
247 IP4=7
248 IF(X1(3).GT.0.0) IP4=11
249 X1(1)=X1(1)+X1(2)*PARS(IP4)+X1(3)*PARS(2)
250 X1(2)=X1(2)+X1(3)*PARS(3)
251 IP4=7
252 IF(X2(3).GT.0.0) IP4=11
253 X2(1)=X2(1)+X2(2)*PARS(IP4)+X2(3)*PARS(2)
254 X2(2)=X2(2)+X2(3)*PARS(3)
255 100 CONTINUE
256C
257C ROTATE.
258C
259 JROT=LQ(JROTM-IROT)
260 XT1(1)=X1(1)
261 XT1(2)=X1(2)
262 XT1(3)=X1(3)
263 IF(IROT.NE.0) CALL GINROT(X1,Q(JROT+1),XT1)
264 XT2(1)=X2(1)
265 XT2(2)=X2(2)
266 XT2(3)=X2(3)
267 IF(IROT.NE.0) CALL GINROT(X2,Q(JROT+1),XT2)
268C
269C NOW WE HAVE THE END POINTS IN THE MASTER SYSTEM.
270C FIND THE TANGENT POINT TO THE SET OF CONES OF CONSTANT
271C THETA.
272C
273 DS=(XT2(1)-XT1(1))**2+(XT2(2)-XT1(2))**2+(XT2(3)-XT1(3))**2
274 IF(DS.LE.0.0) GO TO 120
275C
276 DS=SQRT(DS)
277C
278 X0=(XT2(1)+XT1(1))/2.0+DX(1)
279 Y0=(XT2(2)+XT1(2))/2.0+DX(2)
280 Z0=(XT2(3)+XT1(3))/2.0+DX(3)
281 ALX=(XT2(1)-XT1(1))/DS
282 ALY=(XT2(2)-XT1(2))/DS
283 ALZ=(XT2(3)-XT1(3))/DS
284C
285 T=90.0
286 IF(Z0.EQ.0.0.AND.ALZ.EQ.0.0) GO TO 110
287C
288 IF(ALX.EQ.0.0.AND.ALY.EQ.0.0) GO TO 120
289C THIS CHECKS WHETHER THE LINE IS PARALLEL TO THE
290C Z AXIS IN WHICH CASE THERE IS NO TANGENT AND
291C THE END POINTS DETERMINE THE THETA RANGE.
292C
293 SNUM=(X0*Z0*ALX+Y0*Z0*ALY+X0*X0*ALZ+Y0*Y0*ALZ)
294 SDEN=(Z0*ALX*ALX-X0*ALX*ALZ+Z0*ALY*ALY-Y0*ALY*ALZ)
295C
296 IF(ABS(SNUM).GT.0.5*DS*ABS(SDEN)) GO TO 120
297C
298C TANGENT EXIST BETWEEN THE TWO ENDS.
299C
300 S = SNUM/SDEN
301 X0=X0+S*ALX
302 Y0=Y0+S*ALY
303 Z0=Z0+S*ALZ
304C
305 R=X0*X0+Y0*Y0
306 RT=R+Z0*Z0
307C
308 IF(RT.LT.1.0E-10) GO TO 999
309C
310 IF(R.GT.0.0) R=SQRT(R)
311 T=90.0
312 IF(ABS(Z0).LT.1.0E-06) GO TO 110
313C
314 T=ATAN(R/Z0)*RADDEG
315 IF(T.EQ.0.0.AND.Z0.LT.0.0) T=180.0
316 IF(T.LT.0.0) T=T+180.0
317C
318 110 CONTINUE
319C
320 IF(T.LT.TL) TL=T
321 IF(T.GT.TH) TH=T
322C
323C CHECK FOR THE POSSIBILITY OF STRADDLING T=0.0 &/OR 180.0.
324C
325 C=X0*DX(1)+Y0*DX(2)
326 IF(C.GT.0.0) GO TO 120
327C
328C CHECK IF SAME SIGN OF Z.
329C
330 IF(Z0*DX(3).LT.0.0) GO TO 999
331C
332 T=0.0
333 IF(Z0.LT.0.0) T=180.0
334 IF(T.LT.TL) TL=T
335 IF(T.GT.TH) TH=T
336C
337 120 CONTINUE
338C
339C DONE SET THE LIMITS.
340C
341 CL=TL
342 CH=TH
343C
344 GO TO 999
345C
346 130 CONTINUE
347C
348C TUBES, CONES ETC.
349C
350 IF(IROT.NE.0.OR.DXS.GT.1.0E-05) GO TO 180
351C
352C UNROTATED AND CENTERED.
353C
354 IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13
355 + .AND.ISH.NE.14)GO TO 170
356C
357C TUBES AND CONES.
358C
359 IERR=0
360 DZ=PARS(3)
361 RMN=PARS(1)
362 RMX=PARS(2)
363 IF(ISH.EQ.13) THEN
364**
365** approxime to a cylinder whit radius
366** equal to the ellipse major axis
367**
368 RMN=0.0
369 IF(PARS(1).GT.RMX) RMX=PARS(1)
370 GOTO 140
371 ENDIF
372 IF(ISH.EQ.14) THEN
373C not really sure of the function of these... keep RMN=PARS(1)
374 RMX = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
375 GO TO 140
376 ENDIF
377 IF(ISH.EQ.NSCTUB) THEN
378 S1 = (1.0-PARS(8))*(1.0+PARS(8))
379 IF( S1 .GT. 0.0) S1 = SQRT(S1)
380 S2 = (1.0-PARS(11))*(1.0+PARS(11))
381 IF( S2 .GT. 0.0) S2 = SQRT(S2)
382 IF( S2 .GT. S1 ) S1 = S2
383 DZ = DZ+RMX*S1
384 ENDIF
385 IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 140
386 DZ=PARS(1)
387 RMN=PARS(2)
388 RMX=PARS(3)
389 IF(PARS(4).LT.RMN) RMN=PARS(4)
390 IF(PARS(5).GT.RMX) RMX=PARS(5)
391 140 CONTINUE
392C
393 IF(DZ.GT.ABS(DX(3))) GO TO 160
394C
395C ALL FORWARD OR ALL BACK.
396C
397 DSH=DX(3)-DZ
398 DLN=DX(3)+DZ
399 IF(DX(3).GT.0.0) GO TO 150
400 DSS=DSH
401 DSH=DLN
402 DLN=DSS
403 150 CONTINUE
404C
405 CL=90.0
406 CH=90.0
407 IF(DLN.NE.0.0) CL=ATAN(RMN/DLN)*RADDEG
408 IF(DSH.NE.0.0) CH=ATAN(RMX/DSH)*RADDEG
409 IF(DX(3).GT.0.0) GO TO 999
410 CS=CL
411 CL=CH
412 CH=CS
413 IF(CH.EQ.0.0) CH=180.0
414 IF(CL.LT.0.0) CL=CL+180.0
415 IF(CH.LT.0.0) CH=CH+180.0
416C
417 GO TO 999
418C
419 160 CONTINUE
420C
421 CL=90.0
422 CH=90.0
423 IF(DZ+DX(3).NE.0.0) CL=ATAN(RMN/(DZ+DX(3)))*RADDEG
424 IF(-DZ+DX(3).NE.0.0) CH=ATAN(RMN/(-DZ+DX(3)))*RADDEG
425 IF(CH.LE.0.0) CH=CH+180.0
426C
427 GO TO 999
428C
429 170 CONTINUE
430 IF(ISH.GT.9) GO TO 999
431C
432C SPHERE.
433C
434 IERR=0
435 CL=PARS(3)
436 CH=PARS(4)
437C
438 GO TO 999
439 180 CONTINUE
440C
441 IF(ISH.EQ.11.OR.ISH.EQ.12) GOTO 999
442**
443 RM=PARS(2)
444 IF(ISH.EQ.9) GO TO 200
445C
446 DZ=PARS(3)
447 IF(ISH.EQ.NSCTUB) THEN
448 S1 = (1.0-PARS(8))*(1.0+PARS(8))
449 IF( S1 .GT. 0.0) S1 = SQRT(S1)
450 S2 = (1.0-PARS(11))*(1.0+PARS(11))
451 IF( S2 .GT. 0.0) S2 = SQRT(S2)
452 IF( S2 .GT. S1 ) S1 = S2
453 DZ = DZ+RM*S1
454 ENDIF
455**
456 IF(ISH.EQ.13) THEN
457 IF(PARS(1).GT.RM) RM=PARS(1)
458 GOTO 190
459 ENDIF
460**
461 IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 190
462C
463 DZ=PARS(1)
464 RM=PARS(3)
465 IF(PARS(5).GT.RM) RM=PARS(5)
466C
467 190 CONTINUE
468C
469 RM=SQRT(RM**2+DZ**2)
470C
471 200 CONTINUE
472C
473 CL=0.0
474 CH=180.0
475 IERR=0
476 RC=DXS**2+DX(3)**2
477 IF(RC.GT.0.0) RC=SQRT(RC)
478 IF(RM.GE.RC) GO TO 999
479C
480 TC=90.0
481 IF(ABS(DX(3)).GT.0.0) TC=ATAN(DXS/DX(3))*RADDEG
482 IF(TC.LT.0.0) TC=TC+180.0
483C
484 DT=ASIN(RM/RC)*RADDEG
485 CL=TC-DT
486 CH=TC+DT
487 IF(CL.LT.0.0) CL=0.0
488 IF(CH.GT.180.0) CH=180.0
489C
490 999 CONTINUE
491 END