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 GFLTHE(ISH,IROT,DX,PARS,CL,CH,IERR)
14 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 *
23 C. * ==>Called by : GFCLIM *
24 C. * Author A.McPherson ********* *
26 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)
33 C. ----------------------------------------------
37 DXS=DX(1)*DX(1)+DX(2)*DX(2)
38 IF(DXS.GT.0.0) DXS=SQRT(DXS)
40 IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 130
42 C CUBES, TRAPEZOIDS AND PARALLELEPIPEDS.
48 IF(DXS.LE.0.0) GO TO 999
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
58 C THIS IS A LOOP OVER THE 8 CORNERS.
59 C FIRST FIND THE LOCAL COORDINATES.
63 C General twisted trapezoid.
69 X(1)=PARS(I0)+PARS(I0+2)*X(3)
70 X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
79 IF(IP.LE.4) X(3)=-X(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
84 IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
86 IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
88 IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
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
93 IF(MOD(IP,2).EQ.1) X(1)=-X(1)
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)
100 IF(ISH.NE.4) GO TO 20
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)
113 IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
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
120 IF(R.GT.0.0) GO TO 30
127 IF(ABS(Z).LT.1.0E-6) GO TO 40
131 IF(T.EQ.0.0.AND.Z.LT.0.0) T=180.0
132 IF(T.LT.0.0) T=T+180.0
140 C THETA LIMITS SET FROM THE POINTS NOW DO THE EDGES.
144 C FIND THE END POINT NUMBERS FOR EACH EDGE.
150 IF(MOD(IL,2).EQ.1) IPP2=3
159 IF(MOD(IL,2).EQ.1) IPP2=7
169 C NOW GET THE POINTS AND ROTATE THEM.
173 C General twisted trapezoid.
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)
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)
200 IF(IPP1.LE.4) X1(3)=-X1(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
205 IF(ISH.EQ.4.AND.X1(3).GT.0.0) IP2=8
207 IF(MOD(IPP1+3,4).LT.2) X1(2)=-X1(2)
209 IF(ISH.NE.1.AND.ISH.NE.10.AND.X1(3).GT.0.0) IP1=2
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
214 IF(MOD(IPP1,2).EQ.1) X1(1)=-X1(1)
220 IF(IPP2.LE.4) X2(3)=-X2(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
225 IF(ISH.EQ.4.AND.X2(3).GT.0.0) IP2=8
227 IF(MOD(IPP2+3,4).LT.2) X2(2)=-X2(2)
229 IF(ISH.NE.1.AND.ISH.NE.10.AND.X2(3).GT.0.0) IP1=2
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
234 IF(MOD(IPP2,2).EQ.1) X2(1)=-X2(1)
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)
243 IF(ISH.NE.4) GO TO 100
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)
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)
260 IF(IROT.NE.0) CALL GINROT(X1,Q(JROT+1),XT1)
264 IF(IROT.NE.0) CALL GINROT(X2,Q(JROT+1),XT2)
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
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
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
283 IF(Z0.EQ.0.0.AND.ALZ.EQ.0.0) GO TO 110
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.
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)
293 IF(ABS(SNUM).GT.0.5*DS*ABS(SDEN)) GO TO 120
295 C TANGENT EXIST BETWEEN THE TWO ENDS.
305 IF(RT.LT.1.0E-10) GO TO 999
307 IF(R.GT.0.0) R=SQRT(R)
309 IF(ABS(Z0).LT.1.0E-06) GO TO 110
312 IF(T.EQ.0.0.AND.Z0.LT.0.0) T=180.0
313 IF(T.LT.0.0) T=T+180.0
320 C CHECK FOR THE POSSIBILITY OF STRADDLING T=0.0 &/OR 180.0.
323 IF(C.GT.0.0) GO TO 120
325 C CHECK IF SAME SIGN OF Z.
327 IF(Z0*DX(3).LT.0.0) GO TO 999
330 IF(Z0.LT.0.0) T=180.0
336 C DONE SET THE LIMITS.
347 IF(IROT.NE.0.OR.DXS.GT.1.0E-05) GO TO 180
349 C UNROTATED AND CENTERED.
351 IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13
352 + .AND.ISH.NE.14)GO TO 170
362 ** approxime to a cylinder whit radius
363 ** equal to the ellipse major axis
366 IF(PARS(1).GT.RMX) RMX=PARS(1)
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)
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
382 IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 140
386 IF(PARS(4).LT.RMN) RMN=PARS(4)
387 IF(PARS(5).GT.RMX) RMX=PARS(5)
390 IF(DZ.GT.ABS(DX(3))) GO TO 160
392 C ALL FORWARD OR ALL BACK.
396 IF(DX(3).GT.0.0) GO TO 150
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
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
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
427 IF(ISH.GT.9) GO TO 999
438 IF(ISH.EQ.11.OR.ISH.EQ.12) GOTO 999
441 IF(ISH.EQ.9) GO TO 200
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
454 IF(PARS(1).GT.RM) RM=PARS(1)
458 IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 190
462 IF(PARS(5).GT.RM) RM=PARS(5)
474 IF(RC.GT.0.0) RC=SQRT(RC)
475 IF(RM.GE.RC) GO TO 999
478 IF(ABS(DX(3)).GT.0.0) TC=ATAN(DXS/DX(3))*RADDEG
479 IF(TC.LT.0.0) TC=TC+180.0
481 DT=ASIN(RM/RC)*RADDEG
485 IF(CH.GT.180.0) CH=180.0