import csv
import numpy as np
from math import exp, sqrt, sin, cos, pi
"""
惑星シミュレーター: Verlet法による連立2階常微分方程式の解法
このスクリプトは、Verlet法を用いて複数の天体間の重力相互作用に基づく運動をシミュレーションします。
天体の位置と速度の時間発展を計算し、結果をCSVファイルに出力します。
関連リンク: :doc:`diffeq2nd_planet_verlet_usage`
"""
# d2x/dt2 = force(t, x, y), d2y/dt2 = force(t, x, y)
# define function to be integrated
[ドキュメント]
def force(t, x, y):
"""
特定の時間と位置におけるテスト用の力を計算します。
この関数は、Verlet法のデモンストレーションのために用意された、シンプルな調和振動子のような力を返します。
実際の惑星シミュレーションでは使用されません。
:param t: float, 現在の時間。
:param x: float, x座標。
:param y: float, y座標。
:returns: tuple[float, float], x方向とy方向の力成分 (fx, fy)。
"""
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):
"""
特定の時間におけるテスト用の解析解を計算します。
この関数は、Verlet法のデモンストレーションのために用意された、調和振動子問題の解析解を返します。
実際の惑星シミュレーションでは使用されません。
:param t: float, 現在の時間。
:returns: tuple[float, float], 解析解のx座標とy座標 (x_exact, y_exact)。
"""
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_Verlet.csv'
x0, vx0 = 0.0, 1.0
y0, vy0 = 2.0, 0.0
dt = 0.01
nt = 501
iprint_interval = 20
nprint_planets = 3
#===================
# functions
#===================
[ドキュメント]
def Fij(istep, i, j, M, x, y, z):
"""
2つの天体間に働く重力相互作用を計算します。
指定されたステップにおける天体iと天体jの間の重力（引力）の各成分を計算します。
ニュートンの万有引力の法則に基づいています。
:param istep: int, 現在のシミュレーションステップ。
:param i: int, 最初の天体のインデックス。
:param j: int, 2番目の天体のインデックス。
:param M: list[float], 各天体の質量リスト。
:param x: list[list[float]], 各天体のx座標履歴リスト。
:param y: list[list[float]], 各天体のy座標履歴リスト。
:param z: list[list[float]], 各天体のz座標履歴リスト。
:returns: tuple[float, float, float], 天体iとjの間に働く力のx, y, z成分 (fx, fy, fz)。
"""
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[i] * 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[float], 各天体の質量リスト。
:param x: list[list[float]], 各天体のx座標履歴リスト。
:param y: list[list[float]], 各天体のy座標履歴リスト。
:param z: list[list[float]], 各天体のz座標履歴リスト。
:returns: tuple[float, float, float], 天体iに作用する全重力のx, y, z成分 (fxi, fyi, fzi)。
"""
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
return fxi, fyi, fzi
[ドキュメント]
def initialize(planets, M, x, y, z, vx, vy, vz, fx, fy, fz):
"""
惑星の初期状態と初期力を設定します。
与えられた惑星データに基づいて、質量、初期位置、初期速度、
そして最初のステップにおける各天体に作用する力を初期化します。
座標、速度、力の各リストは、この関数内で追加されます。
:param planets: list[dict], 惑星のデータを含む辞書のリスト。
:param M: list[float], 各天体の質量を格納するリスト (この関数内で更新)。
:param x: list[list[float]], 各天体のx座標履歴を格納するリスト (この関数内で更新)。
:param y: list[list[float]], 各天体のy座標履歴を格納するリスト (この関数内で更新)。
:param z: list[list[float]], 各天体のz座標履歴を格納するリスト (この関数内で更新)。
:param vx: list[list[float]], 各天体のx方向速度履歴を格納するリスト (この関数内で更新)。
:param vy: list[list[float]], 各天体のy方向速度履歴を格納するリスト (この関数内で更新)。
:param vz: list[list[float]], 各天体のz方向速度履歴を格納するリスト (この関数内で更新)。
:param fx: list[list[float]], 各天体のx方向の力を格納するリスト (この関数内で更新)。
:param fy: list[list[float]], 各天体のy方向の力を格納するリスト (この関数内で更新)。
:param fz: list[list[float]], 各天体のz方向の力を格納するリスト (この関数内で更新)。
:returns: None
"""
for i in range(0, len(planets)):
M.append(planets[i]['Mass'])
x.append([planets[i]['Revolution Radius']])
y.append([0.0])
z.append([0.0])
vx.append([0.0])
vy.append([planets[i]['Revolution Velocity']])
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 readdb(dbfile):
"""
惑星データベースファイルから惑星データを読み込みます。
指定されたCSVファイルから惑星の名前、質量、初期位置、初期速度などのデータを読み込み、
数値データをfloat型に変換してリストとして返します。
:param dbfile: str, 惑星データベースのCSVファイルパス。
:returns: list[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
#===================
# main routine
#===================
[ドキュメント]
def main():
"""
惑星シミュレーターのメインルーチン。
惑星データベースの読み込み、初期化、そしてVerlet法による惑星の軌道計算（現在コメントアウトされています）を
実行し、結果をCSVファイルに出力します。
:returns: None
"""
global nt
global dt
print("Planet simulator: Solve simulataneous second order diffrential equations by Verlet method")
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)
# 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):
name = planets[i]['Name']
print(" {:^12} {:^12}".format(labellist[i], labellist[i+1]), end = '')
print("")
"""
# Solve the 1st data by Heun method
kvx0, kvy0 = force(0.0, x0, y0)
kvx1, kvy1 = force(0.0+dt, x0+kvy0, y0+kvy0)
vx1 = vx0 + dt * (kvx0 + kvx1) / 2.0
vy1 = vy0 + dt * (kvy0 + kvy1) / 2.0
kx0 = dt * vx0
kx1 = dt * vx1
ky0 = dt * vy0
ky1 = dt * vy1
x1 = x0 + (kx0 + kx1) / 2.0
y1 = y0 + (ky0 + ky1) / 2.0
xexact, yexact = xysolution(0.0)
print("t={:5.2f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}".format(0.0, x0, xexact, y0, yexact))
fout.writerow([0.0, x0, xexact, y0, yexact])
for i in range(1, nt):
t1 = i * dt
fx1, fy1 = force(t1, x1, y1)
x2 = 2.0 * x1 - x0 + dt * dt * fx1
y2 = 2.0 * y1 - y0 + dt * dt * fy1
vx1 = (x2 - x0) / 2.0 / dt
vy1 = (y2 - y0) / 2.0 / dt
xexact, yexact = xysolution(t1)
if i == 1 or i % iprint_interval == 0:
print("t={:5.2f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}".format(t1, x1, xexact, y1, yexact))
fout.writerow([t1, x1, xexact, y1, yexact])
x0 = x1
vx0 = vx1
y0 = y1
vy0 = vy1
x1 = x2
y1 = y2
f.close()
"""
if __name__ == '__main__':
main()