惑星の運動を差分で求め、アニメーションにするスクリプト

地球と木星の運動を差分で求め、アニメーションにします。
chmod +x planets.pyで実行権を付け、./planets.pyで動作させてください。
まず、惑星の軌道をプロットした図が出ます。それを「X」で閉じると、アニメーションファイルを出力します。

#!/usr/bin/env python3

# シンプルな惑星運動の差分（シンプレクティック・リープフロッグ）シミュレーション
# - 単位: 距離 AU, 質量 太陽質量, 時間 年
# - G = 4*pi^2 (AU^3 / (yr^2 * Msun))
# - 例: 太陽・地球・木星の3体シミュレーション
# 実行すると軌道をプロットします。必要ならパラメータを編集してください。

import numpy as np
import matplotlib.pyplot as plt
from math import pi
import matplotlib.animation as animation
import os

# --- パラメータ ---
G = 4 * pi**2  # AU^3 / (yr^2 * Msun)
dt = 0.1    # 年単位のタイムステップ (~0.18日)
years = 12     # シミュレーション年数
steps = int(years / dt)

# 3体の初期条件 (質量, x,y,z, vx,vy,vz)
# 太陽, 地球, 木星 (データは概算)
bodies = [
    {"name":"Sun",   "m":1.0,        "r":np.array([0.0,0.0,0.0]),          "v":np.array([0.0,0.0,0.0])},
    {"name":"Earth", "m":3.003e-6,   "r":np.array([1.0,0.0,0.0]),          "v":np.array([0.0,2*pi,0.0])},  # circular approx
    {"name":"Jupiter","m":0.0009543, "r":np.array([5.2,0.0,0.0]),          "v":np.array([0.0,2*pi/np.sqrt(5.2),0.0])}
]

N = len(bodies)
masses = np.array([b["m"] for b in bodies])
r = np.array([b["r"].astype(float) for b in bodies])
v = np.array([b["v"].astype(float) for b in bodies])

# 保存用配列
traj = np.zeros((steps, N, 3))
time = np.zeros(steps)

def accelerations(positions):
    a = np.zeros_like(positions)
    for i in range(N):
        for j in range(N):
            if i==j: continue
            diff = positions[j] - positions[i]
            dist3 = np.linalg.norm(diff)**3 + 1e-12
            a[i] += G * masses[j] * diff / dist3
    return a

# 初期加速度
a = accelerations(r)

# シンプレクティック・リープフロッグ (velocity Verlet 風)
for k in range(steps):
    # 半ステップで速度更新
    v_half = v + 0.5 * dt * a
    # フルステップで位置更新
    r = r + dt * v_half
    # 新しい加速度
    a_new = accelerations(r)
    # 完全な速度更新
    v = v_half + 0.5 * dt * a_new
    a = a_new
    traj[k] = r
    time[k] = k * dt

# プロット
plt.figure(figsize=(8,8))
for i in range(N):
    plt.plot(traj[:,i,0], traj[:,i,1], label=bodies[i]["name"])
    plt.scatter(traj[-1,i,0], traj[-1,i,1], s=20)  # 最終位置マーカー
plt.gca().set_aspect('equal', 'box')
plt.xlabel("x (AU)")
plt.ylabel("y (AU)")
plt.title(f"3-body simulation (leapfrog), dt={dt} yr, T={years} yr")
plt.legend()
plt.grid(True)
# 画像保存
outpath = "planetary_orbits.png"
plt.savefig(outpath, dpi=150, bbox_inches='tight')
plt.show()
print(f"画像を保存しました: {outpath}")

# 小さなエネルギーチェック（運動+ポテンシャル）
def total_energy(positions, velocities):
    K = 0.5 * np.sum(masses[:,None] * velocities**2)
    U = 0.0
    for i in range(N):
        for j in range(i+1,N):
            r_ij = np.linalg.norm(positions[i]-positions[j])
            U -= G * masses[i] * masses[j] / r_ij
    return K + U

E_initial = total_energy(np.array([b["r"] for b in bodies]), np.array([b["v"] for b in bodies]))
E_final = total_energy(traj[-1], v)
print(f"初期エネルギー: {E_initial:.6e}, 最終エネルギー: {E_final:.6e}, 変化: {((E_final-E_initial)/abs(E_initial))*100:.3f}%")


# === アニメーション設定 ===
fig, ax = plt.subplots(figsize=(7,7))
ax.set_aspect('equal', 'box')
ax.set_xlim(-6, 6)
ax.set_ylim(-6, 6)
ax.set_xlabel("x (AU)")
ax.set_ylabel("y (AU)")
ax.set_title("Planetary Motion Simulation")

lines = []
points = []
for i in range(N):
    line, = ax.plot([], [], '-', lw=1, label=bodies[i]["name"])
    point, = ax.plot([], [], 'o', markersize=6)
    lines.append(line)
    points.append(point)

ax.legend()

# === 更新関数 ===
def init():
    for line, point in zip(lines, points):
        line.set_data([], [])
        point.set_data([], [])
    return lines + points

def update(frame):
    for i, (line, point) in enumerate(zip(lines, points)):
        line.set_data(traj[:frame, i, 0], traj[:frame, i, 1])
        point.set_data([traj[frame, i, 0]], [traj[frame, i, 1]])
    return lines + points

# === アニメーション生成 ===
frames = len(traj)
print("Wait... now making animation.")
ani = animation.FuncAnimation(fig, update, frames=frames, init_func=init, interval=10, blit=True)

# 保存（MP4ファイル）
out_movie = "planetary_motion.mp4"
ani.save(out_movie, writer="ffmpeg", fps=60)

plt.close(fig)
print(f"アニメーションを保存しました: {out_movie}")