Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gvdthe.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1999/05/18 15:55:17  fca
6 * AliRoot sources
7 *
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)
16 C.
17 C.    ******************************************************************
18 C.    *                                                                *
19 C.    *    ROUTINE TO COMPUTE THE THETA LIMITS FOR VOLUME OF SHAPE     *
20 C.    *    ISH ROTATED BY MATRIX IROT AND DISPLACED BY VECTOR DX.      *
21 C.    *    THE VOLUME HAS NPAR PARAMETERS IN THE ARRAY PARS. THE LOWER *
22 C.    *    LIMIT IS RETURNED IN CL THE HIGHER IN CH. IF THE            *
23 C.    *    CALCULATION CANNOT BE MADE IERR IS SET TO 1 OTHERWISE IT    *
24 C.    *    IS SET TO 0.                                                *
25 C.    *                                                                *
26 C.    *    ==>Called by : GVDLIM                                       *
27 C.    *         Author  S.Giani  *********                             *
28 C.    *                                                                *
29 C.    ******************************************************************
30 C.
31 #include "geant321/gcbank.inc"
32 #include "geant321/gconsp.inc"
33 #include "geant321/gcshno.inc"
34       DIMENSION DX(3),PARS(100),X(3),XT(3),X1(3),X2(3),XT1(3),XT2(3)
35 C.
36 C.          ----------------------------------------------
37 C.
38       IERR=1
39 C
40       DXS=DX(1)*DX(1)+DX(2)*DX(2)
41       IF(DXS.GT.0.0) DXS=SQRT(DXS)
42 C
43       IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 130
44 C
45 C            CUBES, TRAPEZOIDS AND PARALLELEPIPEDS.
46 C
47       IERR=0
48       CL=0.0
49       CH=180.0
50 C
51       IF(DXS.LE.0.0) GO TO 999
52 C
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
58 C
59       DO 50 IP=1,8
60 C
61 C           THIS IS A LOOP OVER THE 8 CORNERS.
62 C           FIRST FIND THE LOCAL COORDINATES.
63 C
64       IF(ISH.EQ.28) THEN
65          goto 999
66 C            General twisted trapezoid.
67 C
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
75 C
76       ENDIF
77 C
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)
97 C
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
102 C
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
109 C
110 C          ROTATE.
111 C
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)
117 C
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
122 C
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
127 C
128    30 CONTINUE
129       T=90.0
130       IF(ABS(Z).LT.1.0E-6) GO TO 40
131 C
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
136 C
137    40 CONTINUE
138       IF(T.GT.TH) TH=T
139       IF(T.LT.TL) TL=T
140 C
141    50 CONTINUE
142 C
143 C           THETA LIMITS SET FROM THE POINTS NOW DO THE EDGES.
144 C
145       DO 120 IL=1,12
146 C
147 C           FIND THE END POINT NUMBERS FOR EACH EDGE.
148 C
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
154 C
155       GO TO 80
156    60 CONTINUE
157       IF(IL.LT.9) GO TO 70
158 C
159       IPP1=5
160       IF(IL.GT.10) IPP1=8
161       IPP2=6
162       IF(MOD(IL,2).EQ.1) IPP2=7
163 C
164       GO TO 80
165    70 CONTINUE
166 C
167       IPP1=IL-4
168       IPP2=IL
169 C
170    80 CONTINUE
171 C
172 C           NOW GET THE POINTS AND ROTATE THEM.
173 C
174       IF(ISH.EQ.28) THEN
175 C
176 C            General twisted trapezoid.
177 C
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)
194 C
195          GO TO 100
196 C
197       ENDIF
198 C
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)
218 C
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)
238 C
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
245 C
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
256 C
257 C          ROTATE.
258 C
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)
268 C
269 C           NOW WE HAVE THE END POINTS IN THE MASTER SYSTEM.
270 C           FIND THE TANGENT POINT TO THE SET OF CONES OF CONSTANT
271 C           THETA.
272 C
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
275 C
276       DS=SQRT(DS)
277 C
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
284 C
285       T=90.0
286       IF(Z0.EQ.0.0.AND.ALZ.EQ.0.0) GO TO 110
287 C
288       IF(ALX.EQ.0.0.AND.ALY.EQ.0.0) GO TO 120
289 C             THIS CHECKS WHETHER THE LINE IS PARALLEL TO THE
290 C             Z AXIS IN WHICH CASE THERE IS NO TANGENT AND
291 C             THE END POINTS DETERMINE THE THETA RANGE.
292 C
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)
295 C
296       IF(ABS(SNUM).GT.0.5*DS*ABS(SDEN)) GO TO 120
297 C
298 C           TANGENT EXIST BETWEEN THE TWO ENDS.
299 C
300       S = SNUM/SDEN
301       X0=X0+S*ALX
302       Y0=Y0+S*ALY
303       Z0=Z0+S*ALZ
304 C
305       R=X0*X0+Y0*Y0
306       RT=R+Z0*Z0
307 C
308       IF(RT.LT.1.0E-10) GO TO 999
309 C
310       IF(R.GT.0.0) R=SQRT(R)
311       T=90.0
312       IF(ABS(Z0).LT.1.0E-06) GO TO 110
313 C
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
317 C
318   110 CONTINUE
319 C
320       IF(T.LT.TL) TL=T
321       IF(T.GT.TH) TH=T
322 C
323 C        CHECK FOR THE POSSIBILITY OF STRADDLING T=0.0 &/OR 180.0.
324 C
325       C=X0*DX(1)+Y0*DX(2)
326       IF(C.GT.0.0) GO TO 120
327 C
328 C          CHECK IF SAME SIGN OF Z.
329 C
330       IF(Z0*DX(3).LT.0.0) GO TO 999
331 C
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
336 C
337   120 CONTINUE
338 C
339 C            DONE SET THE LIMITS.
340 C
341       CL=TL
342       CH=TH
343 C
344       GO TO 999
345 C
346   130 CONTINUE
347 C
348 C          TUBES, CONES ETC.
349 C
350       IF(IROT.NE.0.OR.DXS.GT.1.0E-05) GO TO 180
351 C
352 C          UNROTATED AND CENTERED.
353 C
354       IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13
355      +           .AND.ISH.NE.14)GO TO 170
356 C
357 C               TUBES AND CONES.
358 C
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
373 C   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
392 C
393       IF(DZ.GT.ABS(DX(3))) GO TO 160
394 C
395 C          ALL FORWARD OR ALL BACK.
396 C
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
404 C
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
416 C
417       GO TO 999
418 C
419   160 CONTINUE
420 C
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
426 C
427       GO TO 999
428 C
429   170 CONTINUE
430       IF(ISH.GT.9) GO TO 999
431 C
432 C           SPHERE.
433 C
434       IERR=0
435       CL=PARS(3)
436       CH=PARS(4)
437 C
438       GO TO 999
439   180 CONTINUE
440 C
441       IF(ISH.EQ.11.OR.ISH.EQ.12) GOTO 999
442 **
443       RM=PARS(2)
444       IF(ISH.EQ.9) GO TO 200
445 C
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
462 C
463       DZ=PARS(1)
464       RM=PARS(3)
465       IF(PARS(5).GT.RM) RM=PARS(5)
466 C
467   190 CONTINUE
468 C
469       RM=SQRT(RM**2+DZ**2)
470 C
471   200 CONTINUE
472 C
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
479 C
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
483 C
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
489 C
490   999 CONTINUE
491       END