積分と人工衛星のスイングバイ

人工衛星のスイングバイ

前回ボールの投げ上げをシミュレーションしました。もう少し高度な問題についてシミュレーションしてみたいと思います。今回は人工衛星の軌道について計算してみます。
人工衛星の軌道の中でスイングバイという軌道があります。人工衛星の燃料は限られているため惑星などの重力を利用して加速・方向転換を行う軌道です。今回はこのような軌道を行うための速度などをシミュレーションで求めてみます。
jaxaスイングバイ

人工衛星の軌道シミュレーション

前回の記事でボールの投げ上げ問題をシミュレーションで計算してみました。人工衛星の軌道シミュレーションもボールと同じ加速度を使って少し先の未来を予測し続ける問題です。しかし今回の場合下記の２点を注意する必要があります。

　（１）重力の大きさが星と人工衛星の距離に応じて減少する
　（２）重力の方向が人工衛星から見た星の方向になる

そのためこの２つ（距離と方向x,y成分）を求める必要があります。

距離の算出

座標の差異から算出する、地球の位置（Ex,Ey）と人工衛星の位置（x,y）の差から算出できます。

 r=\sqrt{(E_x - x)^2 + (E_y - y)^2}

方向の算出

ベクトルのxy成分を使用して算出する。

\displaylines{
 r_x=(E_x - x)/|r| \\
 r_y=(E_y - y)/|r|
}

重力の大きさ

以上により、地球から人工衛星に働く重力加速度は下記のように定義できる。
GMは重力定数と地球の質量になります。

\displaylines{
 G_x=GM(E_x - x)/|r|^3\\
 G_y=GM(E_y - y)/|r|^3
}

シミュレーションのコード

簡単のため太陽・地球・人工衛星を３つの物体を用いて軌道をシミュレーションします。太陽を中心に固定し、その周りを地球を円運動させた簡易モデルを使用します。そのうえで太陽と地球から重力の影響を受けた衛星の軌道を計算します。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import math

# 定数関連
StepCount = 1000
delta = 0.05
earth_r = 20.0
earth_GM = 10
sun_GM = 200

# 重力による加速度の計算
def cacl_acc(gm, diff_x , diff_y):
    r = math.sqrt(diff_x * diff_x + diff_y * diff_y)
    Ax = gm * diff_x / (r ** 3)
    Ay = gm * diff_y / (r ** 3)
    return Ax, Ay

if __name__ == '__main__':
    position_satellite = np.zeros((StepCount, 2)).astype(np.float32)  # (x, y)の組み合わせのペア
    position_earth = np.zeros((StepCount, 2)).astype(np.float32)  # (x, y)の組み合わせのペア

    #初期(位置・速度)
    x = earth_r + 1.0
    y = 0.0
    Vx = 2.0
    Vy = 2.2

    # 位置のシミュレーション（積分）
    for i in range(StepCount):
        # 地球の現在位置
        position_earth[i, 0] = earth_r * math.cos(i * 2 * math.pi / 360 * 2.0)
        position_earth[i, 1] = earth_r * math.sin(i * 2 * math.pi / 360 * 2.0)

        # 現在位置の設定
        position_satellite[i, 0] = x
        position_satellite[i, 1] = y

        # 地球による加速度の計算
        diff_x = position_earth[i, 0] - position_satellite[i, 0]
        diff_y = position_earth[i, 1] - position_satellite[i, 1]
        Ax_e, Ay_e = cacl_acc(earth_GM, diff_x, diff_y)

        # 太陽による加速度の計算
        diff_x = 0.0 - position_satellite[i, 0]
        diff_y = 0.0 - position_satellite[i, 1]
        Ax_s, Ay_s = cacl_acc(sun_GM, diff_x, diff_y)

        # 加速度の計算
        Ax = Ax_e + Ax_s
        Ay = Ay_e + Ay_s

        # 位置の更新
        x = x + Vx * delta
        y = y + Vy * delta

        # 速度の更新
        Vx = Vx + Ax * delta
        Vy = Vy + Ay * delta

        print('(x,y)=' + str(x) + ',' + str(y))


    fig = plt.figure()
    ims = []
    for i in range(StepCount):
        # (i * delta)秒までの位置を描写
        # 軌跡の描写
        im1 = plt.plot(position_satellite[0:i, 0], position_satellite[0:i, 1], linestyle='None', marker='.', color='grey')
        plt.xlim([-25.0, 25.0])
        plt.ylim([-25.0, 25.0])
        plt.subplot().set_aspect("equal", adjustable="box")

        # 衛星の現在位置
        x = position_satellite[i, 0]
        y = position_satellite[i, 1]
        im2 = plt.plot(x, y, color='green', marker='o')

        # 地球の現在位置
        x = position_earth[i, 0]
        y = position_earth[i, 1]
        im3 = plt.plot(x, y, color='blue', marker='o')

        # 太陽の現在位置
        x = 0.0
        y = 0.0
        im4 = plt.plot(x, y, color='red', marker='o')

        ims.append(im1 + im2 + im3 + im4)

    ani = animation.ArtistAnimation(fig, ims, interval=10)
    # GIFとして保存する場合
    #ani.save('anim.gif', writer='pillow')
    plt.show()

このコードを実行すると人工衛星（緑丸）の軌道をシミュレーションできます。人工衛星の軌道が地球（青丸）の軌道と重なるタイミングで大きく進行方向が変わるのを確認できると思います。

コード内の初期パラメータ（位置・速度）を変えることで人工衛星の軌道を違ったものに変えることができますので色々試してみてください。

    #初期(位置・速度)
    x = earth_r + 1.0
    y = 0.0
    Vx = 2.0
    Vy = 2.2

さいごに

かなりざっくりとした人工衛星の軌道シミュレーションを行いました。JAXAなどで実際に行われている軌道シミュレーションはこれと同じような計算を限りなく現実に近づけたものになります。