太陽系の天体

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

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

やりたいこと

要件は以下の３つ。

太陽系の天体を描画する
どの天体か分かるようにする
10秒おきに再描画する

環境構築

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

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

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

$ ls
bin/		include/	lib/		pyvenv.cfg

仮想環境の有効化

source bin/activate

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

pip install matplotlib astropy numpy

matplotlibとは

データの可視化を行うライブラリです。
グラフやチャートの作成に使用され、多様なプロットスタイルやフォーマットを提供します。

astropyとは

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

コード

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)