Return to 楕円軌道上の時刻ごとの位置を数値的に求める準備

Python で楕円軌道上の時刻ごとの位置を数値的に求めてアニメをつくる

「楕円軌道上の時刻ごとの位置を数値的に求める準備」で下ごしらえした式を Python で数値的に解く。

必要なモジュールの import

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.integrate import solve_ivp

# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']

無次元化された式

「楕円軌道上の時刻ごとの位置を数値的に求める準備」で下ごしらえした式（無次元化済み）。

scipy.integrate.solve_ivp() は $\displaystyle \frac{dy}{dt} = f(t, y)$ の形の式を解くので，変数を揃える。

$ \phi \rightarrow y, \ T \rightarrow t$

\begin{eqnarray}
\frac{dy}{dt} &=& f(t, y) = 2 \pi \frac{(1+e \cos y)^2}{(1-e^2)^{3/2}} \tag{1}\\
R(y) &\equiv& \frac{r}{a} = \frac{1-e^2}{1+e\cos y} \tag{2}\\
X(y) &\equiv& \frac{x}{a} = R(y) \cos y \tag{3}\\
Y(y) &\equiv& \frac{y}{a} = R(y) \sin y \tag{4}
\end{eqnarray}

In [2]:
# 微分方程式の右辺 独立変数は t, y 決め打ち
def f(t, y):
    global e
    return 2*np.pi*(1+e*np.cos(y))**2/((1-e**2)**(3/2))

def R(y):
    global e
    return (1-e**2)/(1+e*np.cos(y))

def Xdaen(y):
    global e
    return R(y)*np.cos(y)

def Ydaen(y):
    global e
    return R(y)*np.sin(y)

# 離心率。かなり大きめ
e = 0.9

solve_ivp() で常微分方程式を数値的に解く

In [3]:
# T の初期値 
t0 = 0
# T の終了値
t1 = 1
t_span = [t0, t1]
# phi の初期値
y_ini = [0]

# 分割数
Ndiv = 36

t_list = np.linspace(t0, t1, Ndiv+1)
answer = solve_ivp(f, t_span, y_ini, 
                   t_eval = t_list, 
                   # 計算精度の設定
                   rtol = 1.e-13, atol = 1.e-15)

# 計算結果は .y として取り出します。
phii = answer.y[0]

X = Xdaen(phii)
Y = Ydaen(phii)

ちょうど1周期 $T = T_{\rm end}$ たてば $\phi = 2 \pi$ になっているはずだから，数値誤差を確認する。

$$\phi(T_{\rm end}) – 2 \pi$$

In [4]:
# [-1] は最後の項

phii[-1] - 2*np.pi
Out[4]:
-7.275335889289636e-11

一定時間間隔ごとの位置をグラフに

一定時間間隔 $\displaystyle \Delta T = \frac{1}{N_{\rm div}}$ ごとの位置：

In [5]:
fig = plt.figure(figsize=[6,6])
# 縦横のアスペクト比を等しく
plt.axes().set_aspect('equal')
# 軸の設定
plt.axis([-2, 1, -1, 1])

plt.scatter(X, Y, c="blue", s=20);

# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);

問題

  1. 離心率 $e$ を変えて数値計算を行う。
  2. 一定時間間隔ごとの位置をアニメにする。
    1. シーン1：軌道が楕円になる様子のアニメ。
    2. シーン2：楕円上に一定時間間隔ごとを位置を示すアニメ。
    3. シーン3：一定時間間隔ごとに掃く扇形を塗りつぶすアニメ。

アニメーション作成

滑らかにするため Ndiv は大きめにとって数値積分

In [6]:
# Ndiv = 720 として数値計算

# T の初期値 
t0 = 0
# T の終了値
t1 = 1
t_span = [t0, t1]
# phi の初期値
y_ini = [0]

# 分割数
Ndiv = 720

t_list = np.linspace(t0, t1, Ndiv+1)
answer = solve_ivp(f, t_span, y_ini, 
                   t_eval = t_list, 
                   # 計算精度の設定
                   rtol = 1.e-13, atol = 1.e-15)

# 計算結果は .y として取り出します。
phii = answer.y[0]

# 後々のことを考えてリストに変換しておく
X = Xdaen(phii).tolist()
Y = Ydaen(phii).tolist()

まずは途中までの軌道を描く

数値計算の結果は X と Y に格納されたので，適宜（間引きしたりして）グラフにする。

In [7]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()

# 外枠と目盛を非表示に
ax.axis('off')

# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')

# 軸の設定
plt.axis([-2, 1, -1, 1])

i = 1
# i*20 番目までの軌道
xi = X[:i*20]
yi = Y[:i*20]
plt.plot(xi, yi)

# i*20 番目の位置に赤丸
plt.plot(X[i*20], Y[i*20], "or", markersize = 3)

# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);

シーン1：FuncAnimation() でアニメーション作成

特定の i の値のところのグラフができたら，当該コードを def func(i): の中にコピぺして FuncAnimation() を呼べば，パラパラアニメをつくることができるんだったよね。

In [8]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()

# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')

def func(i):
    # 前の frame を消す
    plt.cla()

    # 外枠と目盛を非表示に
    ax.axis('off')

    # 軸の設定
    plt.axis([-2, 1, -1, 1])
    
    # i*20 番目までの軌道
    xi = X[:i*20]
    yi = Y[:i*20]
    plt.plot(xi, yi, linewidth = 0.5)

    # i*20 番目の位置に赤丸
    plt.plot(X[i*20], Y[i*20], "or", markersize = 3)

    # x軸 y軸
    plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
    plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);

    plt.title("ケプラーの第1法則：楕円軌道")

frames = 36
ani = FuncAnimation(fig, func,  
        interval = 100, 
        frames = range(frames))

ani.save("Pd2anim01.mp4", dpi=300)
ani.save("Pd2anim01.gif", dpi=150)

楕円軌道全体に一定時間ごとの位置を描く

In [9]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()

# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')

# 外枠と目盛を非表示に
ax.axis('off')

# 軸の設定
plt.axis([-2, 1, -1, 1])

# 楕円軌道は全体
plt.plot(X, Y, linewidth = 0.5)

# i*20 番目までの 20 ごとの位置に赤丸
xi = X[:i*20:10]
yi = Y[:i*20:10]
plt.plot(xi, yi, "or", markersize = 3)

# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);

シーン2：FuncAnimation() でアニメーション作成

特定の i の値のところのグラフができたら，当該コードを def func(i): の中にコピぺして FuncAnimation() を呼べば，パラパラアニメをつくることができるんだったよね。

In [10]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()

# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')

def func(i):
    # 前の frame を消す
    plt.cla()

    # 外枠と目盛を非表示に
    ax.axis('off')

    # 軸の設定
    plt.axis([-2, 1, -1, 1])

    # 楕円軌道は全体
    plt.plot(X, Y, linewidth = 0.5)

    # i*20 番目までの 20 ごとの位置に赤丸
    xi = X[:i*20+1:20]
    yi = Y[:i*20+1:20]
    plt.plot(xi, yi, "or", markersize=3)

    # x軸 y軸
    plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
    plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)

    plt.title("ケプラーの第2法則：一定時間間隔ごとの位置")

frames = 36
ani = FuncAnimation(fig, func,  
        interval = 100, 
        frames = range(frames))
ani.save("Pd2anim02.mp4", dpi=300)
ani.save("Pd2anim02.gif", dpi=150)

Δt の間に掃く扇形を塗りつぶす

In [11]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()

# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')

# 外枠と目盛を非表示に
ax.axis('off')

# 軸の設定
plt.axis([-2, 1, -1, 1])

i = 0
# i 番目と i+1 番目でできる扇形を塗りつぶす
Xougi = [0] + X[i*20:(i+1)*20+1]
Yougi = [0] + Y[i*20:(i+1)*20+1]
plt.fill(Xougi, Yougi, fc = "yellow")

# 原点と i 番目を結ぶ直線
plt.plot([Xougi[0], Xougi[1]], 
         [Yougi[0], Yougi[1]], 'tab:blue', linewidth = 0.5)
# 原点と i+1 番目を結ぶ直線
plt.plot([Xougi[0], Xougi[-1]], 
         [Yougi[0], Yougi[-1]], 'tab:blue', linewidth = 0.5)

# 楕円軌道は全体
plt.plot(X, Y, linewidth = 0.5)

# 20 ごとの全ての位置に赤丸
plt.plot(X[::20], Y[::20], "or", markersize = 3)

# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);

シーン3：FuncAnimation() でアニメーション作成

特定の i の値のところのグラフができたら，当該コードを def func(i): の中にコピぺして FuncAnimation() を呼べば，パラパラアニメをつくることができるんだったよね。

In [12]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()

# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')

def func(i):
    # 前の frame を消す
    plt.cla()

    # 外枠と目盛を非表示に
    ax.axis('off')

    # 軸の設定
    plt.axis([-2, 1, -1, 1])

    # i 番目と i+1 番目でできる扇形を塗りつぶす
    Xougi = [0] + X[i*20:(i+1)*20+1]
    Yougi = [0] + Y[i*20:(i+1)*20+1]
    plt.fill(Xougi, Yougi, fc = "yellow")
    
    # 原点と i 番目を結ぶ直線
    plt.plot([Xougi[0], Xougi[1]], 
             [Yougi[0], Yougi[1]], 'tab:blue', linewidth = 0.5)
    # 原点と i+1 番目を結ぶ直線
    plt.plot([Xougi[0], Xougi[-1]], 
             [Yougi[0], Yougi[-1]], 'tab:blue', linewidth = 0.5)
    
    # 楕円軌道は全体
    plt.plot(X, Y, linewidth = 0.5)
    
    # 20 ごとの全ての位置に赤丸
    plt.plot(X[::20], Y[::20], "or", markersize = 3)

    # x軸 y軸
    plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
    plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)

    plt.title("ケプラーの第2法則：面積速度一定")

frames = 36
ani = FuncAnimation(fig, func,  
        interval = 100, 
        frames = range(frames))
ani.save("Pd2anim03.mp4", dpi=300)
ani.save("Pd2anim03.gif", dpi=150)

面積を比較する静止画パート

面積速度が一定であることをじっくり目視で確認するために，2つの面積を塗りつぶした静止画を最後に追加しておこう。

In [13]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()

# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')

def func(i):
    # 前の frame を消す
    plt.cla()

    # 外枠と目盛を非表示に
    ax.axis('off')

    # 軸の設定
    plt.axis([-2, 1, -1, 1])

    # 近点付近の扇形
    Xougi = [0] +  X[-21:-1] + X[0:21]
    Yougi = [0] +  Y[-21:-1] +Y[0:21]
    plt.fill(Xougi, Yougi, fc = "yellow")

    # 扇形の直線2本
    plt.plot([Xougi[0], Xougi[1]], 
             [Yougi[0], Yougi[1]], 'tab:blue', linewidth = 0.5)
    plt.plot([Xougi[0], Xougi[-1]], 
             [Yougi[0], Yougi[-1]], 'tab:blue', linewidth = 0.5)

    # 遠点付近の扇形
    Xougi = [0] + X[340:381]
    Yougi = [0] + Y[340:381]
    plt.fill(Xougi, Yougi, fc = "yellow")

    # 扇形の直線2本
    plt.plot([Xougi[0], Xougi[1]], 
             [Yougi[0], Yougi[1]], 'tab:blue', linewidth = 0.5)
    plt.plot([Xougi[0], Xougi[-1]], 
             [Yougi[0], Yougi[-1]], 'tab:blue', linewidth = 0.5)

    # 楕円軌道は全体
    plt.plot(X, Y, linewidth = 0.5)

    # 20 ごとの全ての位置に赤丸
    plt.plot(X[::20], Y[::20], "or", markersize = 3)

    # x軸 y軸
    plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
    plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)

    plt.title("ケプラーの第2法則：面積は等しい")

frames = 36
ani = FuncAnimation(fig, func,  
        interval = 100, 
        frames = range(frames))
ani.save("Pd2anim04.mp4", dpi=300)
ani.save("Pd2anim04.gif", dpi=150)

ffmpeg で動画を連結する

gif ファイルを連結する方法がうまくいかず，とりあえずは mp4 ファイルを連結して gif に変換することに。

In [14]:
# mp4 ファイルを連結する
dat = ["file Pd2anim01.mp4", 
       "file Pd2anim02.mp4", 
       "file Pd2anim03.mp4", 
       "file Pd2anim03.mp4", 
       "file Pd2anim04.mp4"]
np.savetxt('input.txt', dat, fmt='%s') # 文字列として書き込む

# すでに（古い）ファイルがある場合は削除
!rm -f Pd2outfile1234.mp4

# input.txt の内容にしたがって連結する
!ffmpeg -hide_banner -loglevel error -f concat -i input.txt -c copy Pd2outfile1234.mp4

# すでに（古い）ファイルがある場合は削除
!rm -f Pd2outfile1234.gif

# mp4 を gif に変換する
!ffmpeg -hide_banner -loglevel error -i Pd2outfile1234.mp4 -vf scale=1280:-1 Pd2outfile1234.gif