Pythonで太陽系の3Dシミュレーションをしてみる

こんにちは！

[Pythonで太陽系の天体位置を描画してみる](https://qiita.com/yam_dev/items/6f5cdd0e5857aee4e889)という記事が面白かったので、次は、3Dシミュレーションしてみます。

# 今回扱う天体
太陽系の天体は８つありますが、全部扱うと、太陽〜火星が観測しにくいです。
そこで今回は以下の天体に絞りました。
水星、金星、地球、火星

# やりたいこと
要件は以下の3つ。
- 3Dシミュレーションする
- 月と季節が分かるようにする
- 様々な角度から見えるようにする

# 環境構築
Pythonの仮想環境を作ります。
```bash
python -m venv solar-system-orbits-3d
cd solar-system-orbits-3d
```
ディレクトリに移動した直後の状態
```bash
$ ls
bin/		include/	lib/		pyvenv.cfg
```
仮想環境の有効化
```bash
source bin/activate
```
必要なライブラリのインストール
```bash
pip install plotly astropy numpy
```
# plotlyとは
データの可視化を行うライブラリです。
matplotlibよりも柔軟な対応ができます。
今回、matplotlibでは実現が難しい部分があったので、こちらを利用しました。
# astropyとは
天文学の計算とデータ分析を行うライブラリです。
星や惑星の位置の計算、時間の変換、天文データの読み書き、単位変換など、天文学に特化した機能を備えてます。

# コード
```Python
import plotly.graph_objs as go
from astropy.coordinates import get_body_barycentric_posvel, solar_system_ephemeris
from astropy.time import Time
import numpy as np
import astropy.units as u
from datetime import datetime

# 現在の日付から経過日数を計算（年初からではなく、前年の同日から）
today = datetime.now()
start_of_last_year = datetime(today.year - 1, today.month, today.day)
elapsed_days = (today - start_of_last_year).days

# アニメーションの設定
total_frames = 72

# 現在の月を取得
current_month = datetime.now().month

# 現在の日付から月の名前を取得
current_month_name = datetime.now().strftime("%B")

# 現在の季節を判定
if 3 <= current_month <= 5:
    current_season = "spring"
elif 6 <= current_month <= 8:
    current_season = "summer"
elif 9 <= current_month <= 11:
    current_season = "autumn"
else:
    current_season = "winter"

# 惑星のリストと色
planets = ['mercury', 'venus', 'earth', 'mars']
planet_colors = ['gray', 'orange', 'blue', 'red']

# 季節のインデックスに基づいて季節の名前を取得するための辞書
season_names = {0: "spring", 1: "summer", 2: "autumn", 3: "winter"}

# 季節ごとの背景色（辞書形式で指定）
season_colors = {
    "spring": 'rgba(255, 120, 120, 0.5)',  # 春 - 桜色（淡いピンク）
    "summer": 'rgba(120, 180, 120, 0.5)',  # 夏 - 新緑色（淡い緑）
    "autumn": 'rgba(255, 93, 56, 0.5)',   # 秋 - 紅葉色（淡いオレンジ）
    "winter": 'rgba(240, 240, 240, 0.5)'   # 冬 - 雪色（淡いグレー）
}

# 初期季節の背景色を設定
initial_bgcolor = season_colors[current_season]

season_titles = {
    "spring": "Spring",
    "summer": "Summer",
    "autumn": "Autumn",
    "winter": "Winter"
}

# フレームごとの季節を決定するための関数
def get_season_for_frame(frame_index, total_frames):
    season_index = (frame_index * 4) // total_frames
    return season_names[season_index]

# 軌道データを計算（前年の同日から現在日まで）
orbit_traces = []
with solar_system_ephemeris.set('builtin'):
    for planet in planets:
        # 現在の位置から過去1年間のデータを計算
        time_step = np.linspace(elapsed_days - 365, elapsed_days, total_frames) * u.day
        positions, _ = zip(*[get_body_barycentric_posvel(planet, Time.now() - t, ephemeris='builtin') for t in time_step])
        xyz = np.array([pos.xyz.to(u.au).value for pos in positions])
        orbit_traces.append(xyz)

# 太陽
sun = go.Scatter3d(x=[0], y=[0], z=[0], mode='markers', name='Sun', marker=dict(size=5, color='yellow'))

# 初期フレームのデータを作成（現在の位置から）
initial_data = [sun] + [
    go.Scatter3d(
        x=orbit[-1:, 0], y=orbit[-1:, 1], z=orbit[-1:, 2],
        mode='lines+markers',
        name=planet,
        line=dict(color=color),
        marker=dict(size=1, color=color)
    )
    for planet, orbit, color in zip(planets, orbit_traces, planet_colors)
]

# 初期季節の背景色とタイトルを設定
initial_bgcolor = season_colors[current_season]
initial_title = f"Solar System Planet Orbits - Animated: {current_month_name} / {season_titles[current_season]}"

# 初期カメラ設定（視点を調整）
initial_camera = dict(
    eye=dict(x=1.5, y=1.5, z=0.1),
    up=dict(x=0, y=0, z=1)
)

# 現在の季節に基づいて開始インデックスを計算
season_start_indexes = {"spring": 0, "summer": total_frames // 4, "autumn": total_frames // 2, "winter": 3 * total_frames // 4}
start_frame_index = season_start_indexes[current_season]

# 現在の季節のインデックス
current_season_index = list(season_names.values()).index(current_season)

# フレームごとに季節と月を更新するための関数
def get_season_and_month_for_frame(frame_index, total_frames):
    # 現在の月から季節を決定する関数
    def get_season_from_month(month):
        if 3 <= month <= 5:
            return "spring"
        elif 6 <= month <= 8:
            return "summer"
        elif 9 <= month <= 11:
            return "autumn"
        else:
            return "winter"

    # 月を計算（現在の月から始める）
    month_index = (current_month - 1 + frame_index * 12 // total_frames) % 12 + 1
    # 月に基づいて季節を決定
    season = get_season_from_month(month_index)

    return season, month_index

# アニメーションの各フレームを作成
frames = [
    go.Frame(
        data=[sun] + [
            go.Scatter3d(
                x=orbit[-i-1:, 0].tolist(),
                y=orbit[-i-1:, 1].tolist(),
                z=orbit[-i-1:, 2].tolist(),
                mode='lines+markers',
                name=planet,
                line=dict(color=color),
                marker=dict(size=1, color=color)
            )
            for planet, orbit, color in zip(planets, orbit_traces, planet_colors)
        ],
        layout=go.Layout(
            scene=dict(
                bgcolor=season_colors[get_season_and_month_for_frame(i, total_frames)[0]]
            ),
            title=f"Solar System Planet Orbits - Animated: {datetime(2023, get_season_and_month_for_frame(i, total_frames)[1], 1).strftime('%B')} / {season_titles[get_season_and_month_for_frame(i, total_frames)[0]]}"
        )
    )
    for i in range(total_frames)
]

# レイアウト設定
layout = go.Layout(
    title=initial_title,
    scene=dict(
        xaxis=dict(title='X AU'),
        yaxis=dict(title='Y AU'),
        zaxis=dict(title='Z AU'),
        bgcolor=initial_bgcolor,
        camera=initial_camera
    ),
    scene_aspectmode='cube',
    updatemenus=[dict(
        type="buttons",
        buttons=[
            dict(
                label="Play",
                method="animate",
                args=[None, {"frame": {"duration": 500, "redraw": True}, "fromcurrent": True, "repeat": True}]
            ),
            dict(
                label="Pause",
                method="animate",
                args=[[None], {"frame": {"duration": 0, "redraw": False}, "mode": "immediate"}]
            )
        ]
    )]
)

# フィギュアの作成
fig = go.Figure(data=initial_data, layout=layout, frames=frames)

# フィギュアの表示
fig.show()

```

実行結果
```bash
python main.py
```
ブラウザが起動します。
こちらは起動直後の状態。
![スクリーンショット 2023-12-19 22.08.16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/407770/faaaa6cb-c9a9-1219-c86f-64ff668ff8e6.png)
太陽と、水星、金星、地球、火星がプロットされてますね。

それでは、Playで再生してみましょう。（以降、Pauseで止めてスクショ）
プロットが始まり、軌道が描画されます。この時点で2月。まだ冬ですね。
![スクリーンショット 2023-12-19 22.08.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/407770/57e1143c-1d28-d4c2-3036-f6ed17d31c39.png)

軌道が少し進んで、5月。
![スクリーンショット 2023-12-19 22.09.20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/407770/6154f8b3-4334-af63-03db-b7251bb21cba.png)

続いて8月。
![スクリーンショット 2023-12-19 22.09.46.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/407770/4ca3f7c1-d6d7-8118-3536-47a2caaf3415.png)

そして11月。シミュレーションはここまでです。
![11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/407770/985e3ca1-1b37-e853-3b9f-4e086a7c25cf.png)


続いて、様々な角度から見てみます。
![11-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/407770/10976104-6bd5-44b6-8b10-69fa4fe29a71.png)
![11-3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/407770/ebced1e7-a8cb-c84c-3d87-7ce3fbe28bb0.png)

立体的に見れると面白いですね！

いずれ、銀河系にチャレンジしようと思います！

# まとめ
今回は、太陽系を3Dシミュレーションしてみました。
まだまだ改善点はありますが、こうやって動くと、すごくワクワクしますね:smiley:

宇宙にはたくさん面白いものがあるので、どんどんチャレンジしてみようと思います！

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