エクセルでシミュレーション ３次元で衛星の軌道計算
ルンゲ＝クッタ法（RK4法）を３次元へ適用 （２階微分方程式：３次元） スイングバイ あかつき
エクセルでどこまでできるか？
衛星の初速度を地球の公転速度程度に設定して軌道計算してみた。
'以下、計算部分を抜粋、（※別途、シート側にてx-y散布図作図が必要） ※データはＲ列〜Ｙ列
及び、ＡＪ列〜ＡＱ列にて印字
'パブリック又はプライベイト変数の宣言 x,y,z,tほか共通化 Privateはプロシージャ内の共通化
Public x,
y, z,
t
As Double '衛星の位置
Public vx, vy,
vz
As Double '衛星の速度
Private x_e_sat,
y_e_sat, z_e_sat As
Double '地球に対する位置
Private vx_e_sat,
vy_e_sat, vz_e_sat As Double
'地球に対する速度
Private x_e, y_e,
z_e
As Double '地球の位置
Private vx_e, vy_e,
vz_e
As Double '地球の速度
Private
r_e2sat
As Double '地球−衛星の距離
Public RTvX, RTvY,
RTvZ
As Double '引数の共通化
Public RTgX,
RTgY,
RTgZ
As Double '引数の共通化
Public dt As Double ' 通常は１時間 ＝ 3600 秒
'各種定数
Public Const
G As Double = 6.673 * 10 ^
(-11) ' 重力定数
Public Const
Msol As Double = 1.9891 * 10 ^
30 ' 太陽質量 kg
Public
Const M_e As Double = 5.9736 * 10 ^
24 ' 地球質量 kg
Public
Const R0 As Double = 1.4959787 * 10 ^
11
' 地球の公転半径 約1億5千万km
Public Const ω As
Double = 2 * 3.1415926535 / (365.256363 * 24 * 60 *
60)
' ２π／Ｔ Ｔ：地球の周期は365日
Public Const hill As Double = 150
* 10 ^
7
' HILL圏 150万km
Sub RK4_メインルーチン()
range("A1")="※データはＲ列〜Ｙ列
及び、ＡＪ列〜ＡＱ列にて印字⇒"
'◆◆◆◆◆◆◆◆◆◆◆◆◆ ルンゲクッタ法で 軌道計算 ◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆
'----------------------------------------------------------------
x
= 1.4959787 * 10 ^ 11 '衛星位置 x成分初期値
y =
0
'衛星位置 y成分初期値
z =
0
'衛星位置 z成分初期値
vx =
0
'衛星速度 x成分初期値
vy = 29.78 * 10 ^ 3 '衛星速度 y成分初期値
vz =
0
'衛星速度 z成分初期値
dt =
43200
'分割数 3600×12秒（テスト用）
cal_time = 12000000 '計算時間
x_e
= 1.4959787 * 10 ^ 11 '地球の位置
y_e =
0
'地球の位置
z_e =
0
'地球の位置
vx_e =
0
'地球の速度
vy_e = 29.78 * 10 ^ 3
'地球の速度 m/sec
vz_e = 7 * 10 ^
3
'地球の速度 （テスト用）
'----------------------------------------------------------------
For
t = 0 To cal_time Step dt 't：時間経過
０〜３６５日 (＝8760時間)
Call 衛星の位置
Call 衛星の速度
Call 地球の位置
Call 地球の速度
Call
ルンゲクッタ '衛星の軌道計算
Call
ルンゲクッタ_地球 '地球の軌道計算
Next t
End Sub
Sub 衛星の位置()
r = Sqr(x ^ 2 +
y ^ 2 + z ^ 2)
Cells(t / dt + 7, 18) = r
Cells(t / dt + 7, 19) = x
'x成分
Cells(t / dt + 7, 20) = y 'y成分
Cells(t / dt + 7, 21) = z 'z成分
End
Sub
Sub 衛星の速度()
v = Sqr(vx ^ 2 + vy ^ 2 + vz ^ 2)
Cells(t / dt + 7, 22)
= v
Cells(t / dt + 7, 23) = vx 'x成分
Cells(t / dt + 7, 24) = vy
'y成分
Cells(t / dt + 7, 25) = vz 'z成分
End Sub
Sub
衛星の地球に対する相対位置()
x_e_sat = x - x_e 'x成分
y_e_sat = y - y_e
'y成分
z_e_sat = z - z_e 'z成分
r_e2sat = Sqr(x_e_sat ^ 2 + y_e_sat ^ 2
+ z_e_sat ^ 2)
Cells(t / dt + 7, 26) = r_e2sat
Cells(t / dt + 7, 27) =
x_e_sat 'x成分
Cells(t / dt + 7, 28) = y_e_sat 'y成分
Cells(t / dt
+ 7, 29) = z_e_sat 'z成分
End Sub
Sub 地球の位置() ' 地球の位置
r_e =
Sqr(x_e ^ 2 + y_e ^ 2 + z_e ^ 2)
Cells(t / dt + 7, 36) = r_e
Cells(t / dt
+ 7, 37) = x_e
Cells(t / dt + 7, 38) = y_e
Cells(t / dt + 7, 39) =
z_e
End Sub
Sub 地球の速度()
v_e = Sqr(vx_e ^ 2 + vy_e ^ 2 + vz_e ^
2)
Cells(t / dt + 7, 40) = v_e
Cells(t / dt + 7, 41) =
vx_e 'x成分
Cells(t / dt + 7, 42) = vy_e
'y成分
Cells(t / dt + 7, 43) = vz_e 'z成分
End Sub
Sub
衛星の地球に対する相対速度()
vx_e_sat = vx - vx_e 'x成分
vy_e_sat = vy -
vy_e 'y成分
vz_e_sat = vz - vz_e 'z成分
v =
Sqr(vx_e_sat ^ 2 + vy_e_sat ^ 2 + vz_e_sat ^ 2)
Cells(t / dt + 7, 30) =
v
Cells(t / dt + 7, 31) = vx_e_sat 'x成分
Cells(t / dt + 7, 32) = vy_e_sat
'y成分
Cells(t / dt + 7, 33) = vz_e_sat 'z成分
End Sub
Sub
ルンゲクッタ()
'●●●●●●ルンゲ＝クッタ法 エンジン部分●●●●●●
Call f1_ENG(t, x, y, z, vx, vy, vz, RTvX, RTvY,
RTvZ)
Call f2_ENG(Msol, t, x, y, z, vx, vy, vz, RTgX, RTgY, RTgZ)
k1r_x =
dt * RTvX
k1r_y = dt * RTvY
k1r_z = dt * RTvZ
k1v_x = dt *
RTgX
k1v_y = dt * RTgY
k1v_z = dt * RTgZ
Call f1_ENG(t + dt / 2, x + k1r_x / 2, y + k1r_y /
2, z + k1r_z / 2, vx + k1v_x / 2, vy + k1v_y / 2, vz + k1v_z / 2, RTvX, RTvY,
RTvZ)
Call f2_ENG(Msol, t + dt / 2, x + k1r_x / 2, y + k1r_y / 2, z + k1r_z /
2, vx + k1v_x / 2, vy + k1v_y / 2, vz + k1v_z / 2, RTgX, RTgY, RTgZ)
k2r_x =
dt * RTvX
k2r_y = dt * RTvY
k2r_z = dt * RTvZ
k2v_x = dt *
RTgX
k2v_y = dt * RTgY
k2v_z = dt * RTgZ
Call f1_ENG(t + dt / 2, x + k2r_x / 2, y + k2r_y /
2, z + k2r_z / 2, vx + k2v_x / 2, vy + k2v_y / 2, vz + k2v_z / 2, RTvX, RTvY,
RTvZ)
Call f2_ENG(Msol, t + dt / 2, x + k2r_x / 2, y + k2r_y / 2, z + k2r_z /
2, vx + k2v_x / 2, vy + k2v_y / 2, vz + k2v_z / 2, RTgX, RTgY, RTgZ)
k3r_x =
dt * RTvX
k3r_y = dt * RTvY
k3r_z = dt * RTvZ
k3v_x = dt *
RTgX
k3v_y = dt * RTgY
k3v_z = dt * RTgZ
Call f1_ENG(t + dt, x + k3r_x, y + k3r_y, z +
k3r_z, vx + k3v_x, vy + k3v_y, vz + k3v_z, RTvX, RTvY, RTvZ)
Call
f2_ENG(Msol, t + dt, x + k3r_x, y + k3r_y, z + k3r_z, vx + k3v_x, vy + k3v_y, vz
+ k3v_z, RTgX, RTgY, RTgZ)
k4r_x = dt * RTvX
k4r_y = dt * RTvY
k4r_z =
dt * RTvZ
k4v_x = dt * RTgX
k4v_y = dt * RTgY
k4v_z = dt *
RTgZ
x = x + 1 / 6 * (k1r_x + 2 * k2r_x + 2 * k3r_x +
k4r_x)
y = y + 1 / 6 * (k1r_y + 2 * k2r_y + 2 * k3r_y + k4r_y)
z = z + 1 /
6 * (k1r_z + 2 * k2r_z + 2 * k3r_z + k4r_z)
vx = vx + 1 / 6 * (k1v_x + 2 * k2v_x + 2 * k3v_x +
k4v_x)
vy = vy + 1 / 6 * (k1v_y + 2 * k2v_y + 2 * k3v_y + k4v_y)
vz = vz +
1 / 6 * (k1v_z + 2 * k2v_z + 2 * k3v_z + k4v_z)
End Sub
Sub
ルンゲクッタ_地球()
'●●●●●●ルンゲ＝クッタ法 エンジン部分●●●●●●
Call f1_ENG(t, x_e, y_e, z_e, vx_e, vy_e, vz_e,
RTvX, RTvY, RTvZ)
Call f2_ENG(Msol, t, x_e, y_e, z_e, vx_e, vy_e, vz_e, RTgX,
RTgY, RTgZ)
k1r_x = dt * RTvX
k1r_y = dt * RTvY
k1r_z = dt *
RTvZ
k1v_x = dt * RTgX
k1v_y = dt * RTgY
k1v_z = dt * RTgZ
Call f1_ENG(t + dt / 2, x_e + k1r_x / 2, y_e +
k1r_y / 2, z_e + k1r_z / 2, vx_e + k1v_x / 2, vy_e + k1v_y / 2, vz_e + k1v_z /
2, RTvX, RTvY, RTvZ)
Call f2_ENG(Msol, t + dt / 2, x_e + k1r_x / 2, y_e +
k1r_y / 2, z_e + k1r_z / 2, vx_e + k1v_x / 2, vy_e + k1v_y / 2, vz_e + k1v_z /
2, RTgX, RTgY, RTgZ)
k2r_x = dt * RTvX
k2r_y = dt * RTvY
k2r_z = dt *
RTvZ
k2v_x = dt * RTgX
k2v_y = dt * RTgY
k2v_z = dt * RTgZ
Call f1_ENG(t + dt / 2, x_e + k2r_x / 2, y_e +
k2r_y / 2, z_e + k2r_z / 2, vx_e + k2v_x / 2, vy_e + k2v_y / 2, vz_e + k2v_z /
2, RTvX, RTvY, RTvZ)
Call f2_ENG(Msol, t + dt / 2, x_e + k2r_x / 2, y_e +
k2r_y / 2, z_e + k2r_z / 2, vx_e + k2v_x / 2, vy_e + k2v_y / 2, vz_e + k2v_z /
2, RTgX, RTgY, RTgZ)
k3r_x = dt * RTvX
k3r_y = dt * RTvY
k3r_z = dt *
RTvZ
k3v_x = dt * RTgX
k3v_y = dt * RTgY
k3v_z = dt * RTgZ
Call f1_ENG(t + dt, x_e + k3r_x, y_e + k3r_y, z_e +
k3r_z, vx_e + k3v_x, vy_e + k3v_y, vz_e + k3v_z, RTvX, RTvY, RTvZ)
Call
f2_ENG(Msol, t + dt, x_e + k3r_x, y_e + k3r_y, z_e + k3r_z, vx_e + k3v_x, vy_e +
k3v_y, vz_e + k3v_z, RTgX, RTgY, RTgZ)
k4r_x = dt * RTvX
k4r_y = dt *
RTvY
k4r_z = dt * RTvZ
k4v_x = dt * RTgX
k4v_y = dt * RTgY
k4v_z =
dt * RTgZ
x_e = x_e + 1 / 6 * (k1r_x + 2 * k2r_x + 2 * k3r_x
+ k4r_x)
y_e = y_e + 1 / 6 * (k1r_y + 2 * k2r_y + 2 * k3r_y + k4r_y)
z_e =
z_e + 1 / 6 * (k1r_z + 2 * k2r_z + 2 * k3r_z + k4r_z)
vx_e = vx_e + 1 / 6 * (k1v_x + 2 * k2v_x + 2 *
k3v_x + k4v_x)
vy_e = vy_e + 1 / 6 * (k1v_y + 2 * k2v_y + 2 * k3v_y +
k4v_y)
vz_e = vz_e + 1 / 6 * (k1v_z + 2 * k2v_z + 2 * k3v_z +
k4v_z)
End Sub
Sub f1_ENG(t, xx, yy, zz, vxx, vyy, vzz,
f1RTx, f1RTy, f1RTz) 'f1関数(t,位置ベクトル,速度ベクトル、戻り値速度ベクトル)
f1RTx = vxx
'単に速度ベクトルを返す
f1RTy = vyy '単に速度ベクトルを返す
f1RTz = vzz
'単に速度ベクトルを返す
End Sub
Sub f2_ENG(MM, t, xx, yy, zz, vxx, vyy, vzz, f2RTx,
f2RTy, f2RTz)
'f2関数(質量,t,位置ベクトル,速度ベクトル、戻り値加速度ベクトル)
'dv/dt=-Ｇ・Ｍ・(→r)／r^3 万有引力の法則
rr =
Sqr(xx ^ 2 + yy ^ 2 + zz ^ 2)
f2RTx = -G * MM / rr ^ 2 * (xx / rr)
'x成分 加速度ベクトル
f2RTy = -G * MM / rr ^ 2 * (yy / rr) 'y成分 加速度ベクトル
f2RTz
= -G * MM / rr ^ 2 * (zz / rr) 'z成分 加速度ベクトル
End Sub