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