5
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

開発のお役立ち情報Advent Calendar 2023

Day 21

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

Last updated at Posted at 2023-12-19

こんにちは!

Pythonで太陽系の天体位置を描画してみるという記事が面白かったので、次は、3Dシミュレーションしてみます。

今回扱う天体

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

やりたいこと

要件は以下の3つ。

  • 3Dシミュレーションする
  • 月と季節が分かるようにする
  • 様々な角度から見えるようにする

環境構築

Pythonの仮想環境を作ります。

python -m venv solar-system-orbits-3d
cd solar-system-orbits-3d

ディレクトリに移動した直後の状態

$ ls
bin/		include/	lib/		pyvenv.cfg

仮想環境の有効化

source bin/activate

必要なライブラリのインストール

pip install plotly astropy numpy

plotlyとは

データの可視化を行うライブラリです。
matplotlibよりも柔軟な対応ができます。
今回、matplotlibでは実現が難しい部分があったので、こちらを利用しました。

astropyとは

天文学の計算とデータ分析を行うライブラリです。
星や惑星の位置の計算、時間の変換、天文データの読み書き、単位変換など、天文学に特化した機能を備えてます。

コード

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()

実行結果

python main.py

ブラウザが起動します。
こちらは起動直後の状態。
スクリーンショット 2023-12-19 22.08.16.png
太陽と、水星、金星、地球、火星がプロットされてますね。

それでは、Playで再生してみましょう。(以降、Pauseで止めてスクショ)
プロットが始まり、軌道が描画されます。この時点で2月。まだ冬ですね。
スクリーンショット 2023-12-19 22.08.39.png

軌道が少し進んで、5月。
スクリーンショット 2023-12-19 22.09.20.png

続いて8月。
スクリーンショット 2023-12-19 22.09.46.png

そして11月。シミュレーションはここまでです。
11.png

続いて、様々な角度から見てみます。
11-2.png
11-3.png

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

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

まとめ

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

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

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

5
7
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
5
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?