5 * Revision 1.1.1.1 1995/10/24 10:20:48 cernlib
9 #include "geant321/pilot.h"
10 *CMZ : 3.21/02 29/03/94 15.41.28 by S.Giani
12 SUBROUTINE GFLPHI(ISH,IROT,DX,PARS,CL,CH,IERR)
14 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. *
25 C. * ==>Called by : GFCLIM *
26 C. * Author A.McPherson ********* *
28 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)
35 C. -------------------------------------------
39 DXS=DX(1)*DX(1)+DX(2)*DX(2)
40 IF(DXS.GT.0.0) DXS=SQRT(DXS)
42 IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 40
45 C CUBOIDS, TRAPEZOIDS, PARALLELEPIPEDS.
51 C IF IN DOUBT SET IT TO FULL RANGE.
53 IF(DXS.LE.0.0) GO TO 999
56 IF(DX(1).NE.0.)PHC=ATAN2(DX(2),DX(1))*RADDEG
57 IF(PHC.LT.0.0) PHC=PHC+360
63 C THIS IS A LOOP OVER THE 8 CORNERS.
64 C FIRST FIND THE LOCAL COORDINATES.
68 C General twisted trapezoid.
74 X(1)=PARS(I0)+PARS(I0+2)*X(3)
75 X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
84 IF(IP.LE.4) X(3)=-X(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
89 IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
91 IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
93 IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
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
98 IF(MOD(IP,2).EQ.1) X(1)=-X(1)
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)
105 IF(ISH.NE.4) GO TO 20
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)
118 IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
120 XPT=DXS+(DX(1)*XT(1)+DX(2)*XT(2))/DXS
121 YPT=(DX(1)*XT(2)-DX(2)*XT(1))/DXS
123 IF(YPT.EQ.0.0.AND.XPT.EQ.0.0) GO TO 999
125 IF(P.GT.PI) P=P-PI*2.0
133 IF(PH-PL.GT.PI) GO TO 999
138 CL = MOD( ABS(CL),360.0 )
139 IF(SG.LE.0.0) CL=360.-CL
141 CH = MOD( ABS(CH),360.0 )
142 IF(SG.LE.0.0) CH=360.-CH
147 IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14)GO TO 80
155 C WHEN IN DOUBT SET TO FULL RANGE.
158 IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 50
162 ** approxime to a cylinder whit radius
163 ** equal to the ellipse major axis
165 IF(PARS(1).GT.RM) RM=PARS(1)
169 RM = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
173 IF(PARS(5).GT.RM) RM=PARS(5)
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
182 * This is the simple case, no rotation
183 * Compute the position of the limits on
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
194 * Rotated tubes might be more difficult
195 * Just leave it for later
198 IF(ISH.LE.6) GO TO 60
203 IF(IROT.EQ.0) GO TO 65
205 IF(Q(JROT+15).NE.0.0.AND.Q(JROT+15).NE.180.0) GO TO 999
209 IF(PHY.LT.PHX) PHY=PHY+360.0
211 IF(PHY-PHX.GT.180.0) ISPH=-1
214 IF(ISPH.EQ.1) GO TO 65
222 CL = MOD( ABS(CL),360.0 )
223 IF(SG.LE.0.0) CL=360.-CL
225 CH = MOD( ABS(CH),360.0 )
226 IF(SG.LE.0.0) CH=360.-CH
232 C DISPLACEMENT GREATER THAN MAXIMUM RADIUS SO
233 C ASSUME COMPLETE TUBE TO GENERATE 'WORST CASE'.
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
243 ELSEIF(ISH.GT.6.AND.ISH.NE.13.AND.ISH.NE.14) THEN
257 IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
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)
266 XPT=DXS-(COST*DZ+SINT*RM)/DXS
268 IF(XPT.LE.0.0) GO TO 999
269 YPT=(SIN2*RM+COS2*DZ)/DXS
272 P0=ATAN2(DX(2),DX(1))
277 CL = MOD( ABS(CL),360.0 )
278 IF(SG.LE.0.0) CL=360.-CL
280 CH = MOD( ABS(CH),360.0 )
281 IF(SG.LE.0.0) CH=360.-CH
286 IF(ISH.GT.9) GO TO 999
294 IF(IROT.NE.0.OR.DXS.GT.0.0) GO TO 90
296 C UNROTATED AND CENTERED.
302 CL = MOD( ABS(CL),360.0 )
303 IF(SG.LE.0.0) CL=360.-CL
305 CH = MOD( ABS(CH),360.0 )
306 IF(SG.LE.0.0) CH=360.-CH
312 C ROTATED OR NOT CENTERED.
314 IF(DXS.LT.PARS(2)) GO TO 999
315 P0=ATAN2(DX(2),DX(1))
321 CL = MOD( ABS(CL),360.0 )
322 IF(SG.LE.0.0) CL=360.-CL
324 CH = MOD( ABS(CH),360.0 )
325 IF(SG.LE.0.0) CH=360.-CH