地球と超木星の重力シミュレーション

はじめに

このプロジェクトでは、Python を使って重力シミュレーションを開発し、地球と木星の質量の 500 倍の仮想的な「スーパー木星」との相互作用を示します。このシミュレーションは、働いている巨大な重力を考慮して、そのような巨大な天体が地球の運動に与える影響を示すことを目的としています。このプロジェクトは、物理学、天文学、および計算シミュレーションに情熱を持つ学生やホビーist に適しています。これを達成するために、数値計算に NumPy のような Python ライブラリを使用し、天体の動的な動きを視覚化するために Matplotlib を使用します。

👀 プレビュー

🎯 タスク

このプロジェクトでは、以下を学びます。

• 万有引力のニュートンの法則を理解して適用し、天体間の相互作用をモデル化する方法。
• Python プログラミングを使って重力システムの計算モデルを作成する方法。
• Python における効率的な数値計算のために NumPy ライブラリを使用する方法。
• 木星の質量の 500 倍の「スーパー木星」が存在する場合の地球の軌道力学をシミュレートする方法。
• シミュレーションの結果を分析して解釈し、巨大な天体が軌道力学に与える影響を理解する方法。
• Matplotlib を実装してシミュレーションの可視化表現を作成し、天体の軌道と相対位置を示す方法。
• 宇宙の文脈における力、質量、加速度の概念を探求する方法。
• 異なるシナリオに対して質量、距離、タイムステップなどのシミュレーションパラメータを微調整する方法。
• 科学的計算のための Python コードのデバッグと最適化のスキルを開発する方法。

🏆 成果

このプロジェクトを完了すると、以下のことができるようになります。

• 物理学の基本原理、特に万有引力のニュートンの法則を、実際の計算コンテキストで適用すること。
• Python を使って物理学ベースのシミュレーションを作成して実行すること。
• 大規模な数値計算を効率的に処理するための NumPy の使用に精通すること。
• Matplotlib を使って複雑なデータとシミュレーションを可視化し、科学的結果の解釈可能性を高めること。
• 惑星運動のダイナミクスと巨大な天体からの重力の影響を理解すること。
• シミュレーションの結果を分析して解釈し、天体力学に関する意味のある結論を導き出すこと。
• シミュレーションパラメータを調整して実験し、軌道力学の深い理解につながること。
• プログラミング環境、特に科学的計算の文脈における、向上した問題解決とデバッグスキルを示すこと。
• 重力が天体の運動をどのように形成するかに関する基礎知識を示し、天体物理学と計算モデリングのさらなる探求に道を開くこと。

環境のセットアップ

まず、プロジェクトファイルをセットアップし、必要なライブラリをインポートする必要があります。

simulation.py という名前の Python ファイルを作成します。

touch ~/project/simulation.py

次に、以下のライブラリをインポートします。

## 必要なライブラリをインポート
from typing import Tuple
import numpy as np
import pylab as py
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib.lines import Line2D
from tqdm import trange

重力シミュレーションプロジェクトのこの初期段階では、シミュレーションの環境をセットアップするためにいくつかの重要な Python ライブラリをインポートします。

1. **typing**：具体的には、typing モジュールからの Tuple がインポートされます。これは Python における型ヒントに使用され、コードの読みやすさとデバッグに役立ちます。タプルにグループ化できる値の型を指定することができ、関数定義の一貫性と明確性を保証します。
2. **numpy**：np としてインポートされる numpy ライブラリは、Python における科学的計算の基本パッケージです。大規模な多次元配列と行列のサポートを提供し、これらの配列に対して操作するための多数の高度な数学関数が用意されています。シミュレーションにおける数値計算を処理するために不可欠です。
3. **pylab**：py としてインポートされる pylab は、Matplotlib と一緒にインストールされる Matplotlib のモジュールです。MATLAB のようなインターフェイスを提供し、特に対話型の計算とプロットに便利です。ただし、一般的には明示的に matplotlib.pyplot と numpy をインポートすることが推奨されています。
4. **matplotlib.pyplot**：このモジュールは、plt としてインポートされ、Python において静的、アニメーション、対話型の可視化を作成するために使用されます。天体の軌道をプロットするためのこのプロジェクトの重要な要素です。
5. **matplotlib.animation**：Matplotlib のこのサブモジュールは、アニメーションを作成する機能を提供します。重力シミュレーションにおける惑星の動的な動きを視覚化するために不可欠です。
6. **matplotlib.lines.Line2D**：Matplotlib ライブラリのこのクラスは、プロットに描画できる線オブジェクトを作成するために使用されます。シミュレーションにおけるより詳細またはカスタマイズされたプロットに役立ちます。
7. **tqdm.trange**：tqdm モジュールからの trange は、組み込みの進捗メータ付きの Python の range 関数のバリアントです。シミュレーションにおけるループの進捗を表示するために役立ち、特に大規模な計算を行う場合に便利です。

これらのライブラリとモジュールをインポートすることで、重力シミュレーションプロジェクトに必要な複雑な数値計算と可視化の基盤を築きます。このセットアップは、実際の物理とアニメーションロジックが実装される後続のステップにとって重要です。

定数の定義

シミュレーションで使用する定数を定義します。これには、万有引力定数、天文単位、地球、木星、太陽の正規化質量が含まれます。

## 定数
G = 6.673e-11                 ## 万有引力定数
AU = 1.496e11                 ## キロメートル単位の天文単位
YEAR = 365*24*60*60.0         ## 1 年間の秒数
MM = 6e24                     ## 正規化質量
ME = 6e24/MM                  ## 地球の正規化質量
MS = 2e30/MM                  ## 太陽の正規化質量
MJ = 500*1.9e27/MM            ## 木星の正規化質量
GG = (MM*G*YEAR**2)/(AU**3)   ## シミュレーション用の万有引力定数

このステップでは、重力シミュレーションに必要ないくつかの重要な定数を定義します。これらの定数は、シミュレーションにおける天体力学と重力相互作用を正確にモデル化するための基本となります。

1. **万有引力定数 (G)**：SI 単位で 6.673e-11 に設定されており、ニュートンの万有引力の法則において重要な定数です。重力の強さを表し、2 つの質量間の重力を計算するために使用されます。
2. **天文単位 (AU)**：1.496e11 キロメートルと定義され、地球から太陽までの平均距離を表します。この単位は、太陽系内の距離をより扱いやすい形で表すために使用されます。
3. **年 (YEAR)**：365*24*60*60.0 と計算され、1 年間の秒数を表します。これは、時間に関する計算を年に変換するために使用され、天文シミュレーションにとってより直感的な時間単位です。
4. **正規化質量 (MM)**：6e24 に設定されており、シミュレーションで他の質量を正規化するための基準質量として使用されます。おおよそ地球の質量に相当します。
5. **地球の正規化質量 (ME)**：6e24/MM と計算され、正規化単位での地球の質量を表します。正規化質量 (MM) が地球の質量に基づいているため、基本的には 1 に設定されています。
6. **太陽の正規化質量 (MS)**：2e30/MM に設定されており、正規化単位での太陽の質量を表します。太陽の質量は、惑星にかかる重力をシミュレートする際の重要な要素です。
7. **木星の正規化質量 (MJ)**：このシミュレーションでは、木星の質量が 500*1.9e27/MM に拡大されており、「スーパー木星」のシナリオを表しています。このはるかに大きな質量は、はるかに大きな惑星が地球の軌道に与える影響を示すでしょう。
8. **シミュレーション用の万有引力定数 (GG)**：これは、シミュレーションのスケールと単位に調整された万有引力定数で、(MM*G*YEAR**2)/(AU**3) と計算されます。現実世界の万有引力定数を、シミュレーションで使用される正規化単位に適合させます。

これらの定数を定義することで、地球やこの仮想的な「スーパー木星」のような天体の重力と運動を正確にシミュレートするための必要なパラメータが揃います。このステップは、シミュレーションの計算が現実的な天文力学を反映することを確保するために不可欠です。同時に、最終的に異なる結果を観察するために、これらの定数を好きなだけ変更することができます。

重力関数の作成

次に、2 つの物体間の重力を計算する関数を定義します。この関数は、惑星の運動を決定する際に重要です。

## 重力を計算する関数
def gravitational_force(m1: float, m2: float, r: np.ndarray) -> np.ndarray:
    """
    2 つの物体間の重力を計算する。
    """
    F_mag = GG * m1 * m2 / (np.linalg.norm(r) + 1e-20)**2
    theta = np.arctan2(np.abs(r[1]), np.abs(r[0]) + 1e-20)
    F = F_mag * np.array([np.cos(theta), np.sin(theta)])
    F *= -np.sign(r)
    return F

このステップでは、2 つの天体間の重力を計算する関数 gravitational_force を定義します。この関数は、シミュレーションの重要な構成要素であり、空間内の任意の 2 つの質量間に働く力を決定するために、ニュートンの万有引力の法則を適用します。この関数の解説は以下の通りです。

1. 関数定義：関数 gravitational_force は 3 つの引数を取ります。
   • m1 と m2：2 つの物体の質量（浮動小数点数値）。
   • r：2 つの物体間の変位ベクトルを表す NumPy 配列。
2. **力の大きさの計算 (F_mag)**：
   • 関数はまず、式 F_mag = GG * m1 * m2 / (np.linalg.norm(r) + 1e-20)**2 を使用して重力の大きさを計算します。
   • GG はシミュレーション用の万有引力定数です。
   • np.linalg.norm(r) は 2 つの物体間のユークリッド距離を計算します。ゼロ除算を防ぐために小さな数 (1e-20) を加えます。
3. **方向の決定 (theta と F)**：
   • 角度 theta は、極座標における力ベクトルの角度を求めるために np.arctan2 を使用して計算されます。x 方向の変位 (r[0]) がゼロの場合にゼロ除算を避けるために、分母に小さな数 (1e-20) を加えます。
   • その後、力ベクトル F は、F_mag と角度 theta を使用して計算され、x 軸と y 軸に沿った成分 ([np.cos(theta), np.sin(theta)]) を持ちます。
4. 力の方向の調整：
   • 力ベクトル F は、重力の法則に従って常に引力的（他の物体に向かって指す）になるように -np.sign(r) で乗算されます。
5. 力ベクトルの返却：
   • 最後に、関数は重力ベクトル F を返します。

この関数を実装することで、計算環境で物理法則をどのように適用するかを学びます。この関数は、シミュレーションにおける地球と「スーパー木星」間の重力相互作用を計算するために不可欠であり、それらの軌道に影響を与えます。

RK4 ソルバの実装

運動の微分方程式を解くために、4 次ルンゲ・クッタ法（Runge-Kutta 4th order）のソルバを実装します。この方法は、常微分方程式を数値的に解くために広く使用されています。

## RK4 ソルバ
def RK4Solver(t: float, r: np.ndarray, v: np.ndarray, h: float, planet: str, r_other: np.ndarray, v_other: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    """
    惑星運動のための 4 次ルンゲ・クッタソルバ。
    """
    def dr_dt(v: np.ndarray) -> np.ndarray:
        return v

    def dv_dt(r: np.ndarray, planet: str) -> np.ndarray:
        if planet == 'earth':
            return (gravitational_force(ME, MS, r) + gravitational_force(ME, MJ, r - r_other)) / ME
        elif planet == 'jupiter':
            return (gravitational_force(MJ, MS, r) - gravitational_force(MJ, ME, r - r_other)) / MJ

    k11 = dr_dt(v)
    k21 = dv_dt(r, planet)

    k12 = dr_dt(v + 0.5 * h * k21)
    k22 = dv_dt(r + 0.5 * h * k11, planet)

    k13 = dr_dt(v + 0.5 * h * k22)
    k23 = dv_dt(r + 0.5 * h * k12, planet)

    k14 = dr_dt(v + h * k23)
    k24 = dv_dt(r + h * k13, planet)

    y0 = r + h * (k11 + 2 * k12 + 2 * k13 + k14) / 6
    y1 = v + h * (k21 + 2 * k22 + 2 * k23 + k24) / 6

    return y0, y1

このステップでは、RK4Solver という関数を実装します。この関数は、4 次ルンゲ・クッタ（RK4）法を使って微分方程式を解きます。この方法は、重力モデルにおける惑星の運動を正確にシミュレートするために重要です。以下は、RK4Solver 関数の概要です。

1. 関数定義：
   • RK4Solver はいくつかのパラメータを取ります。
     • t：現在の時間。
     • r：惑星の現在の位置ベクトル。
     • v：惑星の現在の速度ベクトル。
     • h：シミュレーションの時間ステップ。
     • planet：シミュレートされている惑星（地球または木星）を示す文字列。
     • r_other と v_other：他の惑星の位置と速度ベクトル。
2. RK4Solver 内のヘルパー関数：
   • dr_dt(v: np.ndarray) -> np.ndarray：この関数は、位置の変化率を返します。これは単に速度 v です。
   • dv_dt(r: np.ndarray, planet: str) -> np.ndarray：この関数は、重力による惑星の加速度を計算します。これは、先に定義された gravitational_force 関数を使って、太陽と他の惑星（地球または木星）によって働く力を計算し、加速度を返します。
3. ルンゲ・クッタの計算：
   • RK4 法では、時間ステップ内の異なる点で 4 つの「傾き」（k11, k21, k12, k22, k13, k23, k14, k24）を計算し、それらを組み合わせて、次の時間ステップでの位置と速度の正確な推定値を得ます。
   • これらの傾きは、現在の位置と速度、および位置の微分（dr_dt）と速度の微分（dv_dt）を使って計算されます。
4. 位置と速度の更新：
   • この関数は、RK4 の公式に従って、これらの傾きの加重平均を使って、惑星の次の位置（y0）と速度（y1）を計算します。
5. 結果の返却：
   • 最後に、この関数は、惑星の更新された位置と速度ベクトル（y0, y1）を返します。

RK4Solver 関数を実装することで、常微分方程式を数値的に解くための効率的で正確な方法を学びます。これは、惑星運動のような複雑なシステムを含むシミュレーションでは特に重要で、現実的な結果を得るために精度が鍵となるからです。RK4 法は、計算効率と精度のバランスをとるのに優れており、多くの科学および工学的応用で人気のある選択肢になっています。

アニメーションのセットアップ

シミュレーションを実行する前に、アニメーションのセットアップが必要です。このステップでは、プロットを作成し、惑星を表す線とマーカーを初期化します。

## アニメーションのセットアップ
def setup_animation() -> Tuple[py.Figure, py.Axes, Line2D, Line2D, py.Text]:
    """
    アニメーションプロットをセットアップする。
    """
    ## プロットの図と軸を作成
    fig, ax = py.subplots()

    ## 軸の範囲と目盛りを設定
    ax.axis('square')
    ax.set_xlim((-7.2, 7.2))
    ax.set_ylim((-7.2, 7.2))
    ax.get_xaxis().set_ticks([])
    ax.get_yaxis().set_ticks([])

    ## 太陽をプロット
    ax.plot(0, 0, 'o', markersize=9, markerfacecolor="#FDB813",
            markeredgecolor="#FD7813")

    ## 地球と木星の線を初期化
    line_earth, = ax.plot([], [], 'o-', color='#d2eeff',
                          markevery=10000, markerfacecolor='#0077BE', lw=2)
    line_jupiter, = ax.plot([], [], 'o-', color='#e3dccb', markersize=8,
                            markerfacecolor='#f66338', lw=2, markevery=10000)
    ## テキストオブジェクトを追加
    ttl = ax.text(0.24, 1.05, '', transform=ax.transAxes, va='center')

    ## コンポーネントを返す
    return fig, ax, line_earth, line_jupiter, ttl

このステップでは、重力シミュレーションのアニメーション用にプロットを準備する setup_animation という関数を定義します。この関数は、シミュレーションの可視コンポーネントをセットアップします。たとえば、プロット領域、惑星のマーカー、初期配置などです。

この関数は、地球と木星のシミュレートされた軌道を可視化するために不可欠です。これは、これらの惑星の位置がフレーム毎に更新されてアニメーションを作成するシミュレーションの動的な部分の舞台を築きます。

アニメーション関数の作成

アニメーションの各フレームにおいて惑星の位置を更新する関数を作成します。

## アニメーション関数
def animate(i: int) -> Tuple[Line2D, Line2D, py.Text]:
    """
    惑星運動のアニメーション関数。
    """
    earth_trail, jupiter_trail = 40, 200
    tm_yr = 'Elapsed time = ' + str(round(t[i], 1)) + ' years'
    ttl.set_text(tm_yr)
    line_earth.set_data(r[i:max(1, i - earth_trail):-1, 0],
                        r[i:max(1, i - earth_trail):-1, 1])
    line_jupiter.set_data(r_jupiter[i:max(
        1, i - jupiter_trail):-1, 0], r_jupiter[i:max(1, i - jupiter_trail):-1, 1])
    return line_earth, line_jupiter, ttl

このステップでは、惑星運動シミュレーションのアニメーション作成の核心となる animate 関数を定義します。この関数は、地球と木星の位置をプロット上で繰り返し更新することで、時間の経過に伴う動きの効果を作り出します。以下は animate 関数で起こることの説明です。

1. 地球と木星の軌跡の長さ：
   • earth_trail と jupiter_trail の変数が設定され、各惑星の過去の位置（軌跡）を表示する数が決まります。これにより、各惑星が通過した経路を視覚的な軌跡効果で表示します。
2. 経過時間の表示更新：
   • 経過時間（tm_yr）が計算され、ttl テキストオブジェクトのテキストとして設定されます。これにより、アニメーション上にシミュレーション時間を年単位で表示します。
3. 地球の位置の更新：
   • line_earth.set_data がアニメーションにおける地球の位置を更新します。これは、位置配列 r のスライスを使って軌跡効果を作成し、地球がその軌道上で最近どこにあったかを表示します。
4. 木星の位置の更新：
   • 同様に、line_jupiter.set_data が r_jupiter 配列のスライスを使って木星の位置を更新します。木星の軌跡の長さが長いことは、その軌道が大きいことを反映しています。
5. 更新されたオブジェクトの返却：
   • 関数は更新された Line2D オブジェクト（line_earth、line_jupiter）とテキストオブジェクト（ttl）を返します。これらはアニメーションの各フレームで変化する要素です。

この animate 関数は、地球と木星が相互の重力の影響下で時間の経過とともにどのように動くかを可視化するために不可欠です。各フレームでプロットを動的に更新することで、シミュレーションを生き生きとさせます。

シミュレーションパラメータの初期化

シミュレーションを開始する前に、時間パラメータ、地球と木星の位置および速度を初期化します。

## 初期化
ti, tf = 0, 120  ## 年単位の初期時刻と最終時刻
N = 100 * tf     ## 1 年あたり 100 点
t = np.linspace(ti, tf, N)  ## 時間配列
h = t[1] - t[0]  ## 時間ステップ

## 位置と速度の初期化
r = np.zeros([N, 2])         ## 地球の位置
v = np.zeros([N, 2])         ## 地球の速度
r_jupiter = np.zeros([N, 2])  ## 木星の位置
v_jupiter = np.zeros([N, 2])  ## 木星の速度

## 初期条件
r[0] = [1496e8 / AU, 0]
r_jupiter[0] = [5.2, 0]
v[0] = [0, np.sqrt(MS * GG / r[0, 0])]
v_jupiter[0] = [0, 13.06e3 * YEAR / AU]

このステップでは、地球と木星の運動のアニメーションを開始するためのシミュレーションの初期条件とパラメータを設定します。以下は行われていることの概要です。

1. 時間の初期化：
   • ti と tf は、シミュレーションの初期時刻と最終時刻として設定され、年単位で測定されます（0 から 120 年）。
   • N は、シミュレーションの総ポイント数として定義され、1 年あたり 100 点として計算されます。
   • t は、np.linspace を使用して作成された配列で、初期時刻から最終時刻までの時間ステップを表します。
   • h は、時間ステップであり、t 配列の最初の 2 要素の差として計算されます。
2. 位置と速度の初期化：
   • 配列 r と v は、各時間ステップで地球の位置と速度を格納するために初期化されます。それらはゼロで初期化され、各時間ステップに対して 2 つの座標（x と y）を収容する形状を持っています。
   • 同様に、r_jupiter と v_jupiter は、木星の位置と速度を格納するために初期化されます。
3. 初期条件の設定：
   • 地球の初期位置（r[0]）は、太陽からの距離を天文単位（AU）で正規化した値に設定されます。
   • 木星の初期位置（r_jupiter[0]）は、太陽から 5.2 AU に設定され、太陽系における実際の位置を反映しています。
   • 地球の初期速度（v[0]）は、太陽によって働く重力に基づいて計算されます。
   • 木星の初期速度（v_jupiter[0]）は、その実際の軌道速度を反映する特定の値に設定されます。

このステップは、シミュレーションの起点を確立するために基本的です。初期位置と速度は、重力による両方の惑星の後続の運動を計算するために不可欠です。

シミュレーションの実行

定義された時間枠でシミュレーションを実行し、RK4 ソルバを使用して惑星の位置を更新します。

## シミュレーションの実行
for i in trange(N - 1, desc="Generating Animation"):
    r[i + 1], v[i + 1] = RK4Solver(t[i], r[i],
                                   v[i], h, 'earth', r_jupiter[i], v_jupiter[i])
    r_jupiter[i + 1], v_jupiter[i +
                                1] = RK4Solver(t[i], r_jupiter[i], v_jupiter[i], h, 'jupiter', r[i], v[i])

このステップでは、重力シミュレーションを実行するためのメインループを実行します。ここでは、RK4Solver 関数を使用して、指定された期間にわたって地球と木星の運動が計算されます。このプロセスの簡単な説明は以下の通りです。

1. シミュレーションループの実行：
   • シミュレーションの時間ステップの総数である N の N - 1 回反復するループが設定されます。このループは、各時間ステップでシミュレーションを進めるために不可欠です。
   • tqdm モジュールの trange 関数を使用して時間ステップを反復します。これは、シミュレーションの進行状況を示すプログレスバー（「アニメーション生成中」と表示されます）を提供します。
2. 地球と木星の位置と速度の更新：
   • ループの各反復では、RK4Solver 関数が 2 回呼び出されます。
     • 1 回目は、次の時間ステップで地球の位置（r[i + 1]）と速度（v[i + 1]）を更新するために呼び出されます。
     • 2 回目は、次の時間ステップで木星の位置（r_jupiter[i + 1]）と速度（v_jupiter[i + 1]）を更新するために呼び出されます。
   • RK4Solver 関数は、惑星の現在の位置と速度、および他の惑星の位置と速度を取り、新しい状態を計算します。
3. 相互作用のシミュレーション：
   • このループは、本質的に、各小さな時間ステップで地球と木星の間の重力相互作用（お互いへの影響と太陽の影響を考慮して）をシミュレートし、時間の経過に伴う彼らの軌道の現実的な表現をもたらします。

このループが終了すると、r、v、r_jupiter、および v_jupiter の配列には、各時間ステップで地球と木星の位置と速度が格納され、プロジェクトの次の段階でアニメーション化する準備が整います。このステップでは、天体力学をシミュレートするために、核心となる物理学と数値方法が適用されます。

アニメーションの表示

最後に、アニメーションを実行して結果を表示します。

## アニメーションのセットアップ
fig, ax, line_earth, line_jupiter, ttl = setup_animation()
## 目盛りとラベルの追加
ax.plot([-6,-5],[6.5,6.5],'r-')
ax.text(-4.5,6.3,r'1 AU = $1.496 \times 10^8$ km')

ax.plot(-6,-6.2,'o', color = '#d2eeff', markerfacecolor = '#0077BE')
ax.text(-5.5,-6.4,'地球')

ax.plot(-3.3,-6.2,'o', color = '#e3dccb',markersize = 8, markerfacecolor = '#f66338')
ax.text(-2.9,-6.4,'超木星 (質量 500 倍)')

ax.plot(5,-6.2,'o', markersize = 9, markerfacecolor = "#FDB813",markeredgecolor ="#FD7813")
ax.text(5.5,-6.4,'太陽')

## アニメーションの作成
anim = animation.FuncAnimation(
    fig, animate, frames=4000, interval=1, blit=False)

## アニメーションの表示
plt.show()

プロジェクトの最後のステップでは、地球と架空の「超木星」の運動を示す重力シミュレーションのアニメーションを表示します。このステップでは、シミュレーション期間中の惑星の計算された位置を可視化することが含まれます。

これですべてのステップが完了したので、次のコマンドを使用してデスクトップ環境でコードを実行できます。

cd ~/project
python simulation.py

このステップでは、大量の「超木星」が地球の軌道にどのような影響を与えるかを観察しながら、シミュレーションの結果を視覚的に鑑賞することができます。

まとめ

このプロジェクトでは、Python を使って地球と木星を対象とした重力シミュレーションを効果的に開発しました。高精度な数値積分のために 4 次ルンゲクッタ法（RK4）を採用し、明瞭で魅力的な可視化のために Matplotlib を使用しました。このプロジェクトは、物理とコーディングの原理を実践的に適用するだけでなく、様々な質量の天体をシミュレートするためにコードを変更することで学生に実験を促し、さまざまな重力相互作用と力学を探求することができます。