import csv
import numpy as np
from math import exp, sqrt, sin, cos, pi
"""
惑星シミュレーター: オイラー法による連立2階常微分方程式の解法
このモジュールは、複数の惑星の運動をシミュレーションするために、
ニュートン重力とオイラー法を用いて2階常微分方程式を解きます。
初期条件は`planet_db.csv`ファイルから読み込まれます。
結果はCSVファイルに出力されます。
:doc:`diffeq2nd_planet_euler_usage`
"""
# d2x/dt2 = force(t, x, y), d2y/dt2 = force(t, x, y)
# define function to be integrated
[ドキュメント]
def force(t, x, y):
"""
質点に作用する力を計算します。
この関数は、与えられた位置(x, y)における質点に作用する力を計算し、
加速度(-x, -y)として返します。これは、単純な調和振動子のようなモデルです。
シミュレーションの主要な重力計算とは異なる、別の簡単なモデルやテストケースで
使用される可能性があります。
:param t: 現在の時刻 (float)。
:param x: 質点のX座標 (float)。
:param y: 質点のY座標 (float)。
:returns: X方向の力 (float), Y方向の力 (float)。
"""
return -x, -y
# solution: x(0.0) = 0.0, y(0.0) = 2.0, vx(0.0) = 1.0, vy(0.0) = 0.0
[ドキュメント]
def xysolution(t):
"""
特定の初期条件における解析解のXおよびY座標を計算します。
この関数は、初期条件 x(0.0)=0.0, y(0.0)=2.0, vx(0.0)=1.0, vy(0.0)=0.0
に対応する調和振動子の解析解 (sin(t), 2.0 * cos(t)) を提供します。
シミュレーション結果と比較するための参照解として利用されることがあります。
:param t: 現在の時刻 (float)。
:returns: X座標 (float), Y座標 (float)。
"""
return sin(t), 2.0 * cos(t)
#===================
# constants
#===================
G = 6.67259e-11 #Nm2/kg2
DayToSecond = 60 * 60 * 24 #s
SecondToDay = 1.0 / DayToSecond
AstronomicalUnit = 1.49597870e11 #m
AU = AstronomicalUnit
G1 = G * DayToSecond * DayToSecond / AU / AU / AU
#===================
# parameters
#===================
dbfile = 'planet_db.csv'
outfile = 'diffeq2nd_Planet_Euler.csv'
x0, vx0 = 0.0, 1.0
y0, vy0 = 2.0, 0.0
z0, vz0 = 0.0, 0.0
dt = 1.0
nt = 1000
iprint_interval = 100
nprint_planets = 4
CSVOutputInterval = 10
CriteriaR2ForIdenticalAtom = 0.001
DifferentialEquationSolver = "Verlet"
#===================
# functions
#===================
[ドキュメント]
def readdb(dbfile):
"""
CSVデータベースファイルから惑星データを読み込みます。
指定されたCSVファイルから惑星の情報を辞書のリストとして読み込みます。
'Name'以外のすべてのキーの値は浮動小数点数に変換されます。
:param dbfile: 惑星データが格納されたCSVファイルのパス (str)。
:returns: 惑星データのリスト (list of dict)。各辞書は惑星の名前、質量、半径、速度などの情報を含みます。
"""
planets = []
f = open(dbfile, "r");
reader = csv.DictReader(f)
for row in reader:
planets.append(row)
keys = list(planets[0].keys())
for d in planets:
for key in keys:
if key != 'Name':
d[key] = float(d[key])
return planets
[ドキュメント]
def sum(array):
"""
数値のリストまたはタプルの合計を計算します。
与えられた数値のイテラブルオブジェクト内の全要素の合計を計算して返します。
:param array: 合計する数値のリストまたはタプル (list or tuple of float/int)。
:returns: 配列内の全要素の合計 (float)。
"""
sum = 0.0
for e in array:
sum += e
return sum
[ドキュメント]
def total_momentum(it, M, vx, vy, vz):
"""
システムの総運動量を計算します。
指定された時刻ステップにおける全惑星の総運動量ベクトル(Px, Py, Pz)を計算します。
:param it: 現在の時刻ステップのインデックス (int)。
:param M: 各惑星の質量を含むリスト (list of float)。
:param vx: 各惑星のX方向速度履歴のリスト (list of list of float)。
:param vy: 各惑星のY方向速度履歴のリスト (list of list of float)。
:param vz: 各惑星のZ方向速度履歴のリスト (list of list of float)。
:returns: X方向総運動量 (float), Y方向総運動量 (float), Z方向総運動量 (float)。
"""
Px, Py, Pz = 0.0, 0.0, 0.0
for i in range(0, len(M)):
Px += M[i] * vx[i][it];
Py += M[i] * vy[i][it];
Pz += M[i] * vz[i][it];
return Px, Py, Pz
[ドキュメント]
def normalize_momentum(it, M, x, y, z, vx, vy, vz, fx, fy, fz):
"""
システムの総運動量をゼロに正規化します。
全惑星の総運動量を計算し、各惑星の速度を調整することで
システムの重心が静止するように総運動量をゼロに設定します。
これにより、シミュレーション開始時のドリフトを防ぎます。
:param it: 現在の時刻ステップのインデックス (int)。
:param M: 各惑星の質量を含むリスト (list of float)。
:param x: 各惑星のX座標履歴のリスト (list of list of float)。
:param y: 各惑星のY座標履歴のリスト (list of list of float)。
:param z: 各惑星のZ座標履歴のリスト (list of list of float)。
:param vx: 各惑星のX方向速度履歴のリスト (list of list of float)。
:param vy: 各惑星のY方向速度履歴のリスト (list of list of float)。
:param vz: 各惑星のZ方向速度履歴のリスト (list of list of float)。
:param fx: 各惑星のX方向力履歴のリスト (list of list of float)。
:param fy: 各惑星のY方向力履歴のリスト (list of list of float)。
:param fz: 各惑星のZ方向力履歴のリスト (list of list of float)。
:returns: 運動量正規化前のX方向総運動量 (float), Y方向総運動量 (float), Z方向総運動量 (float)。
"""
Mtot = sum(M)
Px, Py, Pz = total_momentum(it, M, vx, vy, vz);
for ip in range(0, len(M)):
vx[ip][it] -= Px / Mtot;
vy[ip][it] -= Py / Mtot;
vz[ip][it] -= Pz / Mtot;
return Px, Py, Pz
[ドキュメント]
def initialize(planets, M, x, y, z, vx, vy, vz, fx, fy, fz):
"""
惑星シミュレーションの初期状態を設定します。
読み込んだ惑星データに基づき、質量、初期位置、初期速度、
および初期時点での各惑星にかかる力を設定します。
位置と速度はAUと日単位に正規化されます。
各物理量の履歴リストに最初のデータが追加されます。
:param planets: `readdb`関数から読み込まれた惑星データのリスト (list of dict)。
:param M: 質量を格納する空のリスト (list of float)。この関数内でデータが追加されます。
:param x: X座標履歴を格納する空のリスト (list of list of float)。この関数内で初期位置が追加されます。
:param y: Y座標履歴を格納する空のリスト (list of list of float)。この関数内で初期位置が追加されます。
:param z: Z座標履歴を格納する空のリスト (list of list of float)。この関数内で初期位置が追加されます。
:param vx: X方向速度履歴を格納する空のリスト (list of list of float)。この関数内で初期速度が追加されます。
:param vy: Y方向速度履歴を格納する空のリスト (list of list of float)。この関数内で初期速度が追加されます。
:param vz: Z方向速度履歴を格納する空のリスト (list of list of float)。この関数内で初期速度が追加されます。
:param fx: X方向力履歴を格納する空のリスト (list of list of float)。この関数内で初期力が追加されます。
:param fy: Y方向力履歴を格納する空のリスト (list of list of float)。この関数内で初期力が追加されます。
:param fz: Z方向力履歴を格納する空のリスト (list of list of float)。この関数内で初期力が追加されます。
:returns: None
"""
global AU
global DayToSecond
for i in range(0, len(planets)):
M.append(planets[i]['Mass'])
x.append([planets[i]['Revolution Radius'] / AU])
y.append([0.0])
z.append([0.0])
vx.append([0.0])
vy.append([planets[i]['Revolution Velocity'] * DayToSecond / AU])
vz.append([0.0])
for i in range(0, len(planets)):
fxi, fyi, fzi = Fi(0, i, M, x, y, z)
fx.append([fxi])
fy.append([fyi])
fz.append([fzi])
[ドキュメント]
def Fij(istep, i, j, M, x, y, z):
"""
2つの惑星間に作用する重力ベクトルを計算します。
惑星 `j` が惑星 `i` に及ぼす重力を計算します。
ニュートンの万有引力の法則に基づき、X, Y, Z方向の力を返します。
力は正規化された単位 (AU, Day) で計算されます。
:param istep: 現在の時刻ステップのインデックス (int)。
:param i: 力を受ける惑星のインデックス (int)。
:param j: 力を及ぼす惑星のインデックス (int)。
:param M: 各惑星の質量を含むリスト (list of float)。
:param x: 各惑星のX座標履歴のリスト (list of list of float)。
:param y: 各惑星のY座標履歴のリスト (list of list of float)。
:param z: 各惑星のZ座標履歴のリスト (list of list of float)。
:returns: X方向の力 (float), Y方向の力 (float), Z方向の力 (float)。
"""
dx = x[j][istep] - x[i][istep]
dy = y[j][istep] - y[i][istep]
dz = z[j][istep] - z[i][istep]
r2 = dx*dx + dy*dy + dz*dz
r = sqrt(r2)
g = G1 * M[j]
f = g / r2
fx = f * dx / r
fy = f * dy / r
fz = f * dz / r
return fx, fy, fz
[ドキュメント]
def Fi(istep, i, M, x, y, z):
"""
単一の惑星に作用する全ての重力ベクトルを合計します。
指定された惑星 `i` に、他の全ての惑星から作用する重力の影響を合計し、
その惑星にかかる合力を計算します。
:param istep: 現在の時刻ステップのインデックス (int)。
:param i: 合力を計算する惑星のインデックス (int)。
:param M: 各惑星の質量を含むリスト (list of float)。
:param x: 各惑星のX座標履歴のリスト (list of list of float)。
:param y: 各惑星のY座標履歴のリスト (list of list of float)。
:param z: 各惑星のZ座標履歴のリスト (list of list of float)。
:returns: 惑星 `i` に作用するX方向の合力 (float), Y方向の合力 (float), Z方向の合力 (float)。
"""
fxi = 0.0
fyi = 0.0
fzi = 0.0
for j in range(0, len(M)):
if i == j:
continue
fxj, fyj, fzj = Fij(istep, i, j, M, x, y, z)
fxi += fxj
fyi += fyj
fzi += fzj
# print("f={}, {}, {}".format(fxi, fyi, fzi))
return fxi, fyi, fzi
#===================
# main routine
#===================
[ドキュメント]
def main():
"""
惑星シミュレーションのメインルーチンです。
惑星データを読み込み、初期化を行い、オイラー法を用いて惑星の軌道をシミュレートします。
シミュレーション結果はコンソールに表示され、CSVファイルに保存されます。
"""
global nt
global dt
print("Planet simulator: Solve simulataneous second order diffrential equations by Euler method")
print("G = {} Nm2/kg2".format(G))
print("AU = {:e} m".format(AU))
print("G1 = {}".format(G1))
print("")
print("Planets:")
planets = readdb(dbfile)
keys = list(planets[0].keys())
for d in planets:
print(" ", d['Name'])
for key in keys:
if key != 'Name':
print(" {}: {}".format(key, d[key]))
print("")
M = []
x = []
y = []
z = []
vx = []
vy = []
vz = []
fx = []
fy = []
fz = []
initialize(planets, M, x, y, z, vx, vy, vz, fx, fy, fz)
Px, Py, Pz = normalize_momentum(0, M, x, y, z, vx, vy, vz, fx, fy, fz)
print("")
# make label list for display / csv output
labellist = ['t']
for i in range(0, len(planets)):
labellist.append("x({})".format(planets[i]['Name']))
labellist.append("y({})".format(planets[i]['Name']))
# open outfile to write a csv file
print("Write to [{}]".format(outfile))
f = open(outfile, 'w')
fout = csv.writer(f, lineterminator='\n')
fout.writerow(labellist)
print("{:^5}".format('t'), end = '')
for i in range(1, nprint_planets*2, 2):
print(" {:^12} {:^12}".format(labellist[i], labellist[i+1]), end = '')
print("")
# Solve the 1st data by Heun method
# NOTE: The comment says "Heun method", but the code implements Euler method.
datalist = [0.0]
print("{:^5}".format(0.0), end = '')
for i in range(0, len(planets)):
fx0, fy0, fz0 = Fi(0, i, M, x, y, z)
vx1 = vx[i][0] + dt * fx0
vy1 = vy[i][0] + dt * fy0
vz1 = vz[i][0] + dt * fz0
x1 = x[i][0] + dt * vx[i][0]
y1 = y[i][0] + dt * vy[i][0]
z1 = z[i][0] + dt * vz[i][0]
datalist.append(x[i][0])
datalist.append(y[i][0])
x[i].append(x1)
y[i].append(y1)
z[i].append(z1)
vx[i].append(vx1)
vy[i].append(vy1)
vz[i].append(vz1)
for i in range(1, nprint_planets*2, 2):
print(" {:^12.4f} {:^12.4f}".format(x[i][0], y[i][0]), end = '')
print("")
fout.writerow(datalist)
for it in range(1, nt):
t = it * dt
# print("it={} t={}".format(it, t))
datalist = [t]
if it % iprint_interval == 0:
print("{:^5}".format(t), end = '')
for i in range(0, len(planets)):
fx0, fy0, fz0 = Fi(it, i, M, x, y, z)
vx1 = vx[i][it] + dt * fx0
vy1 = vy[i][it] + dt * fy0
vz1 = vz[i][it] + dt * fz0
x1 = x[i][it] + dt * vx[i][it]
y1 = y[i][it] + dt * vy[i][it]
z1 = z[i][it] + dt * vz[i][it]
datalist.append(x[i][it])
datalist.append(y[i][it])
x[i].append(x1)
y[i].append(y1)
z[i].append(z1)
vx[i].append(vx1)
vy[i].append(vy1)
vz[i].append(vz1)
if it % iprint_interval == 0:
for i in range(1, nprint_planets*2, 2):
print(" {:^12.4g} {:^12.4g}".format(x[i][it], y[i][it]), end = '')
print("")
fout.writerow(datalist)
f.close()
if __name__ == '__main__':
main()