弱重力場中の近似解
シュバルツシルト時空中の粒子（天体，人工衛星等）の軌道は，重力場が弱いという近似のもと，以下のように書けることがわかった。
$$r = \frac{a(1-e^2)}{1 + e \cos(\gamma\phi) }$$
$$\gamma = \sqrt{ 1 -\frac{3 r_g}{a(1-e^2)}} \simeq 1 -\frac{3 r_g}{2a(1-e^2)} \equiv 1 – \Delta$$
$$r^2 \frac{d\phi}{d\tau} = \ell = \mbox{const.}$$
$\ell$ は測地線方程式を解く際に出てきた運動の定数で，角運動量の保存則に対応している。
参考：
ニュートン理論の場合に，楕円軌道上の時刻ごとの位置を数値的に求めてアニメをつくるのは，以下のページでやっているので参照：
これを応用して，近点移動がある場合のアニメーションを作成する。
周期
${\color{blue}{\Phi}} \equiv \gamma \phi$ とおくと，
$$r({\color{blue}{\Phi}}) = \frac{a(1-e^2)}{1 + e \cos({\color{blue}{\Phi}}) }$$
${\color{blue}{\Phi}} = 0$ のときに $r$ が最小値 $r_{\rm min} = a(1-e)$。
次に再び $r=r_{\rm min}$ になるまでの固有時間間隔を $P_{\tau}$ とすると，
\begin{eqnarray}
r^2({\color{blue}{\Phi}}) \frac{d{\color{blue}{\Phi}} }{d\tau} &=& \gamma \ell \\
\int_0^{2\pi} r^2( {\color{blue}{\Phi}} ) \,d {\color{blue}{\Phi}} &=& \gamma \ell \int_0^{P_{\tau}} d\tau \\
2 \pi a^2 \sqrt{1-e^2} &=& \gamma \ell P_{\tau} \\
\therefore \ \ P_{\tau} &=& \frac{2 \pi a^2 \sqrt{1-e^2}}{\gamma \ell}
\end{eqnarray}
無次元化
$P_{\tau}$を使って無次元化された固有時間 $T$ を
$$T \equiv \frac{\tau}{P_{\tau}}$$
とすると，
$$\frac{d{\color{blue}{\Phi}} }{dT} = 2 \pi \frac{(1+e \cos\phi)^2}{(1-e^2)^{3/2}} $$
これを数値的に解き，${\color{blue}{\Phi}} (T)$ がわかったら，$a$ で規格化した座標
\begin{eqnarray}
R(T) &\equiv& \frac{r}{a} = \frac{1-e^2}{1+e\cos{\color{blue}{\Phi}}(T)} \\
X(T) &\equiv& \frac{r}{a} \cos\phi \\
&=& \frac{1-e^2}{1+e\cos{\color{blue}{\Phi}}(T)}\cos\frac{{\color{blue}{\Phi}}(T)}{1 – \Delta} \\
&\simeq&\frac{1-e^2}{1+e\cos{\color{blue}{\Phi}}(T)}\cos\left(\left(1+\Delta\right) {\color{blue}{\Phi}}(T)\right) \\
Y(T) &\equiv& \frac{r}{a} \sin\phi \\
&=& \frac{1-e^2}{1+e\cos{\color{blue}{\Phi}}(T)}\sin\frac{{\color{blue}{\Phi}}(T)}{1 – \Delta} \\
&\simeq&\frac{1-e^2}{1+e\cos{\color{blue}{\Phi}}(T)}\sin\left(\left(1+\Delta\right) {\color{blue}{\Phi}}(T)\right)
\end{eqnarray}
でグラフにする。（厳密に言えば，シュバルツシルト時空なので平坦ではないが，近点移動を視覚的に表現するのには $x = r \cos\phi, \ y = r \sin\phi$ としたデカルト座標で描いてよろしいかと。）
Python でアニメーションを作成する
無次元化された式
scipy.integrate.solve_ivp() は $\displaystyle \frac{dy}{dt} = f(t, y)$ の形の式を解くので，変数を揃える。
$ {\color{blue}{\Phi}} = \gamma\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 \left( \left(1+\Delta \right) {y}\right) \tag{3}\\
Y(y) &\equiv& \frac{y}{a} = R(y) \sin \left( \left(1+\Delta \right) {y}\right) \tag{4}
\end{eqnarray}
必要なモジュールの import
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']
# 微分方程式の右辺 独立変数は 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 Xidou(y, Delta):
global e
return R(y)*np.cos((1+Delta)*y)
def Yidou(y, Delta):
global e
return R(y)*np.sin((1+Delta)*y)
# 離心率
e = 0.6
solve_ivp() で常微分方程式を数値的に解く
# T の初期値
t0 = 0
# T の終了値
t1 = 5
t_span = [t0, t1]
# phi の初期値
y_ini = [0]
# 1周期あたりコマ数
frames = 36
# コマの間の分割数（滑らかな曲線にするため）
ndiv = 10
# 1周期あたり分割数
Ndiv = frames * ndiv
t_list = np.linspace(t0, t1, 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]
近点移動がない楕円軌道の場合
グラフ
# 1周期ごとの色を c=colors[i] で色分けする
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
# 近点移動なし
Delta = 0
X = Xidou(phii, Delta)
Y = Yidou(phii, Delta)
# グラフのサイズの設定
fig = plt.figure(figsize=[6,4])
ax = fig.add_subplot()
# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')
# 外枠と目盛を非表示に
ax.axis('off')
# 軸の設定
plt.axis([-2, 1, -1, 1])
for j in range(1):
ini = Ndiv*j
end = Ndiv*(j+1)+1
# Delta T ごとの位置
plt.scatter(X[ini:end:ndiv], Y[ini:end:ndiv], s=5)
# 1周期分の軌道
plt.plot(X[ini:end], Y[ini:end], lw=2)
# 原点と近点を結ぶ点線
plt.plot([0, X[ini]], [0, Y[ini]],
c=colors[j], lw=0.5, linestyle = "dashed")
# 近点の位置
plt.scatter(X[ini], Y[ini], s=10, c=colors[j])
# x軸 y軸
plt.axhline(0, color='gray', dashes=(3, 3), linewidth=0.5)
plt.axvline(0, color='gray', dashes=(3, 3), linewidth=0.5);
アニメーション
# 近点移動なし
Delta = 0
X = Xidou(phii, Delta)
Y = Yidou(phii, Delta)
flist = []
for scene in range(t1):
# グラフのサイズの設定
fig = plt.figure(figsize=[6,4])
ax = fig.add_subplot()
# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')
def func(i):
global scene, frames, ndiv
# 前の frame を消す
plt.cla()
# 外枠と目盛を非表示に
ax.axis('off')
# 軸の設定
plt.axis([-2, 1, -1, 1])
# 1つ前までのシーンの軌道を色分けして描いておく
for j in range(scene):
ini = j*frames*ndiv
end = (j+1)*frames*ndiv + 1
# 軌道
plt.plot(X[ini:end], Y[ini:end],
lw=2, c=colors[j])
# Delta T ごとの位置
plt.scatter(X[ini:end:ndiv], Y[ini:end:ndiv],
s=5, c=colors[j])
# 原点と近点を結ぶ点線
plt.plot([0, X[ini]], [0, Y[ini]],
c=colors[j], lw=0.5, linestyle = "dashed")
# 近点の位置
plt.scatter(X[ini], Y[ini], s=10, c=colors[j])
# 本シーンの軌道を時々刻々と描く
ini = scene*frames*ndiv
end = ini + i*ndiv + 1
# Delta T ごとの位置
plt.scatter(X[ini:end:ndiv], Y[ini:end:ndiv],
s=5, c=colors[scene])
# i*10 までの軌道
plt.plot(X[ini:end], Y[ini:end],
lw=2, c=colors[scene])
# 原点と近点を結ぶ点線
plt.plot([0, X[ini]], [0, Y[ini]],
c=colors[scene], lw=0.5, linestyle = "dashed")
# 近点の位置
plt.scatter(X[ini], Y[ini], s=10, c=colors[scene])
# x軸 y軸
plt.axhline(0, color='gray', dashes=(3, 3), linewidth=0.5)
plt.axvline(0, color='gray', dashes=(3, 3), linewidth=0.5);
ani = FuncAnimation(fig, func,
interval = 100,
frames = range(frames+1))
fname = 'noidou%02d.mp4' % scene
ani.save(fname, dpi=300)
flist += ['file ' + fname]
ffmpeg で動画を連結する
# flist の内容確認
flist
# flist をファイルに書き込む
np.savetxt('input.txt', flist, fmt='%s') # 文字列として書き込む
# すでに（古い）ファイルがある場合は削除
!rm -f noidou.mp4
# input.txt の内容にしたがって連結する
!ffmpeg -hide_banner -loglevel error -f concat -i input.txt -c copy noidou.mp4
# 連結後は各シーンごとのファイル等を削除する
!rm -f noidou??.mp4
!rm -f input.txt
近点移動がある場合
グラフ
# 近点移動あり
Delta = 0.05
X = Xidou(phii, Delta)
Y = Yidou(phii, Delta)
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()
# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')
# 外枠と目盛を非表示に
ax.axis('off')
# 軸の設定
plt.axis([-2, 1, -2, 1])
for j in range(t1):
ini = Ndiv*j
end = Ndiv*(j+1)+1
# Delta T ごとの位置
plt.scatter(X[ini:end:ndiv], Y[ini:end:ndiv], s=5)
# 1周期分の軌道
plt.plot(X[ini:end], Y[ini:end], lw=2)
# 原点と近点を結ぶ点線
plt.plot([0, X[ini]], [0, Y[ini]],
c=colors[j], lw=0.5, linestyle = "dashed")
# 近点の位置
# zorder を大きめの値にして一番手前に近点を描く
plt.scatter(X[ini], Y[ini], s=20, c=colors[j], zorder=50)
# x軸 y軸
plt.axhline(0, color='gray', dashes=(3, 3), linewidth=0.5)
plt.axvline(0, color='gray', dashes=(3, 3), linewidth=0.5);
アニメーション
# 近点移動あり
Delta = 0.05
X = Xidou(phii, Delta)
Y = Yidou(phii, Delta)
flist = []
for scene in range(t1):
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()
# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')
def func(i):
global scene, frames, ndiv
# 前の frame を消す
plt.cla()
# 外枠と目盛を非表示に
ax.axis('off')
# 軸の設定
plt.axis([-2, 1, -2, 1])
# 1つ前までのシーンの軌道を色分けして描いておく
for j in range(scene):
ini = j*frames*ndiv
end = (j+1)*frames*ndiv + 1
# 軌道
plt.plot(X[ini:end], Y[ini:end],
lw=1, c=colors[j])
# 原点と近点を結ぶ点線
plt.plot([0, X[ini]], [0, Y[ini]],
c=colors[j], lw=0.5, linestyle = "dashed")
# 近点の位置
# zorder を大きめの値にして一番手前に近点を描く
plt.scatter(X[ini], Y[ini], s=10, c=colors[j], zorder=50)
# 本シーンの軌道を時々刻々と描く
ini = scene*frames*ndiv
end = ini + i*ndiv + 1
# Delta T ごとの位置
plt.scatter(X[ini:end:ndiv], Y[ini:end:ndiv],
s=5, c=colors[scene])
# i*10 までの軌道
plt.plot(X[ini:end], Y[ini:end],
lw=2, c=colors[scene])
# 原点と近点を結ぶ点線
plt.plot([0, X[ini]], [0, Y[ini]],
c=colors[scene], lw=0.5, linestyle = "dashed")
# 近点の位置
plt.scatter(X[ini], Y[ini], s=10, c=colors[scene], zorder=50)
# x軸 y軸
plt.axhline(0, color='gray', dashes=(3, 3), linewidth=0.5)
plt.axvline(0, color='gray', dashes=(3, 3), linewidth=0.5);
ani = FuncAnimation(fig, func,
interval = 100,
frames = range(frames+1))
fname = 'idou%02d.mp4' % scene
ani.save(fname, dpi=300)
flist += ['file ' + fname]
ffmpeg で動画を連結する
# flist の内容確認
flist
# flist をファイルに書き込む
np.savetxt('input.txt', flist, fmt='%s') # 文字列として書き込む
# すでに（古い）ファイルがある場合は削除
!rm -f idou.mp4
# input.txt の内容にしたがって連結する
!ffmpeg -hide_banner -loglevel error -f concat -i input.txt -c copy idou.mp4
# 連結後は各シーンごとのファイル等を削除する
!rm -f idou??.mp4
!rm -f input.txt