Pythonで太陽系の天体位置を描画してみる

こんにちは！

[ISSが見える範囲に来たら通知する](https://qiita.com/yam_dev/items/3032ce4f23794fef498b)という記事を書いてみたら面白かったので、もう少し宇宙関連のコードを書いてみようと思います:slight_smile:

# 太陽系の天体
太陽系の天体は、以下の８つ。
水星、金星、地球、火星、木星、土星、天王星、海王星

この8つについて、位置を描画します。

# やりたいこと
要件は以下の３つ。
- 太陽系の天体を描画する
- どの天体か分かるようにする
- 10秒おきに再描画する

# 環境構築
Pythonの仮想環境を作ります。
```bash
python -m venv solar-system-orbits
cd solar-system-orbits
```
ディレクトリに移動した直後の状態
```bash
$ ls
bin/		include/	lib/		pyvenv.cfg
```
仮想環境の有効化
```bash
source bin/activate
```
必要なライブラリのインストール
```bash
pip install matplotlib astropy numpy
```
# matplotlibとは
データの可視化を行うライブラリです。
グラフやチャートの作成に使用され、多様なプロットスタイルやフォーマットを提供します。
# astropyとは
天文学の計算とデータ分析を行うライブラリです。
星や惑星の位置の計算、時間の変換、天文データの読み書き、単位変換など、天文学に特化した機能を備えてます。

# コード
```Python
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from astropy.time import Time
from astropy.coordinates import get_body_barycentric, solar_system_ephemeris
import numpy as np
import astropy.units as u
import itertools

# 惑星と軌道をプロットするための点の数を定義
planets = ['mercury', 'venus', 'earth', 'mars', 'jupiter', 'saturn', 'uranus', 'neptune']
# 各軌道をプロットするための点の数
num_points = 100
# 描画間隔
interval_time = 10000

# アニメーションが実行中かどうかを追跡するフラグ
is_running = True

# 惑星の位置と軌道を更新する関数
def update(frame):
    # ウィンドウが閉じられた場合はアニメーションを停止
    if not plt.get_fignums():
        global is_running
        is_running = False
        return
    
    # 現在の時間を取得
    current_time = Time.now()

    # ライブラリに組み込まれているデータを利用
    with solar_system_ephemeris.set('builtin'):
        # 現在の惑星の位置を取得
        planet_positions = {planet: get_body_barycentric(planet, current_time).xyz for planet in planets}

    # 描画のクリア
    ax.clear()
    # 軸のアスペクト比を等倍に設定
    ax.set_aspect('equal')
    # 太陽をプロット
    ax.plot(0, 0, 'yo', markersize=10, label='Sun')

    # 惑星の軌道と現在位置をプロット
    for planet in planets:
        # 軌道配列の初期化
        orbit = np.zeros((3, num_points))
        # 惑星の過去の位置を計算するための時間ステップを設定
        time_step = current_time - np.arange(num_points) * (1 / num_points) * 365 * 24 * 3600 * u.s
        # 指定された時間ステップでの惑星の位置を計算
        orbit = get_body_barycentric(planet, time_step).xyz.to(u.au).value

        # 惑星の軌道をプロット
        ax.plot(orbit[0], orbit[1], label=f'{planet.capitalize()} Orbit')
        # 惑星の現在位置をプロット
        ax.plot(planet_positions[planet][0].to(u.au).value, planet_positions[planet][1].to(u.au).value, 'o', markersize=5)

    # グラフに凡例を追加
    ax.legend()
    # タイトルをセット
    ax.set_title(f"Planetary Orbits and Positions on {current_time.to_value('iso', 'date')}")
    # グリッドを追加
    ax.grid(True)

frames = itertools.count()

# 図と軸の設定
fig, ax = plt.subplots(figsize=(10, 5))

# アニメーションを実行
ani = FuncAnimation(fig, update, frames=frames, interval=interval_time, blit=False, save_count=100)

try:
    # アニメーションが実行中の間だけメインループを実行
    plt.show()
except KeyboardInterrupt:
    plt.close(fig)
```

# 実行結果
描画までに少し時間がかかりましたが、太陽系が描画されました！
![スクリーンショット 2023-12-16 17.33.43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/407770/356f4613-1565-87d8-f8ee-1146ca828d14.png)

拡大ツールがあるので、拡大してみます。
![スクリーンショット 2023-12-16 17.36.08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/407770/69718018-2d3c-bad2-d340-7c374663c25b.png)
太陽、水星、金星、地球が表示されてますね。

位置が合ってるか、他のサイトも確認してみましょう。

[太陽系の天体の位置- STUDIO KAMADA](https://stdkmd.net/ssg/)
![スクリーンショット 2023-12-16 18.47.18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/407770/7ef77721-2bc9-ac92-d184-0f78e36af4d3.png)

合っているようですね！:smiley:
これだけ短いコードで実現できるのは、ライブラリのおかげですね。

# まとめ
今回は、太陽系の天体について、位置を描画してみました。
色々やってみると面白いので、他のライブラリも使ってみようと思います。

みなさんも一緒に頑張りましょう！