+++ /dev/null
-*
-* $Id$
-*
-* $Log$
-* Revision 1.1.1.1 1999/05/18 15:55:17 fca
-* AliRoot sources
-*
-* Revision 1.1.1.1 1995/10/24 10:20:57 cernlib
-* Geant
-*
-*
-#include "geant321/pilot.h"
-*CMZ : 3.21/02 29/03/94 15.41.31 by S.Giani
-*-- Author :
- SUBROUTINE GVDTHE(ISH,IROT,DX,PARS,CL,CH,IERR)
-C.
-C. ******************************************************************
-C. * *
-C. * ROUTINE TO COMPUTE THE THETA LIMITS FOR VOLUME OF SHAPE *
-C. * ISH ROTATED BY MATRIX IROT AND DISPLACED BY VECTOR DX. *
-C. * THE VOLUME HAS NPAR PARAMETERS IN THE ARRAY PARS. THE LOWER *
-C. * LIMIT IS RETURNED IN CL THE HIGHER IN CH. IF THE *
-C. * CALCULATION CANNOT BE MADE IERR IS SET TO 1 OTHERWISE IT *
-C. * IS SET TO 0. *
-C. * *
-C. * ==>Called by : GVDLIM *
-C. * Author S.Giani ********* *
-C. * *
-C. ******************************************************************
-C.
-#include "geant321/gcbank.inc"
-#include "geant321/gconsp.inc"
-#include "geant321/gcshno.inc"
- DIMENSION DX(3),PARS(100),X(3),XT(3),X1(3),X2(3),XT1(3),XT2(3)
-C.
-C. ----------------------------------------------
-C.
- IERR=1
-C
- DXS=DX(1)*DX(1)+DX(2)*DX(2)
- IF(DXS.GT.0.0) DXS=SQRT(DXS)
-C
- IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 130
-C
-C CUBES, TRAPEZOIDS AND PARALLELEPIPEDS.
-C
- IERR=0
- CL=0.0
- CH=180.0
-C
- IF(DXS.LE.0.0) GO TO 999
-C
- TH=90.
- IF(ABS(DX(3)).LT.1.0E-06)GO TO 5
- TH=ATAN(DXS/DX(3))*RADDEG
- IF(TH.LT.0.0) TH=TH+180.0
- 5 TL=TH
-C
- DO 50 IP=1,8
-C
-C THIS IS A LOOP OVER THE 8 CORNERS.
-C FIRST FIND THE LOCAL COORDINATES.
-C
- IF(ISH.EQ.28) THEN
- goto 999
-C General twisted trapezoid.
-C
-* IL=(IP+1)/2
-* I0=IL*4+11
-* IS=(IP-IL*2)*2+1
-* X(3)=PARS(1)*IS
-* X(1)=PARS(I0)+PARS(I0+2)*X(3)
-* X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
-* GO TO 20
-C
- ENDIF
-C
- IP3=ISH+2
- IF(ISH.EQ.10) IP3=3
- IF(ISH.EQ.4) IP3=1
- X(3)=PARS(IP3)
- IF(IP.LE.4) X(3)=-X(3)
- IP2=3
- IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
- IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
- IF(ISH.EQ.4) IP2=4
- IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
- X(2)=PARS(IP2)
- IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
- IP1=1
- IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
- IF(ISH.EQ.4) IP1=5
- IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
- IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
- X(1)=PARS(IP1)
- IF(MOD(IP,2).EQ.1) X(1)=-X(1)
-C
- IF(ISH.NE.10) GO TO 10
- X(1)=X(1)+X(2)*PARS(4)+X(3)*PARS(5)
- X(2)=X(2)+X(3)*PARS(6)
- 10 CONTINUE
-C
- IF(ISH.NE.4) GO TO 20
- IP4=7
- IF(X(3).GT.0.0) IP4=11
- X(1)=X(1)+X(2)*PARS(IP4)+X(3)*PARS(2)
- X(2)=X(2)+X(3)*PARS(3)
- 20 CONTINUE
-C
-C ROTATE.
-C
- JROT=LQ(JROTM-IROT)
- XT(1)=X(1)
- XT(2)=X(2)
- XT(3)=X(3)
- IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
-C
- Z=DX(3)+XT(3)
- R=(DX(1)+XT(1))**2+(DX(2)+XT(2))**2
- RT=R+(DX(3)+XT(3))**2
- IF(RT.LT.1.0E-10) GO TO 999
-C
- IF(R.GT.0.0) GO TO 30
- IF(Z.GT.0.0) T=0.0
- IF(Z.LT.0.0) T=180.0
- GO TO 40
-C
- 30 CONTINUE
- T=90.0
- IF(ABS(Z).LT.1.0E-6) GO TO 40
-C
- R=SQRT(R)
- T=ATAN(R/Z)*RADDEG
- IF(T.EQ.0.0.AND.Z.LT.0.0) T=180.0
- IF(T.LT.0.0) T=T+180.0
-C
- 40 CONTINUE
- IF(T.GT.TH) TH=T
- IF(T.LT.TL) TL=T
-C
- 50 CONTINUE
-C
-C THETA LIMITS SET FROM THE POINTS NOW DO THE EDGES.
-C
- DO 120 IL=1,12
-C
-C FIND THE END POINT NUMBERS FOR EACH EDGE.
-C
- IF(IL.GT.4) GO TO 60
- IPP1=1
- IF(IL.GT.2) IPP1=4
- IPP2=2
- IF(MOD(IL,2).EQ.1) IPP2=3
-C
- GO TO 80
- 60 CONTINUE
- IF(IL.LT.9) GO TO 70
-C
- IPP1=5
- IF(IL.GT.10) IPP1=8
- IPP2=6
- IF(MOD(IL,2).EQ.1) IPP2=7
-C
- GO TO 80
- 70 CONTINUE
-C
- IPP1=IL-4
- IPP2=IL
-C
- 80 CONTINUE
-C
-C NOW GET THE POINTS AND ROTATE THEM.
-C
- IF(ISH.EQ.28) THEN
-C
-C General twisted trapezoid.
-C
- ILL=IPP1
- IF(IPP1.EQ.3) ILL=4
- IF(IPP1.EQ.4) ILL=3
- I0=ILL*4+11
- X1(3)=PARS(1)
- IF(IPP1.LT.5) X1(3)=-X1(3)
- X1(1)=PARS(I0)+PARS(I0+2)*X1(3)
- X1(2)=PARS(I0+1)+PARS(I0+3)*X1(3)
- ILL=IPP2
- IF(IPP2.EQ.3) ILL=4
- IF(IPP2.EQ.4) ILL=3
- I0=ILL*4+11
- X2(3)=PARS(1)
- IF(IPP2.LT.5) X2(3)=-X2(3)
- X2(1)=PARS(I0)+PARS(I0+2)*X2(3)
- X2(2)=PARS(I0+1)+PARS(I0+3)*X2(3)
-C
- GO TO 100
-C
- ENDIF
-C
- IP3=ISH+2
- IF(ISH.EQ.10) IP3=3
- IF(ISH.EQ.4) IP3=1
- X1(3)=PARS(IP3)
- IF(IPP1.LE.4) X1(3)=-X1(3)
- IP2=3
- IF(ISH.GT.2.AND.X1(3).GT.0.0) IP2=4
- IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
- IF(ISH.EQ.4) IP2=4
- IF(ISH.EQ.4.AND.X1(3).GT.0.0) IP2=8
- X1(2)=PARS(IP2)
- IF(MOD(IPP1+3,4).LT.2) X1(2)=-X1(2)
- IP1=1
- IF(ISH.NE.1.AND.ISH.NE.10.AND.X1(3).GT.0.0) IP1=2
- IF(ISH.EQ.4) IP1=5
- IF(ISH.EQ.4.AND.X1(3).GT.0.0) IP1=IP1+4
- IF(ISH.EQ.4.AND.X1(2).GT.0.0) IP1=IP1+1
- X1(1)=PARS(IP1)
- IF(MOD(IPP1,2).EQ.1) X1(1)=-X1(1)
-C
- IP3=ISH+2
- IF(ISH.EQ.10) IP3=3
- IF(ISH.EQ.4) IP3=1
- X2(3)=PARS(IP3)
- IF(IPP2.LE.4) X2(3)=-X2(3)
- IP2=3
- IF(ISH.GT.2.AND.X2(3).GT.0.0) IP2=4
- IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
- IF(ISH.EQ.4) IP2=4
- IF(ISH.EQ.4.AND.X2(3).GT.0.0) IP2=8
- X2(2)=PARS(IP2)
- IF(MOD(IPP2+3,4).LT.2) X2(2)=-X2(2)
- IP1=1
- IF(ISH.NE.1.AND.ISH.NE.10.AND.X2(3).GT.0.0) IP1=2
- IF(ISH.EQ.4) IP1=5
- IF(ISH.EQ.4.AND.X2(3).GT.0.0) IP1=IP1+4
- IF(ISH.EQ.4.AND.X2(2).GT.0.0) IP1=IP1+1
- X2(1)=PARS(IP1)
- IF(MOD(IPP2,2).EQ.1) X2(1)=-X2(1)
-C
- IF(ISH.NE.10) GO TO 90
- X1(1)=X1(1)+X1(2)*PARS(4)+X1(3)*PARS(5)
- X1(2)=X1(2)+X1(3)*PARS(6)
- X2(1)=X2(1)+X2(2)*PARS(4)+X2(3)*PARS(5)
- X2(2)=X2(2)+X2(3)*PARS(6)
- 90 CONTINUE
-C
- IF(ISH.NE.4) GO TO 100
- IP4=7
- IF(X1(3).GT.0.0) IP4=11
- X1(1)=X1(1)+X1(2)*PARS(IP4)+X1(3)*PARS(2)
- X1(2)=X1(2)+X1(3)*PARS(3)
- IP4=7
- IF(X2(3).GT.0.0) IP4=11
- X2(1)=X2(1)+X2(2)*PARS(IP4)+X2(3)*PARS(2)
- X2(2)=X2(2)+X2(3)*PARS(3)
- 100 CONTINUE
-C
-C ROTATE.
-C
- JROT=LQ(JROTM-IROT)
- XT1(1)=X1(1)
- XT1(2)=X1(2)
- XT1(3)=X1(3)
- IF(IROT.NE.0) CALL GINROT(X1,Q(JROT+1),XT1)
- XT2(1)=X2(1)
- XT2(2)=X2(2)
- XT2(3)=X2(3)
- IF(IROT.NE.0) CALL GINROT(X2,Q(JROT+1),XT2)
-C
-C NOW WE HAVE THE END POINTS IN THE MASTER SYSTEM.
-C FIND THE TANGENT POINT TO THE SET OF CONES OF CONSTANT
-C THETA.
-C
- DS=(XT2(1)-XT1(1))**2+(XT2(2)-XT1(2))**2+(XT2(3)-XT1(3))**2
- IF(DS.LE.0.0) GO TO 120
-C
- DS=SQRT(DS)
-C
- X0=(XT2(1)+XT1(1))/2.0+DX(1)
- Y0=(XT2(2)+XT1(2))/2.0+DX(2)
- Z0=(XT2(3)+XT1(3))/2.0+DX(3)
- ALX=(XT2(1)-XT1(1))/DS
- ALY=(XT2(2)-XT1(2))/DS
- ALZ=(XT2(3)-XT1(3))/DS
-C
- T=90.0
- IF(Z0.EQ.0.0.AND.ALZ.EQ.0.0) GO TO 110
-C
- IF(ALX.EQ.0.0.AND.ALY.EQ.0.0) GO TO 120
-C THIS CHECKS WHETHER THE LINE IS PARALLEL TO THE
-C Z AXIS IN WHICH CASE THERE IS NO TANGENT AND
-C THE END POINTS DETERMINE THE THETA RANGE.
-C
- SNUM=(X0*Z0*ALX+Y0*Z0*ALY+X0*X0*ALZ+Y0*Y0*ALZ)
- SDEN=(Z0*ALX*ALX-X0*ALX*ALZ+Z0*ALY*ALY-Y0*ALY*ALZ)
-C
- IF(ABS(SNUM).GT.0.5*DS*ABS(SDEN)) GO TO 120
-C
-C TANGENT EXIST BETWEEN THE TWO ENDS.
-C
- S = SNUM/SDEN
- X0=X0+S*ALX
- Y0=Y0+S*ALY
- Z0=Z0+S*ALZ
-C
- R=X0*X0+Y0*Y0
- RT=R+Z0*Z0
-C
- IF(RT.LT.1.0E-10) GO TO 999
-C
- IF(R.GT.0.0) R=SQRT(R)
- T=90.0
- IF(ABS(Z0).LT.1.0E-06) GO TO 110
-C
- T=ATAN(R/Z0)*RADDEG
- IF(T.EQ.0.0.AND.Z0.LT.0.0) T=180.0
- IF(T.LT.0.0) T=T+180.0
-C
- 110 CONTINUE
-C
- IF(T.LT.TL) TL=T
- IF(T.GT.TH) TH=T
-C
-C CHECK FOR THE POSSIBILITY OF STRADDLING T=0.0 &/OR 180.0.
-C
- C=X0*DX(1)+Y0*DX(2)
- IF(C.GT.0.0) GO TO 120
-C
-C CHECK IF SAME SIGN OF Z.
-C
- IF(Z0*DX(3).LT.0.0) GO TO 999
-C
- T=0.0
- IF(Z0.LT.0.0) T=180.0
- IF(T.LT.TL) TL=T
- IF(T.GT.TH) TH=T
-C
- 120 CONTINUE
-C
-C DONE SET THE LIMITS.
-C
- CL=TL
- CH=TH
-C
- GO TO 999
-C
- 130 CONTINUE
-C
-C TUBES, CONES ETC.
-C
- IF(IROT.NE.0.OR.DXS.GT.1.0E-05) GO TO 180
-C
-C UNROTATED AND CENTERED.
-C
- IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13
- + .AND.ISH.NE.14)GO TO 170
-C
-C TUBES AND CONES.
-C
- IERR=0
- DZ=PARS(3)
- RMN=PARS(1)
- RMX=PARS(2)
- IF(ISH.EQ.13) THEN
-**
-** approxime to a cylinder whit radius
-** equal to the ellipse major axis
-**
- RMN=0.0
- IF(PARS(1).GT.RMX) RMX=PARS(1)
- GOTO 140
- ENDIF
- IF(ISH.EQ.14) THEN
-C not really sure of the function of these... keep RMN=PARS(1)
- RMX = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
- GO TO 140
- ENDIF
- IF(ISH.EQ.NSCTUB) THEN
- S1 = (1.0-PARS(8))*(1.0+PARS(8))
- IF( S1 .GT. 0.0) S1 = SQRT(S1)
- S2 = (1.0-PARS(11))*(1.0+PARS(11))
- IF( S2 .GT. 0.0) S2 = SQRT(S2)
- IF( S2 .GT. S1 ) S1 = S2
- DZ = DZ+RMX*S1
- ENDIF
- IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 140
- DZ=PARS(1)
- RMN=PARS(2)
- RMX=PARS(3)
- IF(PARS(4).LT.RMN) RMN=PARS(4)
- IF(PARS(5).GT.RMX) RMX=PARS(5)
- 140 CONTINUE
-C
- IF(DZ.GT.ABS(DX(3))) GO TO 160
-C
-C ALL FORWARD OR ALL BACK.
-C
- DSH=DX(3)-DZ
- DLN=DX(3)+DZ
- IF(DX(3).GT.0.0) GO TO 150
- DSS=DSH
- DSH=DLN
- DLN=DSS
- 150 CONTINUE
-C
- CL=90.0
- CH=90.0
- IF(DLN.NE.0.0) CL=ATAN(RMN/DLN)*RADDEG
- IF(DSH.NE.0.0) CH=ATAN(RMX/DSH)*RADDEG
- IF(DX(3).GT.0.0) GO TO 999
- CS=CL
- CL=CH
- CH=CS
- IF(CH.EQ.0.0) CH=180.0
- IF(CL.LT.0.0) CL=CL+180.0
- IF(CH.LT.0.0) CH=CH+180.0
-C
- GO TO 999
-C
- 160 CONTINUE
-C
- CL=90.0
- CH=90.0
- IF(DZ+DX(3).NE.0.0) CL=ATAN(RMN/(DZ+DX(3)))*RADDEG
- IF(-DZ+DX(3).NE.0.0) CH=ATAN(RMN/(-DZ+DX(3)))*RADDEG
- IF(CH.LE.0.0) CH=CH+180.0
-C
- GO TO 999
-C
- 170 CONTINUE
- IF(ISH.GT.9) GO TO 999
-C
-C SPHERE.
-C
- IERR=0
- CL=PARS(3)
- CH=PARS(4)
-C
- GO TO 999
- 180 CONTINUE
-C
- IF(ISH.EQ.11.OR.ISH.EQ.12) GOTO 999
-**
- RM=PARS(2)
- IF(ISH.EQ.9) GO TO 200
-C
- DZ=PARS(3)
- IF(ISH.EQ.NSCTUB) THEN
- S1 = (1.0-PARS(8))*(1.0+PARS(8))
- IF( S1 .GT. 0.0) S1 = SQRT(S1)
- S2 = (1.0-PARS(11))*(1.0+PARS(11))
- IF( S2 .GT. 0.0) S2 = SQRT(S2)
- IF( S2 .GT. S1 ) S1 = S2
- DZ = DZ+RM*S1
- ENDIF
-**
- IF(ISH.EQ.13) THEN
- IF(PARS(1).GT.RM) RM=PARS(1)
- GOTO 190
- ENDIF
-**
- IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 190
-C
- DZ=PARS(1)
- RM=PARS(3)
- IF(PARS(5).GT.RM) RM=PARS(5)
-C
- 190 CONTINUE
-C
- RM=SQRT(RM**2+DZ**2)
-C
- 200 CONTINUE
-C
- CL=0.0
- CH=180.0
- IERR=0
- RC=DXS**2+DX(3)**2
- IF(RC.GT.0.0) RC=SQRT(RC)
- IF(RM.GE.RC) GO TO 999
-C
- TC=90.0
- IF(ABS(DX(3)).GT.0.0) TC=ATAN(DXS/DX(3))*RADDEG
- IF(TC.LT.0.0) TC=TC+180.0
-C
- DT=ASIN(RM/RC)*RADDEG
- CL=TC-DT
- CH=TC+DT
- IF(CL.LT.0.0) CL=0.0
- IF(CH.GT.180.0) CH=180.0
-C
- 999 CONTINUE
- END