5 * Revision 1.1.1.1 1999/05/18 15:55:17 fca
8 * Revision 1.1.1.1 1995/10/24 10:20:57 cernlib
12 #include "geant321/pilot.h"
13 *CMZ : 3.21/04 02/03/95 11.58.00 by S.Giani
15 SUBROUTINE GVDPHI(ISH,IROT,DX,PARS,CL,CH,IERR)
17 C. **********************************************************
19 C. * ROUTINE TO FIND THE PHI LIMITS OF THE OBJECT SHAPE *
20 C. * ISH ROTATED BY MATRIX IROT AND DISPLACED BY VECTOR *
21 C. * DX. IT HAS NPAR PARAMTERS IN THE ARRAY PARS. THE *
22 C. * LOWER LIMIT IS RETURNED IN CL AND THE HIGHER IN CH. *
23 C. * NOTE THE OBJECT IS CONTAINED IN THE RANGE OF *
24 C. * INCREASING PHI FROM CL TO CH THOUGH CL AND CH ARE *
25 C. * FORCED TO LIE IN THE RANGE 0.0 TO 360.0 SO THAT THE *
26 C. * VALUE OF CL CAN BE HIGHER THAN THAT OF CH. *
28 C. * ==>Called by : GVDLIM modified *
29 C. * Author S.Giani ********* *
31 C. **********************************************************
33 #include "geant321/gcbank.inc"
34 #include "geant321/gconsp.inc"
35 #include "geant321/gcshno.inc"
36 DIMENSION DX(3),PARS(100),X(3),XT(3)
38 C. -------------------------------------------
42 DXS=DX(1)*DX(1)+DX(2)*DX(2)
43 IF(DXS.GT.0.0) DXS=SQRT(DXS)
45 IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 40
48 C CUBOIDS, TRAPEZOIDS, PARALLELEPIPEDS.
54 C IF IN DOUBT SET IT TO FULL RANGE.
56 IF(DXS.LE.0.0) GO TO 999
59 IF(DX(1).NE.0.)PHC=ATAN2(DX(2),DX(1))*RADDEG
60 IF(DX(1).EQ.0..AND.DX(2).LT.0.)PHC=-90.
61 IF(PHC.LT.0.0) PHC=PHC+360
67 C THIS IS A LOOP OVER THE 8 CORNERS.
68 C FIRST FIND THE LOCAL COORDINATES.
72 C General twisted trapezoid.
78 X(1)=PARS(I0)+PARS(I0+2)*X(3)
79 X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
88 IF(IP.LE.4) X(3)=-X(3)
90 IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
91 IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
93 IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
95 IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
97 IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
99 IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
100 IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
102 IF(MOD(IP,2).EQ.1) X(1)=-X(1)
104 IF(ISH.NE.10) GO TO 10
105 X(1)=X(1)+X(2)*PARS(4)+X(3)*PARS(5)
106 X(2)=X(2)+X(3)*PARS(6)
109 IF(ISH.NE.4) GO TO 20
111 IF(X(3).GT.0.0) IP4=11
112 X(1)=X(1)+X(2)*PARS(IP4)+X(3)*PARS(2)
113 X(2)=X(2)+X(3)*PARS(3)
122 IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
124 XPT=DXS+(DX(1)*XT(1)+DX(2)*XT(2))/DXS
125 YPT=(DX(1)*XT(2)-DX(2)*XT(1))/DXS
127 IF(YPT.EQ.0.0.AND.XPT.EQ.0.0) GO TO 999
129 IF(P.GT.PI) P=P-PI*2.0
137 IF(PH-PL.GT.PI) GO TO 999
141 *** SG = SIGN(1.0,CL)
142 *** CL = MOD( ABS(CL),360.0 )
143 *** IF(SG.LE.0.0) CL=360.-CL
145 *** CH = MOD( ABS(CH),360.0 )
146 *** IF(SG.LE.0.0) CH=360.-CH
161 DO 145 I=7,IZLAST+2,3
162 IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
163 IF(PARS(I-1).LT.TMPMIN)TMPMIN=PARS(I-1)
167 PHIMAX=PHIMIN+PARS(2)
168 AANG=ABS(PHIMAX-PHIMIN)
170 AATMAX=NANG*360./AANG
173 IF(ALA.GT..5)LATMAX=LATMAX+1
174 AFINV=1./COS(PI/LATMAX)
177 IF(PARS(2).EQ.360)THEN
180 PHIMIN=PARS(1)*DEGRAD
181 PHIMAX=(PARS(1)+PARS(2))*DEGRAD
184 ELSEIF(ISH.EQ.12)THEN
193 DO 146 I=6,IZLAST+2,3
194 IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
195 IF(PARS(I-1).LT.TMPMIN)TMPMIN=PARS(I-1)
199 IF(PARS(2).EQ.360)THEN
202 PHIMIN=PARS(1)*DEGRAD
203 PHIMAX=(PARS(1)+PARS(2))*DEGRAD
206 ELSEIF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14)THEN
216 C WHEN IN DOUBT SET TO FULL RANGE.
218 IF(MYFLAG.EQ.0)RM=PARS(2)
219 IF(ISH.LE.6.OR.ISH.EQ.NSCTUB.OR.MYFLAG.NE.0) GO TO 50
223 ** approxime to a cylinder whit radius
224 ** equal to the ellipse major axis
226 IF(PARS(1).GT.RM) RM=PARS(1)
230 RM = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
234 IF(PARS(5).GT.RM) RM=PARS(5)
238 IF(DXS.GT.RM) GO TO 70
239 IF(ISH.EQ.5.OR.ISH.EQ.7.OR.ISH.EQ.14.OR.MYFLAG.EQ.5) GO TO 999
240 IF(ISH.EQ.13) GOTO 999
241 * Here we treat the CONS and TUBS
242 * This is the simple case, no rotation
243 * Compute the position of the limits on
245 IF(MYFLAG.EQ.0.AND.ISH.LE.6)THEN
246 PHIMIN=PARS(4)*DEGRAD
247 PHIMAX=PARS(5)*DEGRAD
249 ELSEIF(MYFLAG.EQ.0.AND.ISH.EQ.8)THEN
250 PHIMIN=PARS(6)*DEGRAD
251 PHIMAX=PARS(7)*DEGRAD
252 RLOW=MIN(PARS(2),PARS(4))
257 DDX1 = DX(1)+RM*COS(PHIMIN)
258 DDY1 = DX(2)+RM*SIN(PHIMIN)
259 DDX2 = DX(1)+RM*COS(PHIMAX)
260 DDY2 = DX(2)+RM*SIN(PHIMAX)
261 CLA = ATAN2(DDY1,DDX1)*RADDEG
262 CHA = ATAN2(DDY2,DDX2)*RADDEG
263 DDX1 = DX(1)+RLOW*COS(PHIMIN)
264 DDY1 = DX(2)+RLOW*SIN(PHIMIN)
265 DDX2 = DX(1)+RLOW*COS(PHIMAX)
266 DDY2 = DX(2)+RLOW*SIN(PHIMAX)
267 CLB = ATAN2(DDY1,DDX1)*RADDEG
268 CHB = ATAN2(DDY2,DDX2)*RADDEG
269 CL=MIN(CLA,CLB,CHA,CHB)
270 CH=MAX(CLA,CLB,CHA,CHB)
271 IF((CH-CL).GT.(PHIMAX-PHIMIN)*RADDEG)THEN
272 IF(ISH.EQ.6.OR.ISH.EQ.8.OR.MYFLAG.EQ.6)THEN
286 IF(IROT.EQ.0) GO TO 65
287 IF(CL.EQ.0..AND.CH.EQ.360.)GOTO 65
289 IF(Q(JROT+15).NE.0.0.AND.Q(JROT+15).NE.180.0)THEN
297 IF(PHY.LT.PHX) PHY=PHY+360.0
299 IF(PHY-PHX.GT.180.0) ISPH=-1
301 PHI1 = ISPH*PHIMIN+PHX*DEGRAD
302 PHI2 = ISPH*PHIMAX+PHX*DEGRAD
303 DDX1 = DX(1)+RM*COS(PHI1)
304 DDY1 = DX(2)+RM*SIN(PHI1)
305 DDX2 = DX(1)+RM*COS(PHI2)
306 DDY2 = DX(2)+RM*SIN(PHI2)
307 CLA = ATAN2(DDY1,DDX1)*RADDEG
308 CHA = ATAN2(DDY2,DDX2)*RADDEG
309 DDX1 = DX(1)+RLOW*COS(PHI1)
310 DDY1 = DX(2)+RLOW*SIN(PHI1)
311 DDX2 = DX(1)+RLOW*COS(PHI2)
312 DDY2 = DX(2)+RLOW*SIN(PHI2)
313 CLB = ATAN2(DDY1,DDX1)*RADDEG
314 CHB = ATAN2(DDY2,DDX2)*RADDEG
315 CL=MIN(CLA,CLB,CHA,CHB)
316 CH=MAX(CLA,CLB,CHA,CHB)
321 IF(ISPH.EQ.1) GO TO 65
329 *** CL = MOD( ABS(CL),360.0 )
330 *** IF(SG.LE.0.0) CL=360.-CL
332 *** CH = MOD( ABS(CH),360.0 )
333 *** IF(SG.LE.0.0) CH=360.-CH
339 C DISPLACEMENT GREATER THAN MAXIMUM RADIUS SO
340 C ASSUME COMPLETE TUBE TO GENERATE 'WORST CASE'.
342 IF(MYFLAG.EQ.0)DZ=PARS(3)
343 IF(ISH.EQ.NSCTUB) THEN
344 S1 = (1.0-PARS(8))*(1.0+PARS(8))
345 IF( S1 .GT. 0.0) S1 = SQRT(S1)
346 S2 = (1.0-PARS(11))*(1.0+PARS(11))
347 IF( S2 .GT. 0.0) S2 = SQRT(S2)
348 IF( S2 .GT. S1 ) S1 = S2
350 ELSEIF(ISH.GT.6.AND.ISH.NE.13.AND.ISH.NE.14.AND.MYFLAG.EQ.0)THEN
364 IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
366 COST=ABS(DX(1)*XT(1)+DX(2)*XT(2))
367 COS2=ABS(DX(1)*XT(2)-DX(2)*XT(1))
368 SINT=(DXS+COST)*(DXS-COST)
369 SIN2=(DXS+COS2)*(DXS-COS2)
370 IF(SINT.GT.0.0) SINT=SQRT(SINT)
371 IF(SIN2.GT.0.0) SIN2=SQRT(SIN2)
373 XPT=DXS-(COST*DZ+SINT*RM)/DXS
375 IF(XPT.LE.0.0) GO TO 999
376 YPT=(SIN2*RM+COS2*DZ)/DXS
379 P0=ATAN2(DX(2),DX(1))
384 *** CL = MOD( ABS(CL),360.0 )
385 *** IF(SG.LE.0.0) CL=360.-CL
387 *** CH = MOD( ABS(CH),360.0 )
388 *** IF(SG.LE.0.0) CH=360.-CH
393 IF(ISH.GT.9.AND.MYFLAG.EQ.0) GO TO 999
401 IF(IROT.NE.0.OR.DXS.GT.0.0) GO TO 90
403 C UNROTATED AND CENTERED.
409 *** CL = MOD( ABS(CL),360.0 )
410 *** IF(SG.LE.0.0) CL=360.-CL
412 *** CH = MOD( ABS(CH),360.0 )
413 *** IF(SG.LE.0.0) CH=360.-CH
419 C ROTATED OR NOT CENTERED.
421 IF(DXS.LT.PARS(2)) GO TO 999
422 P0=ATAN2(DX(2),DX(1))
428 *** CL = MOD( ABS(CL),360.0 )
429 *** IF(SG.LE.0.0) CL=360.-CL
431 *** CH = MOD( ABS(CH),360.0 )
432 *** IF(SG.LE.0.0) CH=360.-CH