小惑星の軌道要素の可視化

前回の続きです。

読み込んだ小惑星のデータセットが妥当であるか、またモデルから生成した小惑星がそれらしいかを確かめるため、小惑星の軌道を可視化する簡単なプログラムを書きました。

https://gist.github.com/msakuta/da721b58ca08c3564075

小惑星番号順に最初の1000個までの軌道をプロットしてみます。赤、緑、青の直線はそれぞれx, y, z軸を表し、原点は太陽の位置です。

軌道を消し、ある時刻における小惑星の位置だけをプロットすると下図となります。

軌道要素から時間発展を計算するには、このへんを参照しながら、以下のような計算を用いました。

まず周期 $T$ を以下のように求めます。

$$T = 2 \pi \sqrt{ \frac{a^3}{GM_{sun}} }$$

ここで、 $a$ は軌道長半径、 $G$ は万有引力定数、 $M_{sun}$ は太陽質量です。

ここから、平均近点角 $M$ を求めます。

$$M =  M_0 + (t - t_0) \frac{ 2 \pi }{T}$$

ここで、 $M_0$ は元期(epoch)における平均近点角、 $t$ は位置を求めたい時刻、 $t_0$ は元期の時刻です。

描画上は平均近点角を離心近点角 $E$ に変換する必要がありますが、離心率  $e$ を使って表される両者の関係

$$M = E - e \cdot \sin(E)$$

は $E$ について陽に解くことができませんので、漸化式 $E_{i+1} = M + e \sin(E_i)$ の初項まで使って

$$E \approx M + e \cdot \sin(M)$$

としました。粗い近似ですが、軌道要素から計算している以上、誤差が時間的に積算していくことはなく、1周期回れば近点にて正確な位置に戻ります。


次回は統計モデル（確率密度推定）から生成した小惑星を同じようにプロットしてみて、自然に見えるかどうかを確かめてみようと思います。