This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gflphi.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 GFLPHI(ISH,IROT,DX,PARS,CL,CH,IERR)
13 C.
14 C.    **********************************************************
15 C.    *                                                        *
16 C.    *    ROUTINE TO FIND THE PHI LIMITS OF THE OBJECT SHAPE  *
17 C.    *    ISH ROTATED BY MATRIX IROT AND DISPLACED BY VECTOR  *
18 C.    *    DX. IT HAS NPAR PARAMTERS IN THE ARRAY PARS. THE    *
19 C.    *    LOWER LIMIT IS RETURNED IN CL AND THE HIGHER IN CH. *
20 C.    *    NOTE THE OBJECT IS CONTAINED IN THE RANGE OF        *
21 C.    *    INCREASING PHI FROM CL TO CH THOUGH CL AND CH ARE   *
22 C.    *    FORCED TO LIE IN THE RANGE 0.0 TO 360.0 SO THAT THE *
23 C.    *    VALUE OF CL CAN BE HIGHER THAN THAT OF CH.          *
24 C.    *                                                        *
25 C.    *    ==>Called by : GFCLIM                               *
26 C.    *         Author  A.McPherson  *********                 *
27 C.    *                                                        *
28 C.    **********************************************************
29 C.
30 #include "geant321/gcbank.inc"
31 #include "geant321/gconsp.inc"
32 #include "geant321/gcshno.inc"
33       DIMENSION DX(3),PARS(11),X(3),XT(3)
34 C.
35 C.          -------------------------------------------
36 C.
37       IERR=1
38 C
39       DXS=DX(1)*DX(1)+DX(2)*DX(2)
40       IF(DXS.GT.0.0) DXS=SQRT(DXS)
41 C
42       IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 40
43 C
44 C
45 C           CUBOIDS, TRAPEZOIDS, PARALLELEPIPEDS.
46 C
47       IERR=0
48       CL=0.0
49       CH=360.0
50 C
51 C           IF IN DOUBT SET IT TO FULL RANGE.
52 C
53       IF(DXS.LE.0.0) GO TO 999
54 C
55       PHC=90.
56       IF(DX(1).NE.0.)PHC=ATAN2(DX(2),DX(1))*RADDEG
57       IF(PHC.LT.0.0) PHC=PHC+360
58       PL=0.0
59       PH=0.0
60 C
61       DO 30 IP=1,8
62 C
63 C           THIS IS A LOOP OVER THE 8 CORNERS.
64 C           FIRST FIND THE LOCAL COORDINATES.
65 C
66       IF(ISH.EQ.28) THEN
67 C
68 C            General twisted trapezoid.
69 C
70          IL=(IP+1)/2
71          I0=IL*4+11
72          IS=(IP-IL*2)*2+1
73          X(3)=PARS(1)*IS
74          X(1)=PARS(I0)+PARS(I0+2)*X(3)
75          X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
76          GO TO 20
77 C
78       ENDIF
79 C
80       IP3=ISH+2
81       IF(ISH.EQ.10) IP3=3
82       IF(ISH.EQ.4) IP3=1
83       X(3)=PARS(IP3)
84       IF(IP.LE.4) X(3)=-X(3)
85       IP2=3
86       IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
87       IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
88       IF(ISH.EQ.4) IP2=4
89       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
90       X(2)=PARS(IP2)
91       IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
92       IP1=1
93       IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
94       IF(ISH.EQ.4) IP1=5
95       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
96       IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
97       X(1)=PARS(IP1)
98       IF(MOD(IP,2).EQ.1) X(1)=-X(1)
99 C
100       IF(ISH.NE.10) GO TO 10
101       X(1)=X(1)+X(2)*PARS(4)+X(3)*PARS(5)
102       X(2)=X(2)+X(3)*PARS(6)
103    10 CONTINUE
104 C
105       IF(ISH.NE.4) GO TO 20
106       IP4=7
107       IF(X(3).GT.0.0) IP4=11
108       X(1)=X(1)+X(2)*PARS(IP4)+X(3)*PARS(2)
109       X(2)=X(2)+X(3)*PARS(3)
110    20 CONTINUE
111 C
112 C          ROTATE.
113 C
114       JROT=LQ(JROTM-IROT)
115       XT(1)=X(1)
116       XT(2)=X(2)
117       XT(3)=X(3)
118       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
119 C
120       XPT=DXS+(DX(1)*XT(1)+DX(2)*XT(2))/DXS
121       YPT=(DX(1)*XT(2)-DX(2)*XT(1))/DXS
122 C
123       IF(YPT.EQ.0.0.AND.XPT.EQ.0.0) GO TO 999
124       P=ATAN2(YPT,XPT)
125       IF(P.GT.PI) P=P-PI*2.0
126       IF(P.LT.PL) PL=P
127       IF(P.GT.PH) PH=P
128 C
129 C
130    30 CONTINUE
131 C
132 C
133       IF(PH-PL.GT.PI) GO TO 999
134       CL=PHC+PL*RADDEG
135       CH=PHC+PH*RADDEG
136 C
137       SG = SIGN(1.0,CL)
138       CL = MOD( ABS(CL),360.0 )
139       IF(SG.LE.0.0) CL=360.-CL
140       SG=SIGN(1.0,CH)
141       CH = MOD( ABS(CH),360.0 )
142       IF(SG.LE.0.0) CH=360.-CH
143 C
144       GO TO 999
145 C
146    40 CONTINUE
147       IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14)GO TO 80
148 C
149 C             TUBES AND CONES.
150 C
151       IERR=0
152       CL=0.0
153       CH=360.0
154 C
155 C             WHEN IN DOUBT SET TO FULL RANGE.
156 C
157       RM=PARS(2)
158       IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 50
159 **
160       IF(ISH.EQ.13) THEN
161 **
162 **       approxime to a cylinder whit radius
163 **       equal to the ellipse major axis
164 **
165          IF(PARS(1).GT.RM) RM=PARS(1)
166          GOTO 50
167       ENDIF
168       IF(ISH.EQ.14) THEN
169         RM = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
170         GO TO 50
171       ENDIF
172       RM=PARS(3)
173       IF(PARS(5).GT.RM) RM=PARS(5)
174 C
175    50 CONTINUE
176 C
177       IF(DXS.GT.RM) GO TO 70
178       IF(ISH.EQ.5.OR.ISH.EQ.7.OR.ISH.EQ.14) GO TO 999
179       IF(ISH.EQ.13) GOTO 999
180 *                 Here we treat the CONS
181       IF(IROT.EQ.0) THEN
182 *                 This is the simple case, no rotation
183 *                 Compute the position of the limits on
184 *                 the X-Y plane.
185       PHIMIN=PARS(4)*DEGRAD
186       PHIMAX=PARS(5)*DEGRAD
187       DDX1 = DX(1)+RM*COS(PHIMIN)
188       DDY1 = DX(2)+RM*SIN(PHIMIN)
189       DDX2 = DX(1)+RM*COS(PHIMAX)
190       DDY2 = DX(2)+RM*SIN(PHIMAX)
191       CL = ATAN2(DDY1,DDX1)*RADDEG
192       CH = ATAN2(DDY2,DDX2)*RADDEG
193       ELSE
194 *                 Rotated tubes might be more difficult
195 *                 Just leave it for later
196       CONTINUE
197       ENDIF
198       IF(ISH.LE.6) GO TO 60
199       CL=PARS(6)
200       CH=PARS(7)
201    60 CONTINUE
202 C
203       IF(IROT.EQ.0) GO TO 65
204       JROT=LQ(JROTM-IROT)
205       IF(Q(JROT+15).NE.0.0.AND.Q(JROT+15).NE.180.0) GO TO 999
206 C
207       PHX=Q(JROT+12)
208       PHY=Q(JROT+14)
209       IF(PHY.LT.PHX) PHY=PHY+360.0
210       ISPH=1
211       IF(PHY-PHX.GT.180.0) ISPH=-1
212       CL=ISPH*CL+PHX
213       CH=ISPH*CH+PHX
214       IF(ISPH.EQ.1) GO TO 65
215       CHT=CH
216       CH=CL
217       CL=CHT
218 C
219    65 CONTINUE
220 C
221       SG=SIGN(1.0,CL)
222       CL = MOD( ABS(CL),360.0 )
223       IF(SG.LE.0.0) CL=360.-CL
224       SG=SIGN(1.0,CH)
225       CH = MOD( ABS(CH),360.0 )
226       IF(SG.LE.0.0) CH=360.-CH
227 C
228       GO TO 999
229 C
230    70 CONTINUE
231 C
232 C            DISPLACEMENT GREATER THAN MAXIMUM RADIUS SO
233 C            ASSUME COMPLETE TUBE TO GENERATE 'WORST CASE'.
234 C
235       DZ=PARS(3)
236       IF(ISH.EQ.NSCTUB) THEN
237         S1 = (1.0-PARS(8))*(1.0+PARS(8))
238         IF( S1 .GT. 0.0) S1 = SQRT(S1)
239         S2 = (1.0-PARS(11))*(1.0+PARS(11))
240         IF( S2 .GT. 0.0) S2 = SQRT(S2)
241         IF( S2 .GT. S1 ) S1 = S2
242         DZ = DZ+RM*S1
243       ELSEIF(ISH.GT.6.AND.ISH.NE.13.AND.ISH.NE.14) THEN
244         DZ=PARS(1)
245       ENDIF
246 C
247       X(1)=0.0
248       X(2)=0.0
249       X(3)=1.0
250 C
251 C                    LOCAL Z AXIS.
252 C
253       JROT=LQ(JROTM-IROT)
254       XT(1)=X(1)
255       XT(2)=X(2)
256       XT(3)=X(3)
257       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
258 C
259       COST=ABS(DX(1)*XT(1)+DX(2)*XT(2))
260       COS2=ABS(DX(1)*XT(2)-DX(2)*XT(1))
261       SINT=(DXS+COST)*(DXS-COST)
262       SIN2=(DXS+COS2)*(DXS-COS2)
263       IF(SINT.GT.0.0) SINT=SQRT(SINT)
264       IF(SIN2.GT.0.0) SIN2=SQRT(SIN2)
265 C
266       XPT=DXS-(COST*DZ+SINT*RM)/DXS
267 C
268       IF(XPT.LE.0.0) GO TO 999
269       YPT=(SIN2*RM+COS2*DZ)/DXS
270       DP=ATAN(YPT/XPT)
271 C
272       P0=ATAN2(DX(2),DX(1))
273       CL=(P0-DP)*RADDEG
274       CH=(P0+DP)*RADDEG
275 C
276       SG=SIGN(1.0,CL)
277       CL = MOD( ABS(CL),360.0 )
278       IF(SG.LE.0.0) CL=360.-CL
279       SG=SIGN(1.0,CH)
280       CH = MOD( ABS(CH),360.0 )
281       IF(SG.LE.0.0) CH=360.-CH
282 C
283       GO TO 999
284 C
285    80 CONTINUE
286       IF(ISH.GT.9) GO TO 999
287 C
288 C               SPHERE.
289 C
290       IERR=0
291       CL=0.0
292       CH=360.0
293 C
294       IF(IROT.NE.0.OR.DXS.GT.0.0) GO TO 90
295 C
296 C          UNROTATED AND CENTERED.
297 C
298       CL=PARS(5)
299       CH=PARS(6)
300 C
301       SG=SIGN(1.0,CL)
302       CL = MOD( ABS(CL),360.0 )
303       IF(SG.LE.0.0) CL=360.-CL
304       SG=SIGN(1.0,CH)
305       CH = MOD( ABS(CH),360.0 )
306       IF(SG.LE.0.0) CH=360.-CH
307 C
308       GO TO 999
309 C
310    90 CONTINUE
311 C
312 C            ROTATED OR NOT CENTERED.
313 C
314       IF(DXS.LT.PARS(2)) GO TO 999
315       P0=ATAN2(DX(2),DX(1))
316       DP=ASIN(PARS(2)/DXS)
317       CL=(P0-DP)*RADDEG
318       CH=(P0+DP)*RADDEG
319 C
320       SG=SIGN(1.0,CL)
321       CL = MOD( ABS(CL),360.0 )
322       IF(SG.LE.0.0) CL=360.-CL
323       SG=SIGN(1.0,CH)
324       CH = MOD( ABS(CH),360.0 )
325       IF(SG.LE.0.0) CH=360.-CH
326 C
327   999 CONTINUE
328       END