This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gvdphi.F
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)
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 : GVDLIM modified                      *
26 C.    *         Author  S.Giani  *********                     *
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(50),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(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
61 C
62       DO 30 IP=1,8
63 C
64 C           THIS IS A LOOP OVER THE 8 CORNERS.
65 C           FIRST FIND THE LOCAL COORDINATES.
66 C
67       IF(ISH.EQ.28) THEN
68 C
69 C            General twisted trapezoid.
70 C
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
78 C
79       ENDIF
80 C
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)
100 C
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
105 C
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
112 C
113 C          ROTATE.
114 C
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)
120 C
121       XPT=DXS+(DX(1)*XT(1)+DX(2)*XT(2))/DXS
122       YPT=(DX(1)*XT(2)-DX(2)*XT(1))/DXS
123 C
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
129 C
130 C
131    30 CONTINUE
132 C
133 C
134       IF(PH-PL.GT.PI) GO TO 999
135       CL=PHC+PL*RADDEG
136       CH=PHC+PH*RADDEG
137 C
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
144 C
145       GO TO 999
146 C
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
206 C
207 C             TUBES AND CONES.
208 C
209       IERR=0
210       CL=0.0
211       CH=360.0
212 C
213 C             WHEN IN DOUBT SET TO FULL RANGE.
214 C
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)
232 C
233    50 CONTINUE
234 C
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
280 C
281    60 CONTINUE
282 C
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
291 C
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
322 C
323    65 CONTINUE
324 C
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
331 C
332       GO TO 999
333 C
334    70 CONTINUE
335 C
336 C            DISPLACEMENT GREATER THAN MAXIMUM RADIUS SO
337 C            ASSUME COMPLETE TUBE TO GENERATE 'WORST CASE'.
338 C
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
350 C
351       X(1)=0.0
352       X(2)=0.0
353       X(3)=1.0
354 C
355 C                    LOCAL Z AXIS.
356 C
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)
362 C
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)
369 C
370       XPT=DXS-(COST*DZ+SINT*RM)/DXS
371 C
372       IF(XPT.LE.0.0) GO TO 999
373       YPT=(SIN2*RM+COS2*DZ)/DXS
374       DP=ATAN(YPT/XPT)
375 C
376       P0=ATAN2(DX(2),DX(1))
377       CL=(P0-DP)*RADDEG
378       CH=(P0+DP)*RADDEG
379 C
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
386 C
387       GO TO 999
388 C
389    80 CONTINUE
390       IF(ISH.GT.9.AND.MYFLAG.EQ.0) GO TO 999
391 C
392 C               SPHERE.
393 C
394       IERR=0
395       CL=0.0
396       CH=360.0
397 C
398       IF(IROT.NE.0.OR.DXS.GT.0.0) GO TO 90
399 C
400 C          UNROTATED AND CENTERED.
401 C
402       CL=PARS(5)
403       CH=PARS(6)
404 C
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
411 C
412       GO TO 999
413 C
414    90 CONTINUE
415 C
416 C            ROTATED OR NOT CENTERED.
417 C
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
423 C
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
430 C
431   999 CONTINUE
432       END