[Fortran] 纯文本查看 复制代码
Module M_rk4_ods
Implicit None
Integer :: M = 36000
Contains
Subroutine Solve(Func, T0, Tt, Y0, N)
Implicit None
External Func
Integer :: N
Real *8 :: Y0(N), Y(N)
Real *8 :: K1(N), K2(N), K3(N), K4(N)
Real *8 :: H, Tt, T0, T
Integer :: I
H = (Tt-T0)/M
T = T0
Y = Y0
Do I = 1, M
Call Func(K1, T, Y, N)
Call Func(K2, T+H/2, Y+H/2*K1, N)
Call Func(K3, T+H/2, Y+H/2*K2, N)
Call Func(K4, T+H, Y+H*K3, N)
Y = Y + (K1+2*K2+2*K3+K4)*H/6
T = T0 + I*H
Write (11, 100) T, Y(1:6)
End Do
100 Format (1X, F12.6, 6F15.5)
End Subroutine Solve
Subroutine Fun1(F, T, Y, N)
Implicit None
Integer :: N
Real *8 :: F(N), Y(N)
Real *8 :: T, T0, Tt
Real *8 :: W, R, Fai, Fai0
Real *8 :: Ex, Ey, Ez, Bx, By, Bz
Real Q, W0, K
Complex I
I = (0, 1)
W0 = 200
Q = 100
Fai0 = 0.811
W = W0*Sqrt(1+((2*Y(6))/(W0**2))**2)
R = Y(6)*(1+((W0**2)/(2*Y(6)))**2)
Fai = Atan((W0**2)/(2*Y(6)))
Ex = (W0/W)*Exp((-(Y(4)**2+Y(5)**2)/(W**2)))*Exp(I*(T-Y(6)-Fai-Fai0-((Y(4)**2+Y(5)**2)/(2*R))))
Ey = 0
Ez = ((2*Y(4)**2)/(W**2*R))*Ex
Bx = ((4*Y(4)**2*Y(5)**2)/(W**4*R**2))*Ex
By = ((1-((W0**2)/(2*Y(6)*R))+((8*Y(6)**3*(Y(4)**2+Y(5)**2)*R)/((R-2*Y(6))**2))-((4*Y(4)**4)/(W**4*R**2)))+I*(((4*Y(6))/(W0**2*W**2))+((4*Y(4))/(W**2*R))-((8*Y(6)*(Y(4)**2+Y(5)**2))/(W0**2*W**4))))*Ex
Bz = -((2*Y(5)**2)/(W**2*R))*Ex
F(1) = Q*(Ex+(Y(2)/Sqrt(1+Y(2)**2))*Bz-(Y(3)/Sqrt(1+Y(3)**2))*By)
F(2) = Q*(Ey+(Y(3)/Sqrt(1+Y(3)**2))*Bx-(Y(1)/Sqrt(1+Y(1)**2))*Bz)
F(3) = Q*(Ez+(Y(1)/Sqrt(1+Y(1)**2))*By-(Y(2)/Sqrt(1+Y(2)**2))*Bx)
F(4) = Y(1)/Sqrt(1+Y(1)**2)
F(5) = Y(2)/Sqrt(1+Y(2)**2)
F(6) = Y(3)/Sqrt(1+Y(3)**2)
End Subroutine Fun1
End Module M_rk4_ods
Program Main
Use M_rk4_ods
Implicit None
Integer :: N
Real *8 :: Y(6), Y0(6)
Real *8 :: T0, Tt
Real Px0, Pz0, X0, Z0
Px0 = 6
Pz0 = 60
X0 = -1800
Z0 = -18000
Open (Unit=11, File='result0.txt')
N = 6
T0 = 0.0
Tt = 7200.0
Y0 = (/ Px0, 0, Pz0, X0, 0, Z0 /)
Call Solve(Fun1, T0, Tt, Y0, N)
End Program Main