戻る

PROGRAM MAIN
C ******************************************************************** C
C C
C Error Analysis Software C
C C
C for 3D - Convectin & Diffusion Equation C
C ＝＝ C
C ( for FEMs & FDMs ) C
C C
C ( V.3.5 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CST/ PAI
COMMON /DEB/ IDEB(10)
COMMON /NUM/ MTS, NT, NRX, NBX, NKX, NRY, NBY, NKY, NRZ, NBZ, NKZ
COMMON /ABC/ RVX(10), RVY(10), RVZ(10), BVX(10), BVY(10), BVZ(10),
2 KMX(14), KMY(14), KMZ(14)
C
COMMON /MTV/ MT(50), TV(14)
C
COMMON /OPTC/ OPTE( 20,50), OPTT( 20,50), OPTR( 20,50),
1 OPTG( 20,50), OPET( 50), OPT ( 50),
2 OPR ( 50), OPGZ( 50)
C
C --------------------------------------------
C
OPEN ( 7, FILE='EA_CVDF3D.OUT' )
C
C --------------------------------------------
C
CALL DIPT
C
C -------------------------------------------------------------------
DO 100 I = 1, MTS
C
C -------------------------
C
CALL PCON2 (I)
C
C -------------------------
KT = 0
RCU = 1.D0
C
C ----- For GRT / GCRT / MGRT / EXRT / EXCRT ---
C
IF ( MT(I).EQ.4 .OR. MT(I).EQ.8 .OR. MT(I).EQ.12
1 .OR. MT(I).EQ.18 .OR. MT(I).EQ.22)
2 RCU = ( 4.D0* DSQRT( 6.D0)/9.D0)** ( 1.D0/ 3.D0)
C
NTT = NT
IF ( MT(I).GE.15) NTT = 1
C
C ---------------------------------------------------------
DO 200 J = 1, NTT
C
T = TV( J)
KT = KT + 1
IF ( MT(I).GE.15) GO TO 10
WRITE(6,2000) T
WRITE(7,2000) T
2000 FORMAT(/,3X,'** Time Scheme Par. = ',F6.3,' **'/)
C
10 KR = 0
NBXYZ = NBX* NBY* NBZ
NKXYZ = NKX* NKY* NKZ
C
TSGA = 0.D0
TSGC = 0.D0
AB3 = 0.D0
C
C ------------------------------------------------
DO 300 IR = 1, NRX
C
RX = RVX( IR)* ( RCU** 2)
RY = RVY( IR)* ( RCU** 2)
RZ = RVZ( IR)* ( RCU** 2)
C
SGA = 0.D0
SGC = 0.D0
AB2 = 0.D0
RMX = 0.D0
C
C -----------------
C
CALL INIT
C
C -----------------
C
C --------------------------------------
DO 400 IBX = 1, NBX
C
BX = BVX( IBX)* RCU

C --------------------------------------
DO 400 IBY = 1, NBY
C
BY = BVY( IBY)* RCU
IF ( IDEB(2).EQ.0) BY = BX
C
C -------------------------------------
DO 400 IBZ = 1, NBZ
C
BZ = BVZ(IBZ)* RCU
IF ( IDEB(2).EQ.0) BZ = BX
BSA = 0.D0
BSC = 0.D0
AB1= 0.D0
RMXB= 0.D0
C
C ----------------------------
DO 500 IKX = NKX, 1, -1
C
KX = KMX( IKX)
VAX = 2.D0* PAI/ DFLOAT( KX)* RCU
C
C ----------------------------
DO 500 IKY = NKY, 1, -1
C
KY = KMY( IKY)
IF ( IDEB(3).EQ.0) KY = KX
VAY = 2.D0* PAI/ DFLOAT( KY)* RCU
C
C ----------------------------
DO 500 IKZ = NKZ, 1, -1
C
KZ = KMZ( IKZ)
IF ( IDEB(3).EQ.0) KZ = KX
VAZ = 2.D0* PAI/ DFLOAT( KZ)* RCU
C
C ----- Numerical Methods ------------------------------------------
C
CALL SMTD ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, I, T)
C
C ------------------------------------------------------------------
C
CALL CMP1 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, A, C, AB0, BA, BC, TS, AB1,
2 BSA, BSC, RMXB, SGA, SGC, RMX)
C
C ------------------------------------------------------------------
IF ( IDEB(6).GT.0)
1 CALL OUT1 ( RX, RY, RZ, BX, BY, BZ, KX, KY, KZ, A,
2 BA, BC, TS, AB0)
C
C ------------------------------------------------------------------
C
500 CONTINUE
C ----------------------------
C
CALL CMP2 ( NKXYZ, AB1, AB2, BSA, BSC, TBS )
C
C ----------------------------------------------------
C
IF ( IDEB(7).GT.0)
1 CALL OUT2 ( RX, RY, RZ, BX, BY, BZ, RMXB, BSA, BSC,
2 TBS, AB1 )
C
400 CONTINUE
C ------------------------------------------------
C
CALL CMP3 ( NBXYZ, NKXYZ, AB2, AB3, SGA, SGC, TSGA, TSGC,TVR )
C
C --------------------------------------
C
IF ( MT(I).LE.14) CALL COP1 ( KT, KR, T, RX, RMX, TVR )
C
C --------------------------------------
C
IF ( IDEB(8).GT.0) CALL OUT3 ( RX, RY, RZ, RMX, SGA, SGC, TVR,
1 AB2)
C
C --------------------------------------
C
300 CONTINUE
C ------------------------------------------------
C
CALL CMP4 ( NRX, NBXYZ, NKXYZ, AB3, TSGA, TSGC, SUMT)
C
C ------------------------------------------------------------
C
IF ( IDEB(9).GT.0) CALL OUT4 ( TSGA, TSGC, SUMT, AB3 )
C
200 CONTINUE
C ----------------------------------------------------------
C
IF ( MT(I).LE.14) CALL COP2 ( KT, KR )
C
IF ( MT(I).LE.14 .AND. IDEB(10).GT.0) ! Opt.THEM
C
1 CALL OUT5 ( KR, OPET, OPR, OPT, OPGZ) ! Opt. THETA
C
100 CONTINUE
C --------------------------------------------------------------------
WRITE(6,2100)
WRITE(7,2100)
2100 FORMAT(/10X,'*** END ***')
C
CLOSE( 7)
C
STOP
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,*) ' '
WRITE(6,*) '*****************************************************'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Error Analysis Software for C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C 3-Dimensional Convection- Diffusion Equation C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( FEMs & FDMs ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( V.3.5 ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Copyright 2011 : Yasuhiro MATSUDA C'
WRITE(6,*) 'C C'
WRITE(6,*) '*****************************************************'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMP1 ( P1, P2, P3, P4, RX, RY, RZ, BX,
1 BY, BZ, VAX, VAY, VAZ, AA, CC, AB0,
2 BA, BC, TS, AB1, BSA, BSC, RMXB, SGA,
3 SGC, RMX )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CST/ PAI
C
AE = ( P1* P3 + P2* P4)/( P3* P3 + P4* P4)
BE = ( P2* P3 - P1* P4)/( P3* P3 + P4* P4)
AA = DSQRT( AE*AE + BE* BE)
C
IF ( DABS( BE).LT.0.00001D0.AND.AE.GE.0.D0) TA = 0.D0
TA = DATAN2( BE, AE )
C
CC = ( TA + VAX* BX + VAY* BY + VAZ* BZ)
IF ( DABS( CC).GT.PAI ) CC = 2.D0* PAI - DABS( CC)
CC = CC + 1.D0
C
VAL = VAX* BX + VAY* BY + VAZ* BZ
C
* IF ( VAL.LE.0.000001) WRITE(6,*) ' ** VAL=', VAL
* CC = - TA /( VAX* BX + VAY* BY + VAZ* BZ)
C
C -----------------------------------------------------------
C
AP = DEXP( - VAX* VAX* RX - VAY* VAY* RY - VAZ* VAZ* RZ)
AE = AP* DCOS( VAX* BX + VAY* BY + VAZ* BZ) - AE
BE = - AP* DSIN( VAX* BX + VAY* BY + VAZ* BZ) - BE
C
AB1 = AB1 + AE*AE + BE* BE
AB0 = DSQRT( AE*AE + BE* BE )
C
IF ( AA.GE.RMXB) RMXB = AA
BSA = BSA + ( AA - AP)* ( AA - AP)
BSC = BSC + ( CC - 1.D0)* ( CC - 1.D0)
C
** BSC = BSC + TA*TA* ( CC- 1.D0)* (CC- 1.D0)
C
BA = DSQRT(( AA - AP)* ( AA - AP))
BC = DSQRT(( CC - 1.D0)* ( CC - 1.D0))
TS = DSQRT( BA* BA + BC* BC )
C
IF ( AA.GE.RMX) RMX = AA
C
SGA = SGA + ( AA - AP )* ( AA - AP)
SGC = SGC + ( CC - 1.D0)* ( CC - 1.D0)
C
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMP2 ( NKXYZ, AB1, AB2, BSA, BSC, TBS )
C
IMPLICIT REAL*8(A-H,O-Z)
C
AB2 = AB2 + AB1
C
AB1 = DSQRT( AB1/ DFLOAT( NKXYZ))
BSA = DSQRT( BSA/ DFLOAT( NKXYZ))
BSC = DSQRT( BSC/ DFLOAT( NKXYZ))
C
TBS = DSQRT( BSA* BSA + BSC * BSC )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMP3 ( NBXYZ, NKXYZ, AB2, AB3, SGA, SGC, TSGA, TSGC,
1 TVR )
C
IMPLICIT REAL*8(A-H,O-Z)
C
AB3 = AB2 + AB3
AB2 = DSQRT( AB2/( DFLOAT( NBXYZ)* DFLOAT( NKXYZ)))
TSGA = TSGA + SGA
C
SGA = DSQRT( SGA/( DFLOAT( NBXYZ)* DFLOAT( NKXYZ)))
TSGC = TSGC + SGC
C
SGC = DSQRT( SGC/( DFLOAT( NBXYZ)* DFLOAT( NKXYZ)))
TVR = DSQRT( SGA* SGA + SGC* SGC )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMP4 ( NRX, NBXYZ, NKXYZ, AB3, TSGA, TSGC, SUMT)
C
IMPLICIT REAL*8(A-H,O-Z)
C
AB3 = DSQRT( AB3/ ( DFLOAT( NRX)* DFLOAT( NBXYZ)* DFLOAT( NKXYZ)))
C
TSGA = DSQRT( TSGA/
1 ( DFLOAT( NRX)* DFLOAT( NBXYZ)* DFLOAT( NKXYZ)))
C
TSGC = DSQRT( TSGC/
1 ( DFLOAT( NRX)* DFLOAT( NBXYZ)* DFLOAT( NKXYZ)))
C
SUMT = DSQRT( TSGA* TSGA + TSGC* TSGC )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE COP1 ( KT, KR, T, RX, RMX, TVR )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /OPTC/ OPTE( 20,50), OPTT( 20,50), OPTR( 20,50),
1 OPTG( 20,50), OPET(50), OPT (50),
2 OPR (50), OPGZ(50)
C
KR = KR + 1
C
OPTE( KT,KR) = TVR
OPTT( KT,KR) = T
OPTR( KT,KR) = RX
OPTG( KT,KR) = RMX
C
RETURN
END
C **********************************************************************
C
SUBROUTINE COP2 ( KT, KR )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /OPTC/ OPTE( 20,50), OPTT( 20,50), OPTR( 20,50),
1 OPTG( 20,50), OPET(50), OPT (50),
2 OPR (50), OPGZ(50)
C
C ------------------------------------------------
DO 100 J = 1, KR
C
L = 0
C
C --------------------------------------
DO 200 I = 1, KT
C
IF ( OPTG( I,J).LE.1.D0) L = 1
200 CONTINUE
C
C --------------------------------------
DO 300 I = 1, KT
C
IF ( L.EQ.0) GO TO 10
IF ( OPTG( I,J).LE.1.D0 .AND. OPTE( I,J).LT.OPET( J)) GO TO 20
GO TO 300
C
10 IF ( OPTG( I,J).LT.OPGZ( J)) GO TO 20
GO TO 300
C
20 CONTINUE
OPET( J) = OPTE( I,J)
OPR ( J) = OPTR( I,J)
OPT ( J) = OPTT( I,J)
OPGZ( J) = OPTG( I,J)
C
300 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CRNSN ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ------------------------------
F1( X) = 1.D0 - DCOS( X)
F3( X) = DSIN( X)
C ------------------------------
C
S1 = 2.D0
S2 = - 2.D0* ( RX* F1( VAX) + RY* F1( VAY) + RZ* F1( VAZ))
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
S6 = - ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ))
C
S7 = 2.D0
S8 = 2.D0* ( RX* F1( VAX) + RY* F1( VAY) + RZ* F1( VAZ))
S9 = 0.D0
S10= 0.D0
S11= 0.D0
S12= BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ)
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXCHH ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = 2.D0 + DCOS( X)
F3( X) = DSIN( X)
C ------------------------------
C
RXF = 0.5D0* (( BX + BY + BZ)* ( BX + BY + BZ) - 1.D0)/ 3.D0
RYF = 0.5D0* (( BX + BY + BZ)* ( BX + BY + BZ) - 1.D0)/ 3.D0
RZF = 0.5D0* (( BX + BY + BZ)* ( BX + BY + BZ) - 1.D0)/ 3.D0
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 + ( BX + BY + BZ)* ( BX + BY + BZ)
1 /( 2.D0* ( RX + RY + RZ)) - 1.D0/( 2.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
C
S1 = F2( VAX)* F2( VAY)* F2( VAZ)
S2 = - 6.D0* ( RXF* F1( VAX)* F2( VAY)* F2( VAZ)
1 + RYF* F2( VAX)* F1( VAY)* F2( VAZ)
2 + RZF* F2( VAX)* F2( VAY)* F1( VAZ) )
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
S6 = - 3.D0* G* ( BX* F3( VAX)* F2( VAY)* F2( VAZ)
1 + BY* F2( VAX)* F3( VAY)* F2( VAZ)
2 + BZ* F2( VAX)* F2( VAY)* F3( VAZ))
C
S7 = 27.D0
S8 = 0.D0
S9 = 0.D0
S10= 0.D0
S11= 0.D0
S12= 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXCRT ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C -----------------------------
F1( X) = DCOS( X)
F2( X) = DCOS( X/ SQ3)
F3( X) = DSIN( X)
F4( X) = DSIN( X/ SQ3)
C -----------------------------
C
SQ3 = DSQRT( 3.D0)
SQ2 = DSQRT( 2.D0)
FII = 10.D0* ( 11.D0 + 4.D0* SQ3)/ 219.D0
C
RXF = FII* ( 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)
1 - 0.4D0)/ 3.D0
RYF = FII* ( 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)
1 - 0.4D0)/ 3.D0
RZF = FII* ( 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)
1 - 0.4D0)/ 3.D0
C
G = 10.D0* SQ2* ( 3.D0* SQ3 - 10.D0)/ 73.D0
C
IF (( RX + RY + RZ).GE.0.001 ) THEN
F = FII* ( 1.D0 + ( BX + BY + BZ)* ( BX + BY + BZ)
1 /( 2.D0* ( RX + RY + RZ)) - 2.D0/( 5.D0* ( RX + RY+RZ)))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 10.D0* SQ2* ( 3.D0* SQ3 - 10.D0)/ 73.D0
END IF
C
S1 = 2.D0 + F1( VAX)* F2( VAY) + F1( VAY)* F2( VAZ)
1 + F1( VAZ)* F2( VAX)
C
E1 = 24.D0 - 5.D0* SQ3
E2 = 4.D0 - SQ3
E3 = 5.D0 - 6.D0* SQ3
C
S2 = ( 3.D0/ 4.D0)* ( - 3.D0* ( 11.D0 - 4.D0* SQ3)
1 * ( RX + RY + RZ)
2 + ( E1* RXF + E2* RYF + E3* RZF)* F1( VAX)* F2( VAY)
3 + ( E3* RXF + E1* RYF + E2* RZF)* F1( VAY)* F2( VAZ)
4 + ( E2* RXF + E3* RYF + E1* RZF)* F1( VAZ)* F2( VAX))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
D1 = 4.D0
D2 = 3.D0 + 2.D0* SQ3
C
S6 = ( DSQRT( 6.D0)/ 4.D0)* G
1 * ( D1* BX* F4( VAX)* F1( VAZ) + D2* BZ* F2( VAX)*F3( VAZ)
2 + D2* BY* F3( VAY)* F2( VAZ) + D1* BZ* F1( VAY)*F4( VAZ)
3 + D2* BX* F3( VAX)* F2( VAY) + D1* BY* F1( VAX)*F4( VAY))
C
S7 = 5.D0
S8 = 0.D0
S9 = 0.D0
C
S10= 0.D0
S11= 0.D0
S12= 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXCTH1 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C -------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C -------------------------------
C
RXF = 0.25D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0 - 1.D0/ 20.D0
RYF = 0.25D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0 - 1.D0/ 20.D0
RZF = 0.25D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0 - 1.D0/ 20.D0
C
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001D0) THEN
F = 0.5D0* ( 1.D0 + ( BX + BY + BZ)* ( BX + BY + BZ)/
1 ( 2.D0* ( RX + RY + RZ))) - 3.D0/( 20.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
C
S1 = 4.D0 + 2.D0* ( F2( VAX) + F2( VAY) + F2( VAZ))
S2 = - 40.D0* ( RXF* F1( VAX) + RYF* F1( VAY) + RZF* F1( VAZ))
S3 = 0.D0
C
S4 = 0.D0
S5 = 0.D0
S6 = - 10.D0* G* ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ))
C
S7 = 10.D0
S8 = 0.D0
S9 = 0.D0
C
S10= 0.D0
S11= 0.D0
S12= 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXCTH2 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C -----------------------------------
F1( X,Y) = DCOS( X)* DCOS( Y)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C -----------------------------------
C
RXF = 0.25D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0 - 9.D0/ 40.D0
RYF = 0.25D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0 - 9.D0/ 40.D0
RZF = 0.25D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0 - 9.D0/ 40.D0
C
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 0.5D0* ( 1.D0 + ( BX + BY + BZ)* ( BX + BY + BZ)
1 /( 2.D0* ( RX + RY + RZ))) - 27.D0/( 40.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
C
S1 = 10.D0 + ( F2( VAX) + F2( VAY) + F2( VAZ))
1 + 4.D0* ( F1( VAX,VAY) + F1( VAY,VAZ) + F1( VAZ,VAX))
C
S2 = - 10.D0* ( 3.D0* ( RXF + RYF + RZF)
1 - 2.D0* ( RXF* F2( VAX) + RYF* F2( VAY) + RZF* F2( VAZ))
C
2 - ( RXF* ( F1( VAX,VAY) - F1( VAY,VAZ) + F1( VAZ,VAX))
3 + RYF* ( F1( VAX,VAY) + F1( VAY,VAZ) - F1( VAZ,VAX))
4 + RZF* ( - F1( VAX,VAY) + F1( VAY,VAZ) + F1( VAZ,VAX))))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
S6 = - 5.D0* G* ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ)
1 + ( BX+BY)* F3( VAX+VAY) + ( BY+BZ)* F3( VAY+VAZ)
2 + ( BZ+BX)* F3( VAZ+VAX)
C
3 + ( -BX+BY)* F3( -VAX+VAY) + ( -BY+BZ)* F3( -VAY+VAZ)
4 + ( -BZ+BX)* F3( -VAZ+VAX))
C
S7 = 25.D0
S8 = 0.D0
S9 = 0.D0
C
S10 = 0.D0
S11 = 0.D0
S12 = 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXFU ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ, I )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /MTV/ MT(50), TV(14)
C
IF ( MOD( MT(I),2).NE.0) GO TO 10
C
RXF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
1 - 0.5D0* ( BX + BY + BZ)/ 3.D0
C
RYF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
1 - 0.5D0* ( BX + BY + BZ)/ 3.D0
C
RZF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
1 - 0.5D0* ( BX + BY + BZ)/ 3.D0
C
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 + (( BX + BY + BZ)* ( BX + BY + BZ))
1 /( 2.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D0
END IF
GO TO 20
C
10 RXF = RX
RYF = RY
RZF = RZ
G = 1.D00
C
20 CONTINUE
S1 = 1.D0
C
S2 = - 2.D0* ( RXF* ( 1.D0 - DCOS( VAX))
1 + RYF* ( 1.D0 - DCOS( VAY))
2 + RZF* ( 1.D0 - DCOS( VAZ)))
C
3 - G* ( BX* ( 1.D0 - DCOS( VAX))
4 + BY* ( 1.D0 - DCOS( VAY))
5 + BZ* ( 1.D0 - DCOS( VAZ)))
C
S3 = - G* ( BX* DSIN( VAX) + BY* DSIN( VAY) + BZ* DSIN( VAZ))
C
S4 = 1.D0
S5 = 0.D0
S6 = 0.D0
C
P1 = S1 + S2
P2 = S3
P3 = S4 + S5
P4 = S6
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXFU3 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ, I )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /PRA/ ALFA
COMMON /MTV/ MT(50), TV(14)
C
C ------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C ------------------------------
C
IF ( MOD( MT(I),2).NE.0) GO TO 10
RXF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RYF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RZF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 + ( BX + BY + BZ)* ( BX + BY + BZ)
1 /( 2.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
GO TO 20
C
10 RXF = RX
RYF = RY
RZF = RZ
C
G = 1.D00
C
20 CONTINUE
S1 = 1.D0
S2 = - 2.D0* ( RXF* F1( VAX) + RYF* F1( VAY) + RZF* F1( VAZ))
C
S3 = G* ALFA* ( BX* ( 4.D0* F2( VAX) - F2( 2.D0* VAX) - 3.D0)
1 + BY* ( 4.D0* F2( VAY) - F2( 2.D0* VAY) - 3.D0)
2 + BZ* ( 4.D0* F2( VAZ) - F2( 2.D0* VAZ) - 3.D0))/ 6.D0
C
S4 = 0.D0
S5 = 0.D0
C
S6 = G* ( BX* ( F3( 2.D0* VAX) - 8.D0* F3( VAX))
1 + BY* ( F3( 2.D0* VAY) - 8.D0* F3( VAY))
2 + BZ* ( F3( 2.D0* VAZ) - 8.D0* F3( VAZ)))/ 6.D0
C
S7 = 1.D0
C
S8 = 0.D0
S9 = 0.D0
S10 = 0.D0
S11 = 0.D0
S12 = 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXHH ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = 2.D0 + DCOS( X)
F3( X) = DSIN( X)
C ------------------------------
C
S1 = F2( VAX)* F2( VAY)* F2( VAZ)
C
S2 = - 6.D0* ( RX* F1( VAX)* F2( VAY)* F2( VAZ)
1 + RY* F2( VAX)* F1( VAY)* F2( VAZ)
2 + RZ* F2( VAX)* F2( VAY)* F1( VAZ))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
S6 = - 3.D0* ( BX* F3( VAX)* F2( VAY)* F2( VAZ)
1 + BY* F2( VAX)* F3( VAY)* F2( VAZ)
2 + BZ* F2( VAX)* F2( VAY)* F3( VAZ))
C
S7 = 27.D0
C
S8 = 0.D0
S9 = 0.D0
S10= 0.D0
S11= 0.D0
S12= 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXRT ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C --------------------------------
F1( X) = DCOS( X)
F2( X) = DCOS( X/ SQ3)
F3( X) = DSIN( X)
F4( X) = DSIN( X/ SQ3)
C -------------------------------
C
SQ3 = DSQRT( 3.D0)
C
S1 = 2.D0 + F1( VAX)* F2( VAY) + F1( VAY)* F2( VAZ)
1 + F1( VAZ)* F2( VAX)
C
E1 = 24.D0 - 5.D0* SQ3
E2 = 4.D0 - SQ3
E3 = 5.D0 - 6.D0* SQ3
C
S2 = ( 3.D0/ 4.D0)* ( - 3.D0* ( 11.D0 - 4.D0* SQ3 )
1 * ( RX + RY + RZ)
2 + ( E1* RX + E2* RY + E3* RZ)* F1( VAX)* F2( VAY)
3 + ( E3* RX + E1* RY + E2* RZ)* F1( VAY)* F2( VAZ)
4 + ( E2* RX + E3* RY + E1* RZ)* F1( VAZ)* F2( VAX))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
D1 = 4.D0
D2 = 3.D0 + 2.D0* SQ3
C
S6 = ( DSQRT( 6.D0)/ 4.D0)
1 * ( D1* BX* F4( VAX)* F1( VAZ) + D2* BZ* F2( VAX)* F3( VAZ)
2 + D2* BY* F3( VAY)* F2( VAZ) + D1* BZ* F1( VAY)* F4( VAZ)
3 + D2* BX* F3( VAX)* F2( VAY) + D1* BY* F1( VAX)* F4( VAY))
C
S7 = 5.D0
C
S8 = 0.D0
S9 = 0.D0
S10 = 0.D0
S11 = 0.D0
S12 = 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXTH1 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C -------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C -------------------------------
C
S1 = 4.D0 + 2.D0* ( F2( VAX) + F2( VAY) + F2( VAZ))
S2 = - 40.D0* ( RX* F1( VAX) + RY* F1( VAY) + RZ* F1( VAZ))
S3 = 0.D0
C
S4 = 0.D0
S5 = 0.D0
S6 = - 10.D0* ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ))
C
S7 = 10.D0
C
S8 = 0.D0
S9 = 0.D0
S10= 0.D0
S11= 0.D0
S12= 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXTH2 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----------------------------------
F1(X,Y) = DCOS( X)* DCOS( Y)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C ----------------------------------
C
S1 = 10.D0 + ( F2( VAX) + F2( VAY) + F2( VAZ))
1 + 4.D0* ( F1( VAX,VAY) + F1( VAY,VAZ) + F1( VAZ,VAX))
C
S2 = - 10.D0* ( 3.D0* ( RX + RY + RZ)
1 - 2.D0* ( RX* F2( VAX) + RY* F2( VAY) + RZ* F2( VAZ))
2 - ( RX* ( F1( VAX,VAY) - F1( VAY,VAZ) + F1( VAZ,VAX))
3 + RY* ( F1( VAX,VAY) + F1( VAY,VAZ) - F1( VAZ,VAX))
4 + RZ* ( - F1( VAX,VAY) + F1( VAY,VAZ) + F1( VAZ,VAX))))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
S6 = - 5.D0* ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ)
1 + ( BX+BY)* F3( VAX+VAY) + ( BY+BZ)* F3( VAY+VAZ)
2 + ( BZ+BX)* F3( VAZ+VAX) + (-BX+BY)* F3(-VAX+VAY)
3 + (-BY+BZ)* F3(-VAY+VAZ) + (-BZ+BX)* F3(-VAZ+VAX))
C
S7 = 25.D0
C
S8 = 0.D0
S9 = 0.D0
S10 = 0.D0
S11 = 0.D0
S12 = 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FTCS ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ, VAX,
1 VAY, VAZ, I )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /MTV/ MT(50), TV(14)
C
IF ( MOD( MT(I),2).NE.0) GO TO 10
C
RXF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RYF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RZF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
C
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 + (( BX + BY + BZ)* ( BX + BY + BZ))
1 /( 2.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
GO TO 20
C
10 RXF = RX
RYF = RY
RZF = RZ
G = 1.D00
C
20 CONTINUE
S1 = 1.D0
S2 = - 2.D0* ( RXF* ( 1.D0 - DCOS( VAX))
1 + RYF* ( 1.D0 - DCOS( VAY))
2 + RZF* ( 1.D0 - DCOS( VAZ)))
C
S3 = - G* ( BX* DSIN( VAX) + BY* DSIN( VAY) + BZ* DSIN( VAZ))
C
S4 = 1.D0
S5 = 0.D0
S6 = 0.D0
C
P1 = S1 + S2
P2 = S3
P3 = S4 + S5
P4 = S6
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FULLY ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ, I )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /MTV/ MT(50), TV(14)
C
IF ( MOD( MT(I),2).NE.0) GO TO 10
C
RXF = - 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RYF = - 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RZF = - 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 - (( BX + BY + BZ)* ( BX + BY + BZ))
1 /( 2.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
GO TO 20
C
10 RXF = RX
RYF = RY
RZF = RZ
G = 1.D00
C
20 CONTINUE
S1 = 1.D0
S2 = 0.D0
S3 = 0.D0
S4 = 1.D0
C
S5 = 2.D0* ( RXF* ( 1.D0 - DCOS( VAX))
1 + RYF* ( 1.D0 - DCOS( VAY))
2 + RZF* ( 1.D0 - DCOS( VAZ)))
C
S6 = G* ( BX* DSIN( VAX) + BY* DSIN( VAY) + BZ* DSIN( VAZ))
C
P1 = S1 + S2
P2 = S3
P3 = S4 + S5
P4 = S6
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GCHH ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = 2.D0 + DCOS( X)
F3( X) = DSIN( X)
C ------------------------------
C
RXF = 0.5D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
C
RYF = 0.5D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
C
RZF = 0.5D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
C
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 + ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)* ( BX + BY + BZ)
1 /( 2.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
C
S1 = F2( VAX)* F2( VAY)* F2( VAZ)
C
S2 = - 6.D0* ( 1.D0 - T)* ( RXF* F1( VAX)* F2( VAY)* F2( VAZ)
1 + RYF* F2( VAX)* F1( VAY)* F2( VAZ)
2 + RZF* F2( VAX)* F2( VAY)* F1( VAZ))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
S6 = - 3.D0* ( 1.D0 - T)* G* ( BX* F3( VAX)* F2( VAY)* F2( VAZ)
1 + BY* F2( VAX)* F3( VAY)* F2( VAZ)
2 + BZ* F2( VAX)* F2( VAY)* F3( VAZ))
C
S7 = S1
C
S8 = 6.D0* T* ( RXF* F1( VAX)* F2( VAY)* F2( VAZ)
1 + RYF* F2( VAX)* F1( VAY)* F2( VAZ)
2 + RZF* F2( VAX)* F2( VAY)* F1( VAZ))
C
S9 = 0.D0
S10 = 0.D0
S11 = 0.D0
C
S12 = 3.D0* T* G* ( BX* F3( VAX)* F2( VAY)* F2( VAZ)
1 + BY* F2( VAX)* F3( VAY)* F2( VAZ)
2 + BZ* F2( VAX)* F2( VAY)* F3( VAZ))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GCRT ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ------------------------------
F1( X) = DCOS( X)
F2( X) = DCOS( X/ SQ3)
F3( X) = DSIN( X)
F4( X) = DSIN( X/ SQ3)
C ------------------------------
C
SQ2 = DSQRT( 2.D0)
SQ3 = DSQRT( 3.D0)
FII = 10.D0* (11.D0 + 4.D0* SQ3)/ 219.D0
C
RXF = FII* 0.5D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
RYF = FII* 0.5D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
RZF = FII* 0.5D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
C
G = 10.D0* SQ2* ( 3.D0* SQ3 - 10.D0)/ 73.D0
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = FII* ( 1.D0 + ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/( 2.D0* ( RX + RY + RZ)))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = SQ2* ( 3.D0* SQ3 - 10.D0)/ 73.D0
END IF
C
S1 = 2.D0 + F1( VAX)* F2( VAY) + F1( VAY)* F2( VAZ)
1 + F1( VAZ)* F2( VAX)
C
E1 = 24.D0 - 5.D0* SQ3
E2 = 4.D0 - SQ3
E3 = 5.D0 - 6.D0* SQ3
C
S2 = ( 3.D0/ 4.D0)* ( 1.D0 - T)
1 * ( - 3.D0* ( 11.D0 - 4.D0* SQ3 )* ( RXF + RYF + RZF)
2 + ( E1* RXF + E2* RYF + E3* RZF)* F1( VAX)* F2( VAY)
3 + ( E3* RXF + E1* RYF + E2* RZF)* F1( VAY)* F2( VAZ)
4 + ( E2* RXF + E3* RYF + E1* RZF)* F1( VAZ)* F2( VAX))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
D1 = 4.D0
D2 = 3.D0 + 2.D0* SQ3
C
S6 = ( DSQRT( 6.D0)/ 4.D0)* ( 1.D0 - T)* G
1 * ( D1* BX* F4( VAX)* F1( VAZ) + D2* BZ* F2( VAX)* F3( VAZ)
2 + D2* BY* F3( VAY)* F2( VAZ) + D1* BZ* F1( VAY)* F4( VAZ)
3 + D2* BY* F3( VAX)* F2( VAY) + D1* BY* F1( VAX)* F4( VAY))
C
S7 = S1
C
S8 = ( - 3.D0/ 4.D0)* T* ( - 3.D0* ( 11.D0 - 4.D0* SQ3 )
1 * ( RXF + RYF + RZF)
2 + ( E1* RXF + E2* RYF + E3* RZF)* F1( VAX)* F2( VAY)
3 + ( E3* RXF + E2* RYF + E2* RZF)* F1( VAY)* F2( VAZ)
4 + ( E2* RXF + E3* RYF + E1* RZF)* F1( VAZ)* F2( VAX))
C
S9 = 0.D0
S10= 0.D0
S11= 0.D0
C
S12= ( DSQRT( 6.D0)/ 4.D0)* ( - T)* G
1 * ( D1* BX* F4( VAX)* F1( VAZ) + D2* BZ* F2( VAX)* F3( VAZ)
2 + D2* BY* F3( VAY)* F2( VAZ) + D1* BZ* F1( VAY)* F4( VAZ)
3 + D2* BX* F3( VAX)* F2( VAY) + D1* BY* F1( VAX)* F4( VAY))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GCTH1 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C --------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C --------------------------------
C
RXF = 0.25D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
C
RYF = 0.25D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
C
RZF = 0.25D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
C
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 0.5D0* ( 1.D0 + (1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/( 2.D0* ( RX + RY + RZ)))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
C
S1 = 4.D0 + 2.D0* F2( VAX) + 2.D0* F2( VAY) + 2.D0* F2( VAZ)
S2 = - 40.D0* ( 1.D0 - T)* ( RXF* F1( VAX)
1 + RYF* F1( VAY) + RZF* F1( VAZ))
S3 = 0.D0
C
S4 = 0.D0
S5 = 0.D0
S6 = - 10.D0* ( 1.D0 - T)* G*
1 ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ))
C
S7 = S1
S8 = 40.D0* T* ( RXF* F1( VAX)
1 + RYF* F1( VAY) + RZF* F1( VAZ))
S9 = 0.D0
C
S10 = 0.D0
S11 = 0.D0
S12 = 10.D0* T* G* ( BX* F3( VAX)
1 + BY* F3( VAY) + BZ* F3( VAZ))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GCTH2 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----------------------------------
F1(X,Y) = DCOS( X)* DCOS( Y)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C ----------------------------------
C
RXF = 0.625D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
C
RYF = 0.625D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
C
RZF = 0.625D0* ( 1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/ 3.D0
C
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.25D0* ( 1.D0 + (1.D0 - 2.D0* T)* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/( 2.D0* ( RX + RY + RZ)))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
C
S1 = 10.D0 + ( F2( VAX) + F2( VAY) + F2( VAZ))
1 + 4.D0* ( F1( VAX,VAY) + F1( VAY,VAZ) + F1( VAZ,VAX))
C
S2 = - 10.D0* ( 1.D0 - T)* ( 3.D0* ( RXF + RYF + RZF)
1 - 2.D0* ( RXF* F2( VAX) + RYF* F2( VAY) + RZF* F2( VAZ))
2 - ( RXF* F1( VAX,VAY) - RXF* F1( VAY,VAZ) + RXF* F1( VAZ,VAX))
C
3 - ( RYF* F1( VAX,VAY) + RYF* F1( VAY,VAZ) - RYF* F1( VAZ,VAX))
4 - (-RZF* F1( VAX,VAY) + RZF* F1( VAY,VAZ) + RZF* F1( VAZ,VAX))
5 )
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
S6 = - 5.D0* ( 1.D0 - T)* G* ( BX* F3( VAX) + BY* F3( VAY)
1 + BZ* F3( VAZ) + ( BX+BY)* F3( VAX+VAY)
2 + ( BY+BZ)* F3( VAY+VAZ) + ( BZ+BX)* F3( VAZ+VAX)
C
3 + ( -BX+BY)* F3(-VAX+VAY) + (-BY+BZ)* F3(-VAY+VAZ)
4 + ( -BZ+BX)* F3(-VAZ+VAX))
C
S7 = S1
S8 = 10.D0* T* ( 3.D0* ( RXF + RYF + RZF)
1 - 2.D0* ( RXF* F2( VAX) + RYF* F2( VAY) + RZF* F2( VAZ))
C
2 - ( RXF* F1( VAX,VAY) - RXF* F1( VAY,VAZ) + RXF* F1( VAZ,VAX))
3 - ( RYF* F1( VAX,VAY) + RYF* F1( VAY,VAZ) - RYF* F1( VAZ,VAX))
4 - (-RZF* F1( VAX,VAY) + RZF* F1( VAY,VAZ) + RZF* F1( VAZ,VAX))
5 )
C
S9 = 0.D0
S10 = 0.D0
S11 = 0.D0
C
S12 = 5.D0* T* G* ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ)
1 + ( BX+BY)* F3( VAX+VAY) + ( BY+BZ)* F3( VAY+VAZ)
2 + ( BZ+BX)* F3( VAZ+VAX) + (-BX+BY)* F3(-VAX+VAY)
3 + ( -BY+BZ)* F3(-VAY+VAZ) + (-BZ+BX)* F3(-VAZ+VAX))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GHH ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C -------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = 2.D0 + DCOS( X)
F3( X) = DSIN( X)
C -------------------------------
C
S1 = F2( VAX)* F2( VAY)* F2( VAZ)
S2 = - 6.D0* ( 1.D0 - T)
1 * ( RX* F1( VAX)* F2( VAY)* F2( VAZ)
2 + RY* F2( VAX)* F1( VAY)* F2( VAZ)
3 + RZ* F2( VAX)* F2( VAY)* F1( VAZ))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
S6 = - 3.D0* ( 1.D0 - T)* ( BX* F3( VAX)* F2( VAY)* F2( VAZ)
1 + BY* F2( VAX)* F3( VAY)* F2( VAZ)
2 + BZ* F2( VAX)* F2( VAY)* F3( VAZ))
C
S7 = S1
S8 = 6.D0* T* ( RX* F1( VAX)* F2( VAY)* F2( VAZ)
1 + RY* F2( VAX)* F1( VAY)* F2( VAZ)
2 + RZ* F2( VAX)* F2( VAY)* F1( VAZ))
C
S9 = 0.D0
S10 = 0.D0
S11 = 0.D0
C
S12 = 3.D0* T* ( BX* F3( VAX)* F2( VAY)* F2( VAZ)
1 + BY* F2( VAX)* F3( VAY)* F2( VAZ)
2 + BZ* F2( VAX)* F2( VAY)* F3( VAZ) )
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GRT ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----------------------------
F1( X) = DCOS( X)
F2( X) = DCOS( X/ SQ3)
F3( X) = DSIN( X)
F4( X) = DSIN( X/ SQ3)
C ----------------------------
C
SQ3 = DSQRT( 3.D0)
C
S1 = 2.D0 + F1( VAX)* F2( VAY) + F1( VAY)* F2( VAZ)
1 + F1( VAZ)* F2( VAX)
C
E1 = 24.D0 - 5.D0* SQ3
E2 = 4.D0 - SQ3
E3 = 5.D0 - 6.D0* SQ3
C
S2 = ( 3.D0/ 4.D0)* ( 1.D0 - T)* ( - 3.D0* ( 11.D0 - 4.D0* SQ3 )
1 * ( RX + RY + RZ)
2 + ( E1* RX + E2* RY + E3* RZ) * F1( VAX)* F2( VAY)
3 + ( E3* RX + E1* RY + E2* RZ) * F1( VAY)* F2( VAZ)
4 + ( E2* RX + E3* RY + E1* RZ) * F1( VAZ)* F2( VAX))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
D1 = 4.D0
D2 = 3.D0 + 2.D0* SQ3
C
S6 = ( DSQRT( 6.D0)/ 4.D0)* ( 1.D0 - T)
1 * ( D1* BX* F4( VAX)* F1( VAZ) + D2* BZ* F2( VAX)* F3( VAZ)
2 + D2* BY* F3( VAY)* F2( VAZ) + D1* BZ* F1( VAY)* F4( VAZ)
3 + D2* BY* F3( VAX)* F2( VAY) + D1* BY* F1( VAX)* F4( VAY))
C
S7 = S1
S8 = ( - 3.D0/ 4.D0)* T* ( - 3.D0* ( 11.D0 - 4.D0* SQ3 )
1 * ( RX + RY + RZ)
C
2 + ( E1* RX + E2* RY + E3* RZ)* F1( VAX)* F2( VAY)
3 + ( E3* RX + E2* RY + E2* RZ)* F1( VAY)* F2( VAZ)
4 + ( E2* RX + E3* RY + E1* RZ)* F1( VAZ)* F2( VAX))
C
S9 = 0.D0
S10 = 0.D0
S11 = 0.D0
C
S12 = ( DSQRT( 6.D0)/ 4.D0)* ( - T)
1 * ( D1* BX* F4( VAX)* F1( VAZ) + D2* BZ* F2( VAX)* F3( VAZ)
2 + D2* BY* F3( VAY)* F2( VAZ) + D1* BZ* F1( VAY)* F4( VAZ)
3 + D2* BX* F3( VAX)* F2( VAY) + D1* BY* F1( VAX)* F4( VAY))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GTH1 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C -------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C -------------------------------
C
S1 = 4.D0 + 2.D0* F2( VAX) + 2.D0* F2( VAY) + 2.D0* F2( VAZ)
S2 = - 40.D0* ( 1.D0 - T)* ( RX* F1( VAX)
1 + RY* F1( VAY) + RZ* F1( VAZ))
S3 = 0.D0
C
S4 = 0.D0
S5 = 0.D0
S6 = - 10.D0* ( 1.D0 - T)* ( BX* F3( VAX) + BY* F3( VAY)
1 + BZ* F3( VAZ))
C
S7 = S1
S8 = 40.D0* T* ( RX* F1( VAX) + RY* F1( VAY) + RZ* F1( VAZ))
S9 = 0.D0
C
S10 = 0.D0
S11 = 0.D0
S12 = 10.D0* T* ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GTH2 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C -----------------------------------
F1( X,Y) = DCOS( X)* DCOS( Y)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C -----------------------------------
C
S1 = 10.D0 + ( F2( VAX) + F2( VAY) + F2( VAZ))
1 + 4.D0* ( F1( VAX,VAY) + F1( VAY,VAZ) + F1( VAZ,VAX))
C
S2 = - 10.D0* ( 1.D0 - T)* ( 3.D0* ( RX + RY + RZ)
1 - 2.D0* ( RX* F2( VAX) + RY* F2( VAY) + RZ* F2( VAZ))
C
2 -( RX* F1( VAX,VAY) - RX* F1( VAY,VAZ) + RX* F1( VAZ,VAX))
3 -( RY* F1( VAX,VAY) + RY* F1( VAY,VAZ) - RY* F1( VAZ,VAX))
4 -(-RZ* F1( VAX,VAY) + RZ* F1( VAY,VAZ) + RZ* F1( VAZ,VAX)))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
S6 = - 5.D0* ( 1.D0 - T)
1 * ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ)
C
2 + ( BX+BY)* F3( VAX+VAY) + ( BY+BZ)* F3( VAY+VAZ)
3 + ( BZ+BX)* F3( VAZ+VAX) + (-BX+BY)* F3(-VAX+VAY)
4 + (-BY+BZ)* F3(-VAY+VAZ) + (-BZ+BX)* F3(-VAZ+VAX))
C
S7 = S1
C
S8 = 10.D0* T* ( 3.D0* ( RX + RY + RZ)
1 - 2.D0* ( RX* F2( VAX) + RY* F2( VAY) + RZ* F2( VAZ))
C
2 - ( RX* F1( VAX,VAY) - RX* F1( VAY,VAZ) + RX* F1( VAZ,VAX))
3 - ( RY* F1( VAX,VAY) + RY* F1( VAY,VAZ) - RY* F1( VAZ,VAX))
4 - (-RZ* F1( VAX,VAY) + RZ* F1( VAY,VAZ) + RZ* F1( VAZ,VAX)))
C
S9 = 0.D0
S10 = 0.D0
S11 = 0.D0
C
S12= 5.D0* T* ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ)
1 + ( BX+BY)* F3( VAX+VAY) + ( BY+BZ)* F3( VAY+VAZ)
2 + ( BZ+BX)* F3( VAZ+VAX) + (-BX+BY)* F3(-VAX+VAY)
3 + ( -BY+BZ)* F3(-VAY+VAZ) + (-BZ+BX)* F3(-VAZ+VAX))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE IMFU ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, I )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /MTV/ MT(50), TV(14)
C
C ------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C ------------------------------
C
IF ( MOD( MT(I),2).NE.0) GO TO 10
C
RXF = - 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
1 - 0.5D0* ( BX + BY + BZ)/ 3.D0
C
RYF = - 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
1 - 0.5D0* ( BX + BY + BZ)/ 3.D0
C
RZF = - 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
1 - 0.5D0* ( BX + BY + BZ)/ 3.D0
C
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 - 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)
1 /( RX + RY + RZ)
2 - 0.5D0* ( BX + BY + BZ)/( RX + RY + RZ)
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
GO TO 20
C
10 RXF = RX
RYF = RY
RZF = RZ
G = 1.D00
C
20 CONTINUE
S1 = 1.D0
S2 = 0.D0
S3 = 0.D0
S4 = 0.D0
C
S5 = 0.D0
S6 = 0.D0
S7 = 1.D0
C
S8 = 2.D0* ( RXF* F1( VAX) + RYF* F1( VAY) + RZF* F1( VAZ))
S9 = G* ( BX* F1( VAX) + BY* F1( VAY) + BZ* F1( VAZ))
C
S10= 0.D0
S11= 0.D0
S12= G* ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INIT
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /OPTC/ OPTE( 20,50), OPTT( 20,50), OPTR( 20,50),
1 OPTG( 20,50), OPET(50), OPT (50),
2 OPR (50), OPGZ(50)
C
DO 100 I = 1, 50
C
OPET(I) = 100.D0
OPGZ(I) = 100.D0
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE KOSLA ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, I )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /MTV/ MT(50), TV(14)
C
C ------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C ------------------------------
C
IF ( MOD( MT(I),2).NE.0) GO TO 10
C
RXF = - 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RYF = - 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RZF = - 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 - ( BX + BY + BZ)* ( BX + BY + BZ)
1 /( 2.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
GO TO 20
C
10 RXF = RX
RYF = RY
RZF = RZ
G = 1.D00
C
20 CONTINUE
S1 = 1.D0
S2 = 0.D0
S3 = G* ( BX* F1( VAX) + BY* F1( VAY) + BZ* F1( VAZ))
C
S4 = 0.D0
S5 = 0.D0
S6 = 0.D0
C
S7 = 1.D0
S8 = 2.D0* ( RXF* F1( VAX) + RYF* F1( VAY) + RZF* F1( VAZ))
S9 = G* ( BX * F1( VAX) + BY * F1( VAY) + BZ * F1( VAZ))
C
S10 = 0.D0
S11 = 0.D0
S12 = G* ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE LXWDF ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, I )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /MTV/ MT(50), TV(14)
C
C ------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C ------------------------------
C
IF ( MOD( MT(I),2).NE.0) GO TO 10
C
RXF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
1 - 0.5D0* ( BX* BX + BY* BY + BZ* BZ)/ 3.D0
C
RYF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
1 - 0.5D0* ( BX* BX + BY* BY + BZ* BZ)/ 3.D0
C
RZF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
1 - 0.5D0* ( BX* BX + BY* BY + BZ* BZ)/ 3.D0
C
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 + ( BX + BY + BZ)* ( BX + BY + BZ)
1 /( 2.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
GO TO 20
C
10 RXF = RX
RYF = RY
RZF = RZ
G = 1.D00
C
20 CONTINUE
S1 = 1.D0
S2 = - 2.D0* ( RXF* F1( VAX) + RYF* F1( VAY) + RZF* F1( VAZ))
1 - G* ( BX* BX* F1( VAX) + BY* BY* F1( VAY)
2 + BZ* BZ* F1( VAZ))
S3 = 0.D0
C
S4 = 0.D0
S5 = 0.D0
S6 = - G* ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ) )
C
S7 = 1.D0
S8 = 0.D0
S9 = 0.D0
C
S10= 0.D0
S11= 0.D0
S12= 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MGHH ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = 2.D0 + DCOS( X)
F3( X) = DSIN( X)
C ------------------------------
C
RXF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RYF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RZF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
G = 1.D0
C
IF (( RX + RY + RZ).GE.1.D-10) THEN
F = 1.D0 + 0.5D0* ( BX + BY + BZ)
1 * ( BX + BY + BZ)/( RX + RY + RZ)
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D0
END IF
C
S1 = F2( VAX)* F2( VAY)* F2( VAZ)
C
S2 = - 6.D0* ( 1.D0 - T)
1 * ( RXF* F1( VAX)* F2( VAY)* F2( VAZ)
2 + RYF* F2( VAX)* F1( VAY)* F2( VAZ)
3 + RZF* F2( VAX)* F2( VAY)* F1( VAZ))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
S6 = - 3.D0* G* ( BX* F3( VAX)* F2( VAY)* F2( VAZ)
1 + BY* F2( VAX)* F3( VAY)* F2( VAZ)
2 + BZ* F2( VAX)* F2( VAY)* F3( VAZ))
C
S7 = S1
S8 = 6.D0* T* ( RXF* F1( VAX)* F2( VAY)* F2( VAZ)
1 + RYF* F2( VAX)* F1( VAY)* F2( VAZ)
2 + RZF* F2( VAX)* F2( VAY)* F1( VAZ))
C
S9 = 0.D0
S10 = 0.D0
S11 = 0.D0
S12 = 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MGHHQ1 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ,T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
RXF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RYF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RZF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
C
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 + (( BX + BY + BZ)* ( BX + BY + BZ))
1 /( 2.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
C
S1 = 224.D0
1 + 34.D0* DCOS( 2.D0* VAX + 2.D0* VAY + 2.D0* VAZ)
2 + 34.D0* DCOS( 2.D0* VAX + 2.D0* VAY - 2.D0* VAZ)
3 + 34.D0* DCOS( 2.D0* VAX - 2.D0* VAY + 2.D0* VAZ)
4 + 34.D0* DCOS( 2.D0* VAX - 2.D0* VAY - 2.D0* VAZ)
C
5 - 36.D0* DCOS( 2.D0* VAX + 2.D0* VAY + VAZ)
6 - 36.D0* DCOS( 2.D0* VAX + 2.D0* VAY - VAZ)
7 - 36.D0* DCOS( 2.D0* VAX - 2.D0* VAY + VAZ)
8 - 36.D0* DCOS( 2.D0* VAX - 2.D0* VAY - VAZ)
C
9 - 36.D0* DCOS( 2.D0* VAX + VAY + 2.D0* VAZ)
A - 36.D0* DCOS( 2.D0* VAX + VAY - 2.D0* VAZ)
B - 36.D0* DCOS( 2.D0* VAX - VAY + 2.D0* VAZ)
C - 36.D0* DCOS( 2.D0* VAX - VAY - 2.D0* VAZ)
C
D - 36.D0* DCOS( VAX + 2.D0* VAY + 2.D0* VAZ)
E - 36.D0* DCOS( VAX + 2.D0* VAY - 2.D0* VAZ)
F - 36.D0* DCOS( VAX - 2.D0* VAY + 2.D0* VAZ)
G - 36.D0* DCOS( VAX - 2.D0* VAY - 2.D0* VAZ)
C
H + 80.D0* DCOS( 2.D0* VAX + 2.D0* VAY )
I + 80.D0* DCOS( 2.D0* VAX - 2.D0* VAY )
J + 80.D0* DCOS( 2.D0* VAX + 2.D0* VAZ )
C
S1 = S1
1 + 80.D0* DCOS( 2.D0* VAX - 2.D0* VAZ)
2 + 80.D0* DCOS( 2.D0* VAY + 2.D0* VAZ)
3 + 80.D0* DCOS( 2.D0* VAY - 2.D0* VAZ)
C
4 - 104.D0* DCOS( 2.D0* VAX + VAY)
5 - 104.D0* DCOS( 2.D0* VAX - VAY)
6 - 104.D0* DCOS( VAX + 2.D0* VAY)
7 - 104.D0* DCOS( VAX - 2.D0* VAY)
C
8 - 104.D0* DCOS( 2.D0* VAY + VAZ)
9 - 104.D0* DCOS( 2.D0* VAY - VAZ)
A - 104.D0* DCOS( VAY + 2.D0* VAZ)
B - 104.D0* DCOS( VAY - 2.D0* VAZ)
C
C - 104.D0* DCOS( 2.D0* VAX + VAZ)
D - 104.D0* DCOS( 2.D0* VAX - VAZ)
E - 104.D0* DCOS( VAX + 2.D0* VAZ)
F - 104.D0* DCOS( VAX - 2.D0* VAZ)
C
S1 = S1
1 + 176.D0* DCOS( 2.D0* VAX)
2 + 176.D0* DCOS( 2.D0* VAY)
3 + 176.D0* DCOS( 2.D0* VAZ)
4 - 256.D0* DCOS( VAX) - 256.D0* DCOS( VAY) - 256.D0* DCOS( VAZ)
C
S2 = - 392.D0* ( RXF + RYF + RZF)
1 + ( - 11.D0* RXF - 11.D0* RYF - 11.D0* RZF)
2 * DCOS( 2.D0* VAX + 2.D0* VAY + 2.D0* VAZ)
C
3 + ( - 11.D0* RXF - 11.D0* RYF - 11.D0* RZF)
4 * DCOS( 2.D0* VAX + 2.D0* VAY - 2.D0* VAZ)
5 + ( - 11.D0* RXF - 11.D0* RYF - 11.D0* RZF)
C
6 * DCOS( 2.D0* VAX - 2.D0* VAY + 2.D0* VAZ)
7 + ( - 11.D0* RXF - 11.D0* RYF - 11.D0* RZF)
8 * DCOS( 2.D0* VAX - 2.D0* VAY - 2.D0* VAZ)
C
9 + ( - 24.D0* RXF - 24.D0* RYF + 40.D0* RZF)
A * DCOS( 2.D0* VAX + 2.D0* VAY + VAZ)
B + ( - 24.D0* RXF - 24.D0* RYF + 40.D0* RZF)
C * DCOS( 2.D0* VAX + 2.D0* VAY - VAZ)
C
D + ( - 24.D0* RXF - 24.D0* RYF + 40.D0* RZF)
E * DCOS( 2.D0* VAX - 2.D0* VAY + VAZ)
F + ( - 24.D0* RXF - 24.D0* RYF + 40.D0* RZF)
G * DCOS( 2.D0* VAX - 2.D0* VAY - VAZ)
H + ( - 24.D0* RXF + 40.D0* RYF - 24.D0* RZF)
C
I * DCOS( 2.D0* VAX + VAY + 2.D0* VAZ)
J + ( - 24.D0* RXF + 40.D0* RYF - 24.D0* RZF)
K * DCOS( 2.D0* VAX + VAY - 2.D0* VAZ)
L + ( - 24.D0* RXF + 40.D0* RYF - 24.D0* RZF)
M * DCOS( 2.D0* VAX - VAY + 2.D0* VAZ)
C
N + ( - 24.D0* RXF + 40.D0* RYF - 24.D0* RZF)
O * DCOS( 2.D0* VAX - VAY - 2.D0* VAZ)
P + ( 40.D0* RXF - 24.D0* RYF - 24.D0* RZF)
Q * DCOS( VAX + 2.D0* VAY + 2.D0* VAZ)
C
R + ( 40.D0* RXF - 24.D0* RYF - 24.D0* RZF)
S * DCOS( VAX + 2.D0* VAY - 2.D0* VAZ)
T + ( 40.D0* RXF - 24.D0* RYF - 24.D0* RZF)
U * DCOS( VAX - 2.D0* VAY + 2.D0* VAZ)
C
V + ( 40.D0* RXF - 24.D0* RYF - 24.D0* RZF)
W * DCOS( VAX - 2.D0* VAY - 2.D0* VAZ)
X + ( -68.D0* RXF -68.D0* RYF -58.D0* RZF)
Y * DCOS( 2.D0* VAX + 2.D0* VAY )
C
Z + ( -68.D0* RXF -68.D0* RYF -58.D0* RZF)
1 * DCOS( 2.D0* VAX - 2.D0* VAY )
2 + ( - 36.* RXF+160.* RYF + 48.D0* RZF)
3 * DCOS( 2.D0* VAX + VAY )
C
S2 = S2
1 + ( - 36.D0* RXF + 160.D0* RYF + 48.D0* RZF)
2 * DCOS( 2.D0* VAX - VAY)
3 + ( 160.D0* RXF - 36.D0* RYF + 48.D0* RZF)
4 * DCOS( VAX + 2.D0* VAY)
5 + ( 160.D0* RXF - 36.D0* RYF + 48.D0* RZF)
6 * DCOS( VAX - 2.D0* VAY)
C
7 + ( - 68.D0* RXF - 58.D0* RYF - 68.D0* RZF)
8 * DCOS( 2.D0* VAX + 2.D0* VAZ)
9 + ( - 68.D0* RXF - 58.D0* RYF -68.D0* RZF)
A * DCOS( 2.D0* VAX - 2.D0* VAZ)
B + ( - 36.D0* RXF + 48.D0* RYF+160.D0* RZF)
C * DCOS( 2.D0* VAX + VAZ)
C
D + ( - 36.D0* RXF + 48.D0* RYF + 160.D0* RZF)
E * DCOS( 2.D0* VAX - VAZ)
F + ( 160.D0* RXF + 48.D0* RYF - 36.D0* RZF)
G * DCOS( VAX + 2.D0* VAZ)
C
H + ( 160.D0* RXF + 48.D0* RYF - 36.D0* RZF)
I * DCOS( VAX - 2.D0* VAZ)
J + ( -58.D0* RXF - 68.D0* RYF - 68.D0* RZF)
K * DCOS( 2.D0* VAY + 2.D0* VAZ)
L + ( -58.D0* RXF - 68.D0* RYF - 68.D0* RZF)
M * DCOS( 2.D0* VAY - 2.D0* VAZ)
C
N + ( 48.D0* RXF - 36.D0* RYF + 160.D0* RZF)
O * DCOS( 2.D0* VAY + VAZ)
P + ( 48.D0* RXF - 36.D0* RYF + 160.D0* RZF)
Q * DCOS( 2.D0* VAY - VAZ)
C
R + ( 48.D0* RXF+160.D0* RYF - 36.D0* RZF)
S * DCOS( VAY + 2.D0* VAZ)
T + ( 48.D0* RXF + 160.D0* RYF - 36.D0* RZF)
U * DCOS( VAY - 2.D0* VAZ)
C
V + ( - 248.D0* RXF - 184.D0* RYF - 184.D0* RZF)
W * DCOS( 2.D0* VAX)
X + ( 640.D0* RXF + 72.D0* RYF + 72.D0* RZF)* DCOS( VAX)
C
Y + ( - 184.D0* RXF - 248.D0* RYF - 184.D0* RZF)
Z * DCOS( 2.D0* VAY)
1 + ( 72.D0* RXF + 640.D0* RYF + 72.D0* RZF)* DCOS( VAY)
C
S2 = S2
1 + ( - 184.D0* RXF - 184.D0* RYF - 248.D0* RZF)
2 * DCOS( 2.D0* VAZ)
3 + ( 72.D0* RXF + 72.D0* RYF + 640.D0* RZF)* DCOS( VAZ)
C
S2 = ( 1.D0 - T)* S2
C
S3 =
1 ( - 29.D0* BX - 29.D0* BY - 29.D0* BZ)
2 * DSIN( 2.D0* VAX + 2.D0* VAY + 2.D0* VAZ)
3 + ( - 29.D0* BX - 29.D0* BY + 29.D0* BZ)
4 * DSIN( 2.D0* VAX + 2.D0* VAY - 2.D0* VAZ)
C
5 + ( - 29.D0* BX + 29.D0* BY - 29.D0* BZ)
6 * DSIN( 2.D0* VAX - 2.D0* VAY + 2.D0* VAZ)
7 + ( - 29.D0* BX +29.D0* BY +29.D0* BZ)
8 * DSIN( 2.D0* VAX - 2.D0* VAY - 2.D0* VAZ)
C
9 + ( 44.D0* BX + 44.D0* BY +20.D0* BZ)
A * DSIN( 2.D0* VAX + 2.D0* VAY + VAZ)
B + ( 44.D0* BX + 44.D0* BY - 20.D0* BZ)
C * DSIN( 2.D0* VAX + 2.D0* VAY - VAZ)
C
D + ( 44.D0* BX -44.D0* BY +20.D0* BZ)
E * DSIN( 2.D0* VAX - 2.D0* VAY + VAZ)
F + ( 44.D0* BX -44.D0* BY - 20.D0* BZ)
G * DSIN( 2.D0* VAX - 2.D0* VAY - VAZ)
C
H + ( 44.D0* BX +20.D0* BY + 44.D0* BZ)
I * DSIN( 2.D0* VAX + VAY + 2.D0* VAZ)
J + ( 44.D0* BX + 20.D0* BY + 44.D0* BZ)
K * DSIN( 2.D0* VAX + 2.D0* VAY + 2.D0* VAZ)
C
L + ( 44.D0* BX +20.D0* BY -44.D0* BZ)
M * DSIN( 2.D0* VAX + VAY - 2.D0* VAZ)
N + ( 44.D0* BX - 20.D0* BY + 44.D0* BZ)
O * DSIN( 2.D0* VAX - VAY + 2.D0* VAZ)
P + ( 44.D0* BX - 20.D0* BY -44.D0* BZ)
C
Q * DSIN( 2.D0* VAX - VAY - 2.D0* VAZ)
R + ( 20.D0* BX + 44.D0* BY + 44.D0* BZ)
S * DSIN( VAX + 2.D0* VAY + 2.D0* VAZ)
C
T + ( 20.D0* BX + 44.D0* BY -44.D0* BZ)
U * DSIN( VAX + 2.D0* VAY - 2.D0* VAZ)
V + ( 20.D0* BX -44.D0* BY + 44.D0* BZ)
W * DSIN( VAX - 2.D0* VAY + 2.D0* VAZ)
C
X + ( 20.D0* BX -44.D0* BY -44.D0* BZ)
Y * DSIN( VAX - 2.D0* VAY - 2.D0* VAZ)
Z + ( -52.D0* BX - 52.D0* BY)* DSIN( 2.D0* VAX + 2.D0* VAY)
1 + ( -52.D0* BX + 52.D0* BY)* DSIN( 2.D0* VAX - 2.D0* VAY)
2 + ( 116.D0* BX + 40.D0* BY)* DSIN( 2.D0* VAX + VAY)
C
S3 = S3
1 + ( 116.D0* BX -40.D0* BY )* DSIN( 2.D0* VAX - VAY)
2 + ( 40.D0* BX + 116.D0* BY )* DSIN( VAX + 2.D0* VAY)
3 + ( 40.D0* BX - 116.D0* BY )* DSIN( VAX - 2.D0* VAY)
4 + ( -52.D0* BX - 52.D0* BZ)* DSIN( 2.D0* VAX + 2.D0* VAZ)
5 + ( -52.D0* BX + 52.D0* BZ)* DSIN( 2.D0* VAX - 2.D0* VAZ)
C
6 + ( 116.D0* BX + 40.D0* BZ)* DSIN( 2.D0* VAX + VAZ)
7 + ( 116.D0* BX - 40.D0* BZ)* DSIN( 2.D0* VAX - VAZ)
8 + ( 40.D0* BX + 116.D0* BZ)* DSIN( VAX + 2.D0* VAZ)
9 + ( 40.D0* BX - 116.D0* BZ)* DSIN( VAX - 2.D0* VAZ)
A + ( - 52.D0* BY - 52.D0* BZ)* DSIN( 2.D0* VAY + 2.D0* VAZ)
B + ( - 52.D0* BY + 52.D0* BZ)* DSIN( 2.D0* VAY - 2.D0* VAZ)
C
C + ( 116.D0* BY + 40.D0* BZ)* DSIN( 2.D0* VAY + VAZ)
D + ( 116.D0* BY - 40.D0* BZ)* DSIN( 2.D0* VAY - VAZ)
E + ( 40.D0* BY + 116.D0* BZ)* DSIN( VAY + 2.D0* VAZ)
F + ( 40.D0* BY - 116.D0* BZ)* DSIN( VAY - 2.D0* VAZ)
G + ( -72.D0* BX)* DSIN( 2.D0* VAX )
C
H + ( - 72.D0* BY)* DSIN( 2.D0* VAY)
I + ( - 72.D0* BZ)* DSIN( 2.D0* VAZ)
C
S3 = G* S3
C
S4 = S1
C
S5 = - 392.D0* ( RXF + RYF + RZF)
1 + ( - 11.D0* RXF - 11.D0* RYF - 11.D0* RZF)
2 * DCOS( 2.D0* VAX + 2.D0* VAY + 2.D0* VAZ)
3 + ( - 11.D0* RXF - 11.D0* RYF - 11.D0* RZF)
4 * DCOS( 2.D0* VAX + 2.D0* VAY - 2.D0* VAZ)
5 + ( - 11.D0* RXF - 11.D0* RYF - 11.D0* RZF)
6 * DCOS( 2.D0* VAX - 2.D0* VAY + 2.D0* VAZ)
C
7 + ( - 11.D0* RXF - 11.D0* RYF - 11.D0* RZF)
8 * DCOS( 2.D0* VAX - 2.D0* VAY - 2.D0* VAZ)
9 + ( - 24.D0* RXF - 24.D0* RYF + 40.D0* RZF)
A * DCOS( 2.D0* VAX + 2.D0* VAY + VAZ)
B + ( - 24.D0* RXF - 24.D0* RYF + 40.D0* RZF)
C * DCOS( 2.D0* VAX + 2.D0* VAY - VAZ)
D + ( - 24.D0* RXF - 24.D0* RYF + 40.D0* RZF)
E * DCOS( 2.D0* VAX - 2.D0* VAY + VAZ)
C
F + ( - 24.D0* RXF - 24.D0* RYF + 40.D0* RZF)
G * DCOS( 2.D0* VAX - 2.D0* VAY - VAZ)
H + ( - 24.D0* RXF + 40.D0* RYF - 24.D0* RZF)
I * DCOS( 2.D0* VAX + VAY + 2.D0* VAZ)
J + ( - 24.D0* RXF + 40.D0* RYF - 24.D0* RZF)
K * DCOS( 2.D0* VAX + VAY - 2.D0* VAZ)
L + ( - 24.D0* RXF + 40.D0* RYF - 24.D0* RZF)
M * DCOS( 2.D0* VAX - VAY + 2.D0* VAZ)
C
N + ( - 24.D0* RXF + 40.D0* RYF - 24.D0* RZF)
O * DCOS( 2.D0* VAX - VAY - 2.D0* VAZ)
P + ( 40.D0* RXF - 24.D0* RYF - 24.D0* RZF)
Q * DCOS( VAX + 2.D0* VAY + 2.D0* VAZ)
R + ( 40.D0* RXF - 24.D0* RYF - 24.D0* RZF)
S * DCOS( VAX + 2.D0* VAY - 2.D0* VAZ)
T + ( 40.D0* RXF - 24.D0* RYF - 24.D0* RZF)
U * DCOS( VAX - 2.D0* VAY + 2.D0* VAZ)
V + ( 40.D0* RXF - 24.D0* RYF - 24.D0* RZF)
W * DCOS( VAX - 2.D0* VAY - 2.D0* VAZ)
C
X + ( -68.D0* RXF -68.D0* RYF -58.D0* RZF)
Y * DCOS( 2.D0* VAX + 2.D0* VAY )
Z + ( -68.D0* RXF -68.D0* RYF -58.D0* RZF)
1 * DCOS( 2.D0* VAX - 2.D0* VAY )
2 + ( - 36.* RXF+160.* RYF + 48.D0* RZF)
3 * DCOS( 2.D0* VAX + VAY )
C
S5 = S5
1 + ( - 36.* RXF+160.* RYF + 48.D0* RZF)
2 * DCOS( 2.D0* VAX - VAY)
3 + ( 160.* RXF - 36.* RYF + 48.D0* RZF)
4 * DCOS( VAX + 2.D0* VAY)
5 + ( 160.* RXF - 36.* RYF + 48.D0* RZF)
6 * DCOS( VAX - 2.D0* VAY)
7 + ( -68.D0* RXF -58.D0* RYF -68.D0* RZF)
8 * DCOS( 2.D0* VAX + 2.D0* VAZ)
9 + ( -68.D0* RXF -58.D0* RYF -68.D0* RZF)
A * DCOS( 2.D0* VAX - 2.D0* VAZ)
B + ( - 36.* RXF + 48.D0* RYF+160.* RZF)
C * DCOS( 2.D0* VAX + VAZ)
D + ( - 36.* RXF + 48.D0* RYF+160.* RZF)
E * DCOS( 2.D0* VAX - VAZ)
F + ( 160.* RXF + 48.D0* RYF - 36.* RZF)
G * DCOS( VAX + 2.D0* VAZ)
H + ( 160.* RXF + 48.D0* RYF - 36.* RZF)
I * DCOS( VAX - 2.D0* VAZ)
J + ( -58.D0* RXF -68.D0* RYF -68.D0* RZF)
K * DCOS( 2.D0* VAY + 2.D0* VAZ)
M + ( -58.D0* RXF -68.D0* RYF -68.D0* RZF)
N * DCOS( 2.D0* VAY - 2.D0* VAZ)
O + ( 48.D0* RXF - 36.* RYF+160.* RZF)
P * DCOS( 2.D0* VAY + VAZ)
Q + ( 48.D0* RXF - 36.* RYF+160.* RZF)
R * DCOS( 2.D0* VAY - VAZ)
S + ( 48.D0* RXF+160.* RYF - 36.* RZF)
T * DCOS( VAY + 2.D0* VAZ)
U + ( 48.D0* RXF+160.* RYF - 36.* RZF)
V * DCOS( VAY - 2.D0* VAZ)
C
W + ( - 248.D0* RXF- 184.D0* RYF- 184.D0* RZF)* DCOS( 2.D0* VAX)
X + ( 640.D0* RXF +72.D0* RYF +72.D0* RZF)* DCOS( VAX )
Y + ( - 184.D0* RXF- 248.D0* RYF- 184.D0* RZF)* DCOS( 2.D0* VAY)
Z + ( 72.D0* RXF +640.D0* RYF +72.D0* RZF)* DCOS( VAY)
C
S5 = S5
1 + ( - 184.D0* RXF- 184.D0* RYF- 248.D0* RZF)* DCOS( 2.D0* VAZ)
2 + ( 72.D0* RXF +72.D0* RYF+640.D0* RZF)* DCOS( VAZ)
C
S5 = (-T)*S5
C
P1 = S1 + S2
P2 = S3
P3 = S4 + S5
P4 = 0.0
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MGHHQ2 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
RXF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RYF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RZF = 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 + (( BX + BY + BZ)* ( BX + BY + BZ))
1 /( 2.D0* ( RX + RY + RZ))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
C
S1 = 512.D0
1 - 72.D0* DCOS( 2.D0* VAX + VAY + 2.D0* VAZ)
2 - 72.D0* DCOS( 2.D0* VAX + VAY - 2.D0* VAZ)
3 - 72.D0* DCOS( 2.D0* VAX - VAY + 2.D0* VAZ)
4 - 72.D0* DCOS( 2.D0* VAX - VAY - 2.D0* VAZ)
C
5 + 80.D0* DCOS( 2.D0* VAX + VAY + VAZ)
6 + 80.D0* DCOS( 2.D0* VAX + VAY - VAZ)
7 + 80.D0* DCOS( 2.D0* VAX - VAY + VAZ)
8 + 80.D0* DCOS( 2.D0* VAX - VAY - VAZ)
9 + 80.D0* DCOS( VAX + VAY + 2.D0* VAZ)
A + 80.D0* DCOS( VAX + VAY - 2.D0* VAZ)
B + 80.D0* DCOS( VAX - VAY + 2.D0* VAZ)
C + 80.D0* DCOS( VAX - VAY - 2.D0* VAZ)
C
D - 208.* DCOS( 2.D0* VAX + VAY)
E - 208.* DCOS( 2.D0* VAX - VAY)
F + 320.D0* DCOS( VAX + VAY)
G + 320.D0* DCOS( VAX - VAY)
H + 64.D0* DCOS( 2.D0* VAX + 2.D0* VAZ)
I + 64.D0* DCOS( 2.D0* VAX - 2.D0* VAZ)
C
S1 = S1
1 - 208.* DCOS( VAY + 2.D0* VAZ)
2 - 208.* DCOS( VAY - 2.D0* VAZ)
3 + 320.D0* DCOS( VAY + VAZ)
4 + 320.D0* DCOS( VAY - VAZ)
C
5 + 256.* DCOS( 2.D0* VAX)
6 - 512.D0* DCOS( VAY)
7 + 256.D0* DCOS( 2.D0* VAZ)
C
S2 = ( 384.D0* RXF + 1280.D0* RYF + 384.D0* RZF)
1 + ( 48.D0* RXF - 80.D0* RYF + 48.D0* RZF)
C
2 * DCOS( 2.D0* VAX + VAY + 2.D0* VAZ)
3 + ( 48.D0* RXF -80.D0* RYF + 48.D0* RZF)
4 * DCOS( 2.D0* VAX + VAY - 2.D0* VAZ)
5 + ( 48.D0* RXF -80.D0* RYF + 48.D0* RZF)
C
6 * DCOS( 2.D0* VAX - VAY + 2.D0* VAZ)
7 + ( 48.D0* RXF -80.D0* RYF + 48.D0* RZF)
8 * DCOS( 2.D0* VAX - VAY - 2.D0* VAZ)
C
9 + ( - 120.D0* RXF )* DCOS( 2.D0* VAX + VAY + VAZ)
A + ( - 120.D0* RXF )* DCOS( 2.D0* VAX + VAY - VAZ)
B + ( - 120.D0* RXF )* DCOS( 2.D0* VAX - VAY + VAZ)
C + ( - 120.D0* RXF )* DCOS( 2.D0* VAX - VAY - VAZ)
C
D + ( - 120.D0* RZF)* DCOS( VAX + VAY + 2.D0* VAZ)
E + ( - 120.D0* RZF)* DCOS( VAX + VAY - 2.D0* VAZ)
F + ( - 120.D0* RZF)* DCOS( VAX - VAY + 2.D0* VAZ)
G + ( - 120.D0* RZF)* DCOS( VAX - VAY - 2.D0* VAZ)
C
H + ( 72.D0* RXF- 320.D0* RYF - 96.D0* RZF)
I * DCOS( 2.D0* VAX + VAY)
J + ( 72.D0* RXF- 320.D0* RYF - 96.D0* RZF)
K * DCOS( 2.D0* VAX - VAY)
L + ( 240.D0 * RZF)* DCOS( VAX + VAY)
M + ( 240.D0 * RZF)* DCOS( VAX - VAY)
C
N + ( - 96.D0* RXF+160.* RYF - 96.D0* RZF)
O * DCOS( 2.D0* VAX + 2.D0* VAZ)
P + ( - 96.D0* RXF+160.* RYF - 96.D0* RZF)
Q * DCOS( 2.D0* VAX - 2.D0* VAZ)
R + ( - 96.D0* RXF- 320.D0* RYF +72.D0* RZF)
S * DCOS( VAY + 2.D0* VAZ)
C
S2 = S2
1 + ( - 96.D0* RXF- 320.D0* RYF +72.D0* RZF)
2 * DCOS( VAY - 2.D0* VAZ)
3 + ( 240.D0 * RXF)* DCOS( VAY + VAZ)
4 + ( 240.D0 * RXF)* DCOS( VAY - VAZ)
5 + ( - 384.D0* RXF+640.D0* RYF+192.D0* RZF)* DCOS( 2.D0* VAX)
6 + ( - 144.D0* RXF- 1280.D0* RYF- 144.D0* RZF)* DCOS( VAY)
7 + ( 192.D0* RXF+640.D0* RYF- 384.D0* RZF)* DCOS( 2.D0* VAZ)
C
S2 = ( - 1.D0* ( 1.D0 - T))* S2
C
S3 = ( 88.* BX + 40.D0* BY +88.* BZ)
1 * DSIN( 2.D0* VAX + VAY + 2.D0* VAZ)
2 + ( 88.* BX + 40.D0* BY -88.* BZ)
3 * DSIN( 2.D0* VAX + VAY - 2.D0* VAZ)
4 + ( 88.* BX -40.D0* BY +88.* BZ)
5 * DSIN( 2.D0* VAX - VAY + 2.D0* VAZ)
C
6 + ( 88.* BX -40.D0* BY -88.* BZ)
7 * DSIN( 2.D0* VAX - VAY - 2.D0* VAZ)
8 + ( - 120.D0* BX -80.D0* BY -80.D0* BZ)
9 * DSIN( 2.D0* VAX + VAY + VAZ)
A + ( - 120.D0* BX -80.D0* BY +80.D0* BZ)
B * DSIN( 2.D0* VAX + VAY - VAZ)
C
C + ( - 120.D0* BX +80.D0* BY -80.D0* BZ)
D * DSIN( 2.D0* VAX - VAY + VAZ)
E + ( - 120.D0* BX +80.D0* BY +80.D0* BZ)
F * DSIN( 2.D0* VAX - VAY - VAZ)
G + ( -80.D0* BX -80.D0* BY- 120.D0* BZ)
H * DSIN( VAX + VAY + 2.D0* VAZ)
I + ( -80.D0* BX -80.D0* BY+120.D0* BZ)
J * DSIN( VAX + VAY - 2.D0* VAZ)
C
K + ( -80.D0* BX +80.D0* BY- 120.D0* BZ)
L * DSIN( VAX - VAY + 2.D0* VAZ)
M + ( -80.D0* BX +80.D0* BY+120.D0* BZ)
N * DSIN( VAX - VAY - 2.D0* VAZ)
O + ( 232.D0* BX +80.D0* BY)* DSIN( 2.D0* VAX + VAY)
P + ( 232.D0* BX -80.D0* BY)* DSIN( 2.D0* VAX - VAY)
Q + ( - 320.D0* BX- 320.D0* BY)* DSIN( VAX + VAY)
R + ( - 320.D0* BX+320.D0* BY)* DSIN( VAX - VAY)
S + ( - 96.D0* BX - 96.D0* BZ)* DSIN( 2.D0* VAX + 2.D0* VAZ)
T + ( - 96.D0* BX +96.D0* BZ)* DSIN( 2.D0* VAX - 2.D0* VAZ)
C
S3 = S3
1 + ( 80.D0* BY+232.D0* BZ)* DSIN( VAY + 2.D0* VAZ)
2 + ( 80.D0* BY- 232.D0* BZ)* DSIN( VAY - 2.D0* VAZ)
3 + ( - 320.D0* BY- 320.D0* BZ)* DSIN( VAY + VAZ)
4 + ( - 320.D0* BY+320.D0* BZ)* DSIN( VAY - VAZ)
5 + ( - 384.D0* BX )* DSIN( 2.D0* VAX )
6 + ( - 384.D0* BZ)* DSIN( 2.D0* VAZ)
C
S3 = G* S3
C
S4 = S1
C
S5 = ( 384.D0* RXF + 1280.D0* RYF + 384.D0* RZF)
1 + ( 48.D0* RXF -80.D0* RYF + 48.D0* RZF)
2 * DCOS( 2.D0* VAX + VAY + 2.D0* VAZ)
3 + ( 48.D0* RXF -80.D0* RYF + 48.D0* RZF)
4 * DCOS( 2.D0* VAX + VAY - 2.D0* VAZ)
5 + ( 48.D0* RXF -80.D0* RYF + 48.D0* RZF)
6 * DCOS( 2.D0* VAX - VAY + 2.D0* VAZ)
7 + ( 48.D0* RXF -80.D0* RYF + 48.D0* RZF)
8 * DCOS( 2.D0* VAX - VAY - 2.D0* VAZ)
9 + ( - 120.D0* RXF)* DCOS( 2.D0* VAX + VAY + VAZ)
A + ( - 120.D0* RXF)* DCOS( 2.D0* VAX + VAY - VAZ)
B + ( - 120.D0* RXF)* DCOS( 2.D0* VAX - VAY + VAZ)
C + ( - 120.D0* RXF)* DCOS( 2.D0* VAX - VAY - VAZ)
D + ( - 120.D0* RZF)* DCOS( VAX + VAY + 2.D0* VAZ)
E + ( - 120.D0* RZF)* DCOS( VAX + VAY - 2.D0* VAZ)
F + ( - 120.D0* RZF)* DCOS( VAX - VAY + 2.D0* VAZ)
G + ( - 120.D0* RZF)* DCOS( VAX - VAY - 2.D0* VAZ)
H + ( 72.D0* RXF- 320.D0* RYF - 96.D0* RZF)
J * DCOS( 2.D0* VAX + VAY)
K + ( 72.D0* RXF- 320.D0* RYF - 96.D0* RZF)
L * DCOS( 2.D0* VAX - VAY)
N + ( 240.D0 * RZF)* DCOS( VAX + VAY)
O + ( 240.D0 * RZF)* DCOS( VAX - VAY)
P + ( - 96.D0* RXF+160.* RYF - 96.D0* RZF)
Q * DCOS( 2.D0* VAX + 2.D0* VAZ)
R + ( - 96.D0* RXF+160.* RYF - 96.D0* RZF)
S * DCOS( 2.D0* VAX - 2.D0* VAZ)
T + ( - 96.D0* RXF- 320.D0* RYF +72.D0* RZF)
U * DCOS( VAY + 2.D0* VAZ)
C
S5 = S5
1 + (- 96.D0* RXF- 320.D0* RYF +72.D0* RZF)
2 * DCOS( VAY - 2.D0* VAZ)
3 + ( 240.D0 * RXF)* DCOS( VAY + VAZ)
4 + ( 240.D0 * RXF)* DCOS( VAY - VAZ)
5 + ( - 384.D0* RXF+640.D0* RYF+192.D0* RZF)* DCOS( 2.D0* VAX )
6 + ( - 144.D0* RXF- 1280.D0* RYF- 144.D0* RZF)* DCOS( VAY)
7 + ( 192.D0* RXF+640.D0* RYF- 384.D0* RZF)* DCOS( 2.D0* VAZ)
C
S5 = T* S5
C
P1 = S1 + S2
P2 = S3
P3 = S4 + S5
P4 = 0.D0
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MGRT ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C -----------------------------
F1( X) = DCOS( X)
F2( X) = DCOS( X/ SQ3)
F3( X) = DSIN( X)
F4( X) = DSIN( X/ SQ3)
C -----------------------------
C
SQ3 = DSQRT( 3.D0)
SQ2 = DSQRT( 2.D0)
FII = 10.D0* ( 11.D0 + 4.D0* SQ3 )/ 219.D0
C
RXF = FII* 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RYF = FII* 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RZF = FII* 0.5D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
C
G = 10.D0* SQ2* ( 3.D0* SQ3- 10.D0)/ 73.D0
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = FII* ( 1.D0 +
1 ( BX + BY + BZ)* ( BX + BY + BZ)
2 /( 2.D0* ( RX + RY + RZ)))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 10.D0* SQ2* ( 3.D0* SQ3 - 10.D0)/ 73.D0
END IF
C
S1 = 2 + F1( VAX)* F2( VAY) + F1( VAY)* F2( VAZ)
1 + F1( VAZ)* F2( VAX)
C
E1 = 24.D0 - 5.D0* SQ3
E2 = 4.D0 - SQ3
E3 = 5.D0 - 6.D0*SQ3
C
S2 = ( 3.D0/ 4.D0)* ( 1.D0 - T)
1 * ( - 3.D0* ( 11.D0 - 4.D0* SQ3 )* ( RXF + RYF + RZF)
C
2 + ( E1* RXF + E2* RYF + E3* RZF)* F1( VAX)* F2( VAY)
3 + ( E3* RXF + E1* RYF + E2* RZF)* F1( VAY)* F2( VAZ)
4 + ( E2* RXF + E3* RYF + E1* RZF)* F1( VAZ)* F2( VAX))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
D1 = 4.D0
D2 = 3.D0 + 2.D0* SQ3
C
S6 = ( DSQRT( 6.D0)/ 4.D0)* G
1 * ( D1* BX* F4( VAX)* F1( VAZ) + D2* BZ* F2( VAX)* F3( VAZ)
2 + D2* BY* F3( VAY)* F2( VAZ) + D1* BZ* F1( VAY)* F4( VAZ)
3 + D2* BX* F3( VAX)* F2( VAY) + D1* BY* F1( VAX)* F4( VAY))
C
S7 = S1
C
S8 = ( - 3.D0/ 4.D0)* T
1 * ( - 3.D0* ( 11.D0 - 4.D0* SQ3 )* ( RXF + RYF + RZF)
2 + ( E1* RXF + E2* RYF + E3* RZF) * F1( VAX)* F2( VAY)
3 + ( E3* RXF + E1* RYF + E2* RZF) * F1( VAY)* F2( VAZ)
4 + ( E2* RXF + E3* RYF + E1* RZF) * F1( VAZ)* F2( VAX))
C
S9 = 0.D0
S10 = 0.D0
S11 = 0.D0
S12 = 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MGTH1 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C ------------------------------
C
RXF = 0.25D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RYF = 0.25D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RZF = 0.25D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
C
G = 1.D0
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 0.5D0* ( 1.D0+ ( BX + BY + BZ)* ( BX + BY + BZ)
1 /( 2.D0* ( RX + RY + RZ)))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
C
S1 = 4.D0 + 2.D0* ( F2( VAX) + F2( VAY) + F2( VAZ))
S2 = - 40.D0* ( 1.D0 - T)* ( RXF* F1( VAX)
1 + RYF* F1( VAY) + RZF* F1( VAZ))
S3 = 0.D0
C
S4 = 0.D0
S5 = 0.D0
S6 = - 10.D0* G * ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ))
C
S7 = S1
S8 = 40.D0* T* ( RXF* F1( VAX) + RYF* F1( VAY) + RZF* F1( VAZ))
S9 = 0.D0
C
S10 = 0.D0
S11 = 0.D0
S12 = 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MGTH2 ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ---------------------------------
F1(X,Y) = DCOS( X)* DCOS( Y)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C ---------------------------------
C
RXF = 0.625D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RYF = 0.625D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
RZF = 0.625D0* ( BX + BY + BZ)* ( BX + BY + BZ)/ 3.D0
C
G = 1.D00
C
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.25D0* ( 1.D0 + ( BX + BY + BZ)* ( BX + BY + BZ)
1 /( 2.D0* ( RX + RY + RZ)))
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
C
S1 = 10.D0 + ( F2( VAX) + F2( VAY) + F2( VAZ))
1 + 4.D0* ( F1( VAX,VAY) + F1( VAY,VAZ) + F1( VAZ,VAX))
C
S2 = - 10.D0* ( 1.D0-T)* ( RXF* ( 3.D0 - 2.D0* F2( VAX)
1 - ( F1( VAX,VAY) - F1( VAY,VAZ) + F1( VAZ,VAX)))
2 + RYF* ( 3.D0 - 2.D0* F2( VAY)
C
3 - ( F1( VAX,VAY) + F1( VAY,VAZ) - F1( VAZ,VAX)))
4 + RZF* ( 3.D0 - 2.D0* F2( VAZ)
5 - ( - F1( VAX,VAY) + F1( VAY,VAZ) + F1( VAZ,VAX))))
C
S3 = 0.D0
S4 = 0.D0
S5 = 0.D0
C
S6 = - G* 5.D0* ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ)
1 + ( BX+BY)* F3( VAX+VAY) + ( BY+BZ)* F3( VAY+VAZ)
2 + ( BZ+BX)* F3( VAZ+VAX) + (-BX+BY)* F3(-VAX+VAY)
3 + ( -BY+BZ)* F3(-VAY+VAZ) + (-BZ+BX)* F3(-VAZ+VAX))
C
S7 = S1
C
S8 = 10.D0* T* ( RXF* ( 3.D0 - 2.D0* F2( VAX)
1 - ( F1( VAX,VAY) - F1( VAY,VAZ) + F1( VAZ,VAX)))
2 + RYF* ( 3.D0 - 2.D0* F2( VAY)
C
3 - ( F1( VAX,VAY) + F1( VAY,VAZ) - F1( VAZ,VAX)))
4 + RZF* ( 3.D0 - 2.D0* F2( VAZ)
5 - ( -F1( VAX,VAY) + F1( VAY,VAZ) + F1( VAZ,VAX))))
C
S9 = 0.D0
S10 = 0.D0
S11 = 0.D0
S12 = 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PCON
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /NUM/ MTS, NT, NRX, NBX, NKX, NRY, NBY, NKY, NRZ, NBZ, NKZ
C
COMMON /MTV/ MT(50), TV(14)
C
COMMON /ABC/ RVX(10), RVY(10), RVZ(10),
1 BVX(10), BVY(10), BVZ(10),
2 KMX(14), KMY(14), KMZ(14)
C
WRITE(7,2000) MTS,NT,NRX,NRY,NRZ,NBX,NBY,NBZ,NKX,NKY,NKZ
WRITE(6,2000) MTS,NT,NRX,NRY,NRZ,NBX,NBY,NBZ,NKX,NKY,NKZ
2000 FORMAT(/ 5X,' * No. of Methods = ',I3,5X,' * No. of Theta =',I3,
1 // 5X,' * No. of Rx = ',I3,' N_Ry=',I3,' N_Rz=',I3/
2 5X,' * No. of Bx = ',I3,' N_By=',I3,' N_Bz=',I3/
3 5X,' * No. of Kx = ',I3,' N_Ky=',I3,' N_Kz=',I3/ )
C
WRITE(6,2100) ( MT(I),I = 1, MTS )
WRITE(7,2100) ( MT(I),I = 1, MTS )
2100 FORMAT( 3X,' * METHOD = ',20I3 )
C
WRITE(6,2200) ( TV(I),I = 1, NT)
WRITE(7,2200) ( TV(I),I = 1, NT)
2200 FORMAT(/1X,'THETA = ',9F7.1,/,9X,9F7.1,/,9X,9F7.1,/,9X,
1 9F7.1,/,9X,9F7.1 )
C
WRITE(6,2300) ( RVX(I),I = 1, NRX)
WRITE(6,2400) ( RVY(I),I = 1, NRY)
WRITE(6,2500) ( RVZ(I),I = 1, NRZ)
C
WRITE(7,2300) ( RVX(I),I = 1, NRX)
WRITE(7,2400) ( RVY(I),I = 1, NRY)
WRITE(7,2500) ( RVZ(I),I = 1, NRZ)
2300 FORMAT(/3X,' RX = ',15(F6.3,2X))
2400 FORMAT( 3X,' RY = ',15(F6.3,2X))
2500 FORMAT( 3X,' RZ = ',15(F6.3,2X))
C
WRITE(6,2600) ( BVX(I),I = 1, NBX)
WRITE(6,2700) ( BVY(I),I = 1, NBY)
WRITE(6,2800) ( BVZ(I),I = 1, NBZ)
C
WRITE(7,2600) ( BVX(I),I = 1, NBX)
WRITE(7,2700) ( BVY(I),I = 1, NBY)
WRITE(7,2800) ( BVZ(I),I = 1, NBZ)
2600 FORMAT(/3X,' BX = ',15(F6.3,2X))
2700 FORMAT( 3X,' BY = ',15(F6.3,2X))
2800 FORMAT( 3X,' BZ = ',15(F6.3,2X))
C
WRITE(6,2900) ( KMX(I),I = 1, NKX)
WRITE(6,3000) ( KMY(I),I = 1, NKY)
WRITE(6,3100) ( KMZ(I),I = 1, NKZ)
C
WRITE(7,2900) ( KMX(I),I = 1, NKX)
WRITE(7,3000) ( KMY(I),I = 1, NKY)
WRITE(7,3100) ( KMZ(I),I = 1, NKZ)
2900 FORMAT(/2X,' KX = ',16(I3,2X),/,16X,10(I3,2X))
3000 FORMAT( 2X,' KY = ',16(I3,2X),/,16X,10(I3,2X))
3100 FORMAT( 2X,' KZ = ',16(I3,2X),/,16X,10(I3,2X)///)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PCON2 (I)
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /MTV/ MT(50), TV(14)
C
DIMENSION NMT( 32)
C
CHARACTER*6 NMT
C
C ----- FFEMs & FDMs ---
C
DATA NMT/ 'GHH ','GTH1 ','GTH2 ','GRT ','GCHH ','GCTH1 '
1 ,'GCTH2 ','GCRT ','MGHH ','MGTH1 ','MGTH2 ','MGRT '
2 ,'MGHHQ1','MGHHQ2'
C
3 ,'EXHH ','EXTH1 ','EXTH2 ','EXRT ','EXCHH ','EXCTH1'
4 ,'EXCTH2','EXCRT '
C
4 ,'LXWDF ','SZE ','EXFU ','EXFU3 ','FTCS ','FULLY '
5 ,'CRNSN ','RBWIS ','KOSLA ','IMFU ' /
C
WRITE(6,2000) MT(I), NMT( MT(I))
WRITE(7,2000) MT(I), NMT( MT(I))
2000 FORMAT(/1X,'METHOD = ',I3,2X,'# ',A6,'#')
C
RETURN
END
C **********************************************************************
C
SUBROUTINE RBWIS ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C -------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C -------------------------------
C
S1 = 2.D0
S2 = - 2.D0* ( RX* F1( VAX) + RY* F1( VAY) + RZ* F1( VAZ))
S3 = ( BX* F1( VAX) + BY* F1( VAY) + BZ* F1( VAZ))
C
S4 = 0.D0
S5 = 0.D0
S6 = - ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ))
C
S7 = 2.D0
S8 = 2.D0* ( RX* F1( VAX) + RY* F1( VAY) + RZ* F1( VAZ))
S9 = ( BX* F1( VAX) + BY* F1( VAY) + BZ* F1( VAZ))
C
S10= 0.D0
S11= 0.D0
S12= ( BX* F3( VAX) + BY* F3( VAY) + BZ* F3( VAZ))
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
BLOCK DATA
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /DEB/ IDEB(10)
COMMON /MTV/ MT(50), TV(14)
COMMON /NUM/ MTS, NT, NRX, NBX, NKX, NRY, NBY, NKY, NRZ, NBZ, NKZ
COMMON /ABC/ RVX(10), RVY(10), RVZ(10), BVX(10), BVY(10), BVZ(10),
1 KMX(14), KMY(14), KMZ(14)
C
DATA TV / 0.D0, 0.1D0, 0.2D0, 0.3D0, 0.4D0, 0.5D0, 0.6D0, 0.7D0,
1 0.8D0, 0.9D0, 1.D0, 1.1D0, 1.2D0, 1.3D0/
C
DATA NKX, NKY, NKZ / 14, 14, 14/
C
DATA KMX / -3, -4, -6, -8, - 10, - 20, -40, 3, 4, 6, 8, 10,20,40/
DATA KMY / -3, -4, -6, -8, - 10, - 20, -40, 3, 4, 6, 8, 10,20,40/
DATA KMZ / -3, -4, -6, -8, - 10, - 20, -40, 3, 4, 6, 8, 10,20,40/
C
C ***** (A) Check : OK **************** "'07-1 Version (Expl.)"********
C =====
C
* DATA MTS, NT / 4, 14/ ! Method No. / No. of THETA
C
* DATA IDEB/ 0, 9, 9, 0, 0, 0, 0, 5, 0, 0/
C
* DATA NRX, NRY, NRZ / 6, 6, 6/
* DATA RVX / 0.01, 0.02, 0.04, 0.1, 0.2, 0.3, 0., 0., 0., 0./
* DATA RVY / 0.01, 0.02, 0.04, 0.1, 0.2, 0.3, 0., 0., 0., 0./
* DATA RVZ / 0.01, 0.02, 0.04, 0.1, 0.2, 0.3, 0., 0., 0., 0./
C
* DATA NBX, NBY, NBZ/ 5, 1, 1/
* DATA BVX / 0.1, 0.2, 0.3, 0.4, 0.5D0, 0., 0., 0., 0., 0.D0/
* DATA BVY / 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.D0/
* DATA BVZ / 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.D0/
C
C ****** (B) Check : OK *** Expl. ************************************
C =======
C
DATA MTS, NT / 4, 14/ ! Method No. / No. of THETA
C
DATA IDEB/ 1, 1, 9, 0, 0, 0, 0, 5, 0, 0/
C
DATA NRX, NRY, NRZ / 6, 6, 6/
DATA RVX / 0.01, 0.02, 0.04, 0.1, 0.2, 0.3, 0., 0., 0., 0./
DATA RVY / 0.01, 0.02, 0.04, 0.1, 0.2, 0.3, 0., 0., 0., 0./
DATA RVZ / 0.01, 0.02, 0.04, 0.1, 0.2, 0.3, 0., 0., 0., 0./
C
DATA NBX, NBY, NBZ/ 4, 4, 4/
DATA BVX / 0.05, 0.1, 0.15, 0.2, 0., 0., 0., 0., 0., 0.D0/
DATA BVY / 0.05, 0.1, 0.15, 0.2, 0., 0., 0., 0., 0., 0.D0/
DATA BVZ / 0.05, 0.1, 0.15, 0.2, 0., 0., 0., 0., 0., 0.D0/
C
C ****** (C) Check : OK **** Implicit Method ********************
C =====
C
* DATA MTS, NT / 2, 14 /
C
* DATA IDEB / 0, 9, 9, 0, 0, 0, 0, 5, 0, 1/
C
* DATA NRX, NRY, NRZ / 4, 4, 4/
* DATA RVX / 0.D0, 0.1D0, 0.3D0, 1.D0, 0., 0., 0., 0., 0., 0./
* DATA RVY / 0.D0, 0.1D0, 0.3D0, 1.D0, 0., 0., 0., 0., 0., 0./
* DATA RVZ / 0.D0, 0.1D0, 0.3D0, 1.D0, 0., 0., 0., 0., 0., 0./
C
* DATA NBX, NBY, NBZ/ 7, 1, 1/
* DATA BVX / 0.1D0, 0.2D0, 0.3D0, 0.4D0, 0.5D0, 0.6D0, 0.7D0,
* 1 0., 0., 0.D0/
* DATA BVY / 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.D0/
* DATA BVZ / 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.D0/
C
C --------------------------------------------------------------------
C
END
C ***********************************************************************
C
SUBROUTINE DIPT
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CST/ PAI
COMMON /PRA/ ALFA
COMMON /DEB/ IDEB(10)
COMMON /MTV/ MT(50), TV(14)
COMMON /NUM/ MTS, NT, NRX, NBX, NKX, NRY, NBY, NKY, NRZ, NBZ, NKZ
C
C ----------------
C
CALL TTLE
C
C ----------------
C
PAI = 3.14159 26535 89793 D0
C ---------------------------------
C
C ----- Expl ---
C
MT(1) = 15
MT(2) = 18
MT(3) = 19
MT(4) = 22
C
C ----- Impl. ---
C
* MT(1) = 5
* MT(2) = 8
C ---------------------
C
ALFA = 1.D0
C ---------------------
C
WRITE(6,2000) ALFA
WRITE(7,2000) ALFA
2000 FORMAT(/ 3X,' * Alfa=',F5.0)
C
IF ( IDEB(1).EQ.0) NRY = 1
IF ( IDEB(1).EQ.0) NRZ = 1
C
IF ( IDEB(2).EQ.0) NBY = 1
IF ( IDEB(2).EQ.0) NBZ = 1
C
IF ( IDEB(3).EQ.0) NKY = 1
IF ( IDEB(3).EQ.0) NKZ = 1
C
C -------------------
C
CALL PCON
C
C -------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE OLD_SMTD ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, I, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /MTV/ MT(50), TV(14)
C
* DATA NMT/ 'GHH ','GTH1 ','GTH2 ','GRT ','GCHH ','GCTH1 '
* 1 ,'GCTH2 ','GCRT ','MGHH ','MGTH1 ','MGTH2 ','MGRT '
* 2 ,'MGHHQ1','MGHHQ2','EXHH ','EXTH1 ','EXTH2 ','EXRT '
* 3 ,'EXCHH ','EXCTH1','EXCTH2','EXCRT '
*
* 4 ,'LXWDF ','SZE ','EXFU ','EXFU3 ','FTCS ','FULLY ',
* 5 ,'CRNSN ','RBWIS ','KOSLA ','IMFU ' /
C
C
GO TO ( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
1 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
2 29, 30, 31, 32) , MT(I)
C
C
C ***** FEMs **********************************************
C
C ***** GM ***
C
1 CALL GHH ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
2 CALL GTH1 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
3 CALL GTH2 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
4 CALL GRT ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
C ***** GMC ***
C
5 CALL GCHH ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
6 CALL GCTH1 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
7 CALL GCTH2 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
8 CALL GCRT ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
C *** MGM ***
C
9 CALL MGHH ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
10 CALL MGTH1 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
11 CALL MGTH2 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
12 CALL MGRT ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
13 CALL MGHHQ1( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
14 CALL MGHHQ2( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
C
GO TO 99
C
C ***** EXPLICIT ***
C
15 CALL EXHH ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
C
GO TO 99
C
16 CALL EXTH1 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
C
GO TO 99
C
17 CALL EXTH2 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
C
GO TO 99
C
18 CALL EXRT ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
C
GO TO 99
C
19 CALL EXCHH ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
C
GO TO 99
C
20 CALL EXCTH1( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
C
GO TO 99
C
21 CALL EXCTH2( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
C
GO TO 99
C
22 CALL EXCRT ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
C
GO TO 99
C
C ***** FDMs **********************************************
C
C ----- Explicit FDMs ---
C
23 CALL LXWDF ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
C
GO TO 99
C
24 CALL SZE ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
C
GO TO 99
C
25 CALL EXFU ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
C
GO TO 99
C
26 CALL EXFU3 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
C
GO TO 99
C
27 CALL FTCS ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
C
GO TO 99
C
C ----- Implicit FDMs ---
C
28 CALL FULLY ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
GO TO 99
C
29 CALL CRNSN ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
C
GO TO 99
C
30 CALL RBWIS ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
C
GO TO 99
C
31 CALL KOSLA ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
C
GO TO 99
C
32 CALL IMFU ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
C
C -----------------------------------------------------------------
C
99 RETURN
END
C **********************************************************************
C
SUBROUTINE SMTD ( P1, P2, P3, P4, RX, RY, RZ, BX, BY,
1 BZ, VAX, VAY, VAZ, I, T)
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /MTV/ MT(50), TV(14)
C
C ----- FEM ------------------------------------------------
C
C *** GM ***
C
IF ( MT(I).EQ.1) GO TO 1
IF ( MT(I).EQ.2) GO TO 2
IF ( MT(I).EQ.3) GO TO 3
IF ( MT(I).EQ.4) GO TO 4
C
C *** GMC ***
C
IF ( MT(I).EQ.5) GO TO 5
IF ( MT(I).EQ.6) GO TO 6
IF ( MT(I).EQ.7) GO TO 7
IF ( MT(I).EQ.8) GO TO 8
C
C *** MGM ***
C
IF ( MT(I).EQ.9) GO TO 9
IF ( MT(I).EQ.10) GO TO 10
IF ( MT(I).EQ.11) GO TO 11
IF ( MT(I).EQ.12) GO TO 12
IF ( MT(I).EQ.13) GO TO 13
IF ( MT(I).EQ.14) GO TO 14
C
C *** EXPLICIT ***
C
IF ( MT(I).EQ.15) GO TO 15
IF ( MT(I).EQ.16) GO TO 16
IF ( MT(I).EQ.17) GO TO 17
IF ( MT(I).EQ.18) GO TO 18
IF ( MT(I).EQ.19) GO TO 19
IF ( MT(I).EQ.20) GO TO 20
IF ( MT(I).EQ.21) GO TO 21
IF ( MT(I).EQ.22) GO TO 22
C
C ----- FDM ------------------------------------------------
C
C ----- Explicit -------------
C
* IF ( MT(I).EQ.15.OR. MT(I).EQ.16) GO TO 15
* IF ( MT(I).EQ.17.OR. MT(I).EQ.18) GO TO 17
* IF ( MT(I).EQ.19.OR. MT(I).EQ.20) GO TO 19
* IF ( MT(I).EQ.21.OR. MT(I).EQ.22) GO TO 21
* IF ( MT(I).EQ.23.OR. MT(I).EQ.24) GO TO 23
C
C ----- Implicit -------------
C
* IF ( MT(I).EQ.25.OR. MT(I).EQ.26) GO TO 25
* IF ( MT(I).EQ.27.OR. MT(I).EQ.28) GO TO 27
* IF ( MT(I).EQ.29.OR. MT(I).EQ.30) GO TO 29
* IF ( MT(I).EQ.31.OR. MT(I).EQ.32) GO TO 31
* IF ( MT(I).EQ.33.OR. MT(I).EQ.34) GO TO 33
C
C ----- FEM ----------------------------
C
C *** GM ***
C
1 CALL GHH ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
2 CALL GTH1 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
3 CALL GTH2 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
4 CALL GRT ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
C *** GMC ***
C
5 CALL GCHH ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
6 CALL GCTH1 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
7 CALL GCTH2 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
8 CALL GCRT ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
C *** MGM ***
C
9 CALL MGHH ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
10 CALL MGTH1 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
11 CALL MGTH2 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
12 CALL MGRT ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
13 CALL MGHHQ1 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
14 CALL MGHHQ2 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ,T)
GO TO 99
C
C ----- EXPLICIT ---
C
15 CALL EXHH ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
GO TO 99
C
16 CALL EXTH1 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
GO TO 99
C
17 CALL EXTH2 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
GO TO 99
C
18 CALL EXRT ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
GO TO 99
C
19 CALL EXCHH ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
GO TO 99
C
20 CALL EXCTH1 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
GO TO 99
C
21 CALL EXCTH2 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
GO TO 99
C
22 CALL EXCRT ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
GO TO 99
C
C ----- FDM ------------------------------------------------
C
C ----- Explicit -----------------------
C
* 15 CALL LXWDF ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
* GO TO 99
* 17 CALL SZE ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
* GO TO 99
* 19 CALL EXFU ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
* GO TO 99
* 21 CALL EXFU3 ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
* GO TO 99
* 23 CALL FTCS ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
* GO TO 99
C
C ----- Implicit ----------------------
C
* 25 CALL FULLY ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
* GO TO 99
* 27 CALL CRNSN ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
* GO TO 99
* 29 CALL RBWIS ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ)
* GO TO 99
* 31 CALL KOSLA ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
* GO TO 99
* 33 CALL IMFU ( P1,P2,P3,P4,RX,RY,RZ,BX,BY,BZ,VAX,VAY,VAZ, I )
C
C ----------------------------------------------------------
C
99 RETURN
END
C **********************************************************************
C
SUBROUTINE SZE ( P1, P2, P3, P4, RX, RY, RZ, BX, BY, BZ,
1 VAX, VAY, VAZ, I )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /MTV/ MT(50), TV(14)
C
C ------------------------------
F1( X) = 1.D0 - DCOS( X)
F2( X) = DCOS( X)
F3( X) = DSIN( X)
C ------------------------------
C
IF ( MOD( MT(I),2).NE.0) GO TO 10
RXF = ( BX* BY + BY* BZ + BZ* BX)/ 3.D0
RYF = ( BX* BY + BY* BZ + BZ* BX)/ 3.D0
RZF = ( BX* BY + BY* BZ + BZ* BX)/ 3.D0
C
G = 1.D00
IF (( RX + RY + RZ).GE.0.001) THEN
F = 1.D0 + ( BX*BY + BY*BZ + BZ*BX)/( RX + RY + RZ)
RXF = RX* F
RYF = RY* F
RZF = RZ* F
G = 1.D00
END IF
GO TO 20
C
10 RXF = RX
RYF = RY
RZF = RZ
G = 1.D00
C
20 CONTINUE
S1 = 1.D0
S2 = - 2.D0* ( RXF* F1( VAX) + RYF* F1( VAY) + RZF* F1( VAZ))
C
S3 = G* ( - ( BX* BX + BY* BY + BZ* BZ)/ 4.D0
1 - 3.D0* ( BX + BY + BZ)/ 4.D0
2 + BX* F2( VAX) + BY* F2( VAY) + BZ* F2( VAZ)
C
3 + ( BX* BX - BX)* F2( 2.D0* VAX)/ 4.D0
4 + ( BY* BY - BY)* F2( 2.D0* VAY)/ 4.D0
5 + ( BZ* BZ - BZ)* F2( 2.D0* VAZ)/ 4.D0 )
C
S4 = 0.D0
S5 = 0.D0
C
S6 = G* (( BX* BX - 3.D0* BX)* F3( VAX)/ 2.D0
1 + ( BY* BY - 3.D0* BY)* F3( VAY)/ 2.D0
2 + ( BZ* BZ - 3.D0* BZ)* F3( VAZ)/ 2.D0
C
3 - ( BX* BX - BX)* F3( 2.D0* VAX)/ 4.D0
4 - ( BY* BY - BY)* F3( 2.D0* VAY)/ 4.D0
5 - ( BZ* BZ - BZ)* F3( 2.D0* VAZ)/ 4.D0 )
C
S7 = 1.D0
S8 = 0.D0
S9 = 0.D0
C
S10 = 0.D0
S11 = 0.D0
S12 = 0.D0
C
P1 = S1 + S2 + S3
P2 = S4 + S5 + S6
P3 = S7 + S8 + S9
P4 = S10 + S11 + S12
C
RETURN
END
C **********************************************************************
C
SUBROUTINE OUT1 ( RX, RY, RZ, BX, BY, BZ, KX, KY, KZ, AA, BA,
1 BC, TS, AB0)
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /DEB/ IDEB(10)
C
IF ( IDEB(4).EQ.9) GO TO 30
IF ( IDEB(6).EQ.9) GO TO 20
C
WRITE(6,2000) RX, BX, KX, AA, TS
WRITE(7,2000) RX, BX, KX, AA, TS
2000 FORMAT(' R=',F5.2,3X,'B=',F5.2,3X,'K=',I3,3X,'GZAI=',F9.6,3X,
1 'ET.ERR=' ,F9.6 )
GO TO 50
C
20 WRITE(6,2100) RX, BX, KX, BA, RY, BY, KY, BC, RZ, BZ, KZ, TS
WRITE(7,2100) RX, BX, KX, BA, RY, BY, KY, BC, RZ, BZ, KZ, TS
2100 FORMAT( 5X,'RX =',F5.2,3X,'BX =',F5.2,3X,'KX =',I3,
1 3X,' AMP.ERR. = ',F9.6/5X,'RY =',F5.2,3X,'BY =',F5.2,
2 3X,'KY =',I3,3X,' PHA.ERR. = ',F9.6 /5X,'RZ =',F5.2,3X,
3 'BZ =',F5.2,3X,'KZ =',I3,3X,' ET .ERR. = ',F9.6 / )
GO TO 50
C
30 IF ( IDEB(6).EQ.9) GO TO 40
WRITE(6,2200) RX, BX, KX, AA, AB0
WRITE(7,2200) RX, BX, KX, AA, AB0
2200 FORMAT(' R =',F5.2,3X,'B=',F5.2,3X,'K=',I3,3X,'GZAI=',F9.6,3X,
1 'EM.ERR=' ,F9.6 )
GO TO 50
C
40 WRITE(6,2300) RX, BX, KX, BA, RY, BY, KY, BC, RZ, BZ, KZ, AB0
WRITE(7,2300) RX, BX, KX, BA, RY, BY, KY, BC, RZ, BZ, KZ, AB0
2300 FORMAT( 5X,'RX =',F5.2,3X,'BX =',F5.2,3X,'KX =',I3,
1 3X,' AMP.ERR. = ',F9.6/5X,'RY =',F5.2,3X,'BY =',F5.2,3X,
2 'KY =',I3,3X,' PHA.ERR. = ',F9.6/5X,'RZ =',F5.2,3X,
3 'BZ =',F5.2,3X,'KZ =',I3,3X,' EM .ERR. = ',F9.6/ )
C
50 RETURN
END
C **********************************************************************
C
SUBROUTINE OUT2 ( RX, RY, RZ, BX, BY, BZ, RMXB, BSA, BSC, TBS,
1 AB1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /DEB/ IDEB(10)
C
IF ( IDEB(4).EQ.9) GO TO 30
IF ( IDEB(7).EQ.9) GO TO 20
C
WRITE(6,2000) RX, BX, RMXB, TBS
WRITE(7,2000) RX, BX, RMXB, TBS
2000 FORMAT(' R =',F6.3,3X,'B=',F5.2,3X,'GZAI=',F9.6,3X,'ET.ERR=',
1 F9.6)
GO TO 50
C
20 WRITE(6,2100) RX, BX, RMXB, RY, BY, BSA, RZ, BZ, BSC, TBS
WRITE(7,2100) RX, BX, RMXB, RY, BY, BSA, RZ, BZ, BSC, TBS
2100 FORMAT( 5X,'RX =',F5.2,3X,'BX =',F5.2,6X,'AMP.MAX. = ',F9.6,/
1 5X,'RY =',F5.2,3X,'BY =',F5.2,6X,'AMP.Err.Var. = ',F9.6,/
2 5X,'RZ =',F5.2,3X,'BZ =',F5.2,6X,'PHA.Err.Var. = ',F9.6,/
3 29X,' ET .Err.Var. = ',F9.6 /)
GO TO 50
C
30 IF ( IDEB(7).EQ.9) GO TO 40
WRITE(6,2200) RX, BX, RMXB, AB1
WRITE(7,2200) RX, BX, RMXB, AB1
2200 FORMAT(' R =',F5.2,3X,'B=',F5.2,3X,'GZAI=',F9.6,3X,'EM.ERR=',
1 F9.6)
GO TO 50
C
40 WRITE(6,2300) RX, BX, RMXB, RY, BY, BSA, RZ, BZ, BSC, AB1
WRITE(7,2300) RX, BX, RMXB, RY, BY, BSA, RZ, BZ, BSC, AB1
2300 FORMAT( 5X,'RX =',F5.2,3X,'BX =',F5.2,6X,'AMP.MAX. = ',F9.6,/
1 5X,'RY =',F5.2,3X,'BY =',F5.2,6X,'AMP.Err.Var. = ',F9.6,/
2 5X,'RZ =',F5.2,3X,'BZ =',F5.2,6X,'PHA.Err.Var. = ',F9.6,/
3 29X,' EM .Err.Var. = ',F9.6 )
C
50 RETURN
END
C **********************************************************************
C
SUBROUTINE OUT3 ( RX, RY, RZ, RMX, SGA, SGC, TVR, AB2)
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /DEB/ IDEB(10)
C
IF ( IDEB(4).EQ.9) GO TO 30
IF ( IDEB(8).EQ.9) GO TO 20
C
WRITE(6,2000) RX, RMX, TVR
WRITE(7,2000) RX, RMX, TVR
2000 FORMAT(' R =',F5.2,3X,'GZAI=',F9.6,3X,'ET.ERR=',F9.6)
GO TO 50
C
20 WRITE(6,2100) RX, RMX, RY, SGA, RZ, SGC, TVR
WRITE(7,2100) RX, RMX, RY, SGA, RZ, SGC, TVR
2100 FORMAT(// 5X,'RX =',F5.2,3X,' Amp.Max. = ',F9.6 /
1 5X,'RY =',F5.2,3X,' Amp.Err.Var. = ',F9.6 /
2 5X,'RZ =',F5.2,3X,' Pha.Err.Var. = ',F9.6 /
3 17X,' ET .Err.Var. = ',F9.6 /)
GO TO 50
C
30 IF ( IDEB(8).EQ.9) GO TO 40
WRITE(6,2200) RX, RMX, AB2
WRITE(7,2200) RX, RMX, AB2
2200 FORMAT(' R =',F5.2,3X,'GZAI=',F9.6,3X,'EM.ERR=',F9.6)
GO TO 50
C
40 CONTINUE
WRITE(6,2300) RX, RMX, RY, SGA, RZ, SGC, AB2
WRITE(7,2300) RX, RMX, RY, SGA, RZ, SGC, AB2
2300 FORMAT(/ 5X,'RX =',F5.2,3X,' AMP.MAX. = ',F9.6 /
1 5X,'RY =',F5.2,3X,' AMP.Err.Var. = ',F9.6 /
2 5X,'RZ =',F5.2,3X,' PHA.Err.Var. = ',F9.6 /
3 17X,' EM .Err.Var. = ',F9.6 ///)
C
50 RETURN
END
C **********************************************************************
C
SUBROUTINE OUT4 ( TSGA, TSGC, SUMT, AB3 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /DEB/ IDEB(10)
C
IF ( IDEB(4).EQ.9) GO TO 30
IF ( IDEB(9).EQ.9) GO TO 20
C
WRITE(6,2000) SUMT
WRITE(7,2000) SUMT
2000 FORMAT(' ET.Err.Var. = ',F9.6)
GO TO 900
C
20 WRITE(6,2100) TSGA, TSGC, SUMT
WRITE(7,2100) TSGA, TSGC, SUMT
2100 FORMAT(///15X,'*****************************************' /
1 15X,'** Total Ampl. Err.Var. = ',F9.6,' ** ' /
2 15X,'** Total Phase Err.Var. = ',F9.6,' ** ' /
3 15X,'** Total ET. Err.Var. = ',F9.6,' ** ' /
4 15X,'*****************************************' /)
GO TO 900
C
30 IF ( IDEB(9).EQ.9) GO TO 40
WRITE(6,2200) AB3
WRITE(7,2200) AB3
2200 FORMAT(' EM.Err.Var. = ',F9.6)
GO TO 900
C
40 WRITE(6,2300) TSGA, TSGC, AB3
WRITE(7,2300) TSGA, TSGC, AB3
2300 FORMAT(//15X,'*****************************************' /
1 15X,'** Total Ampl. Err.Var. = ',F9.6,' ** ' /
2 15X,'** Total Phase Err.Var. = ',F9.6,' ** ' /
3 15X,'** Total EM. Err.Var. = ',F9.6,' ** ' /
4 15X,'*****************************************' /)
C
900 RETURN
END
C **********************************************************************
C
SUBROUTINE OUT5 ( KR, OPET, OPR, OPT, OPGZ)
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION OPET(50), OPT(50), OPR(50), OPGZ(50)
C
WRITE(6,2000)
WRITE(7,2000)
2000 FORMAT(/' *** Optimum Time Scheme Ppr. *** ')
C
C --------------------------------------
DO 100 I = 1, KR
C
WRITE(6,2100) OPR(I), OPT(I), OPGZ(I), OPET(I)
WRITE(7,2100) OPR(I), OPT(I), OPGZ(I), OPET(I)
2100 FORMAT(/3X,'R =',F5.2,3X,'OPT. Time Sc. Par.=',F5.2,3X,
1 'GZAI Max.=',F9.6,3X,'ET=',F9.6 )
C
100 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************