初めての投稿です。
はじめに
2018年7月31日、火星が大接近する/したというニュース。
そのなかで、
地球と火星はおよそ２年２か月ごとに接近しますが、地球は太陽を中心に円軌道を回っているのに対し、火星はだ円形の軌道を回っているため、接近しても距離が近い時と遠く離れる時があります。 NHKニュース
ただ、テレビで見たときはこの中の一文、
地球は太陽を中心に円軌道を回っているのに対し、火星はだ円形の軌道を回っているため
何年かおきに接近する、みたいな言い回しじゃなかったかな……(真剣に聞いていたわけではないから今となっては自信がない)。
これを聞いて、「地球は円で、火星は楕円って本当?」「火星の接近って、公転周期が違うのが効いているんじゃないの?」「公転面の傾きが効いているかも」と疑問がでたので、とりあえずWikipediaで調べてみた。地球, 火星
これによると、地球軌道の離心率0.01671022に対して、火星軌道の離心率0.09341233なので、ざっくり一桁違う。かといって、離心率0.1の火星軌道を楕円というなら地球も楕円と言ってもいいのでは…?よくわからん。
ただ、火星の観測 に今回の疑問はたいてい答えてあって、
地球は780日（2年と7週間と1日）ごとに火星を追い越し、そのときの距離は約8000万km（約4光分）まで接近する。しかし、火星軌道が楕円であるために最接近時の距離は変化する。火星の近日点付近で接近すれば接近距離は5600万km程度となるが、遠日点付近で接近すれば1億km程度と2倍近く距離が異なる。(中略)最も観測に適した時期は32年ごとに2回、15年と17年をおいて交互にやってきて「大接近」と呼ばれる。
要するに、火星との距離は主に公転周期が違うのが効いてきて、次に火星の軌道が楕円であることが効いているということらしい。だいたいのイメージはわかったけど、さすがに公転面の傾きについては記載がなかったので、シミュレーションしてみようか、ということになった。
ソースコード
基本的には運動方程式を解くだけなので、細かい説明はいらないと思う。
下準備
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
c_cycle=("#ea5415","#005192","#ffdf00","#1d7a21","#88593b","#737061")
class vec3d:
def __init__(self, x,y,z):
self.x = x
self.y = y
self.z = z
class planet:
def __init__(self, position, mass, velocity, name = ""):
self.position = position
self.mass = mass
self.velocity = velocity
self.name = name
うまく使いこなせていないけどnumpy,描画のためのmatplotlib,3次元なのでAxes3Dをインポート。
3次元ベクトルと惑星くらいはクラスとして定義しておく（関数を持っていないのでクラスというにはおこがましいが）。
c_cycleはなんとなく使っているミッフィーちゃんのカラーリング。
#planet data (position (m), mass (kg), velocity (m/s)
sun = {"position":vec3d(0,0,0), "mass":1.9891e30, "velocity":vec3d(0,0,0)}
mercury = {"position":vec3d(0,5.791e10,0), "mass":3.301e23, "velocity":vec3d(478725,0,0)}
venus = {"position":vec3d(0,1.0820893e11,0), "mass":4.869e24, "velocity":vec3d(35021.4,0,0)}
earth = {"position":vec3d(0,1.4959787e11,0), "mass":5.972e24, "velocity":vec3d(29780,0,0)}
mars = {"position":vec3d(0,2.2793664e11,0), "mass":6.4191e23, "velocity":vec3d(24131,0,0)}
mars2 = {"position":vec3d(0,2.2793664e11,0), "mass":6.4191e23, "velocity":vec3d(24131*np.cos(1.8497*np.pi/180.0),0,24131*np.sin(1.8497*np.pi/180.0))}
jupiter = {"position":vec3d(0,7.78412e11,0), "mass":1.8986e27, "velocity":vec3d(13070,0,0)}
saturn = {"position":vec3d(0,1.4267254e12,0), "mass":5.688e26, "velocity":vec3d(9672,0,0)}
uranus = {"position":vec3d(0,2.87099e12,0), "mass":8.686e25, "velocity":vec3d(6835,0,0)}
neptune = {"position":vec3d(0,4.49506e12,0), "mass":1.02e26, "velocity":vec3d(5477,0,0)}
pluto = {"position":vec3d(0,3.7e12,0), "mass":1.3e22, "velocity":vec3d(4748,0,0)}
solar_system = [
planet( position = sun["position"], mass = sun["mass"], velocity = sun["velocity"], name = "sun"),
planet( position = earth["position"], mass = earth["mass"], velocity = earth["velocity"], name = "earth"),
planet( position = mars["position"], mass = mars["mass"], velocity = mars["velocity"], name = "mars"),
planet( position = venus["position"], mass = venus["mass"], velocity = venus["velocity"], name = "venus"),
planet( position = jupiter["position"], mass = jupiter["mass"], velocity = jupiter["velocity"], name = "jupiter"),
planet( position = saturn["position"], mass = saturn["mass"], velocity = saturn["velocity"], name = "saturn"),
]
Wikipediaからとってきた惑星情報（初期座標、質量、初期速度）をつかって、solar_systemに登録する。
marsは公転面の傾きがないもの、mars2は公転面が現実と同じだけ傾いているもの。marsをsolar_sytemに登録した場合と、mars2をsolar_systemに登録した場合で地球との距離を比較する。
今回のソースコードで一番めんどくさかったのが惑星情報の入力。探してコピーするだけでめんどくさくて、規格化してない。実際にはすべて出来上がった後入力していて、「惑星の絶対座標と速度がないじゃん!!!」となってしまったので、位置は太陽との平均距離、速度は平均速度にして、xy平面上のy方向に並べている。こんなので本当に検証できるのだろうか。
運動方程式を解く
解き方はいろいろあるけど、お手軽に速度ベルレを採用。
def calculate_acceleration(solar_system, target_index):
G_const = 6.67408e-11 #m3 kg-1 s-2
acceleration = vec3d(0,0,0)
target = solar_system[target_index]
for index, external in enumerate(solar_system):
if index != target_index:
r = (target.position.x - external.position.x)**2 + (target.position.y - external.position.y)**2 + (target.position.z - external.position.z)**2
r = np.sqrt(r)
coef = G_const * external.mass / r**3
acceleration.x += coef * (external.position.x - target.position.x)
acceleration.y += coef * (external.position.y - target.position.y)
acceleration.z += coef * (external.position.z - target.position.z)
return acceleration
def solve_equation_of_motion(solar_system, dt =1):
for index, target in enumerate(solar_system):
acceleration = calculate_acceleration(solar_system, index)
target.position.x += target.velocity.x * dt + 0.5 * acceleration.x * dt**2
target.position.y += target.velocity.y * dt + 0.5 * acceleration.y * dt**2
target.position.z += target.velocity.z * dt + 0.5 * acceleration.z * dt**2
acceleration_new = calculate_acceleration(solar_system, index)
target.velocity.x += 0.5 * (acceleration.x + acceleration_new.x) * dt
target.velocity.y += 0.5 * (acceleration.y + acceleration_new.y) * dt
target.velocity.z += 0.5 * (acceleration.z + acceleration_new.z) * dt
非常に素朴なコードです。
運動方程式を解いて、適当な時間間隔でtrajectoryとして保存しておく。
def simulation(solar_system, dt = 1, number_of_steps = 10000, report_freq = 100):
trajectories = []
for current in solar_system:
trajectories.append({"x":[], "y":[], "z":[], "name":current.name})
for i in range(1,number_of_steps+1):
solve_equation_of_motion(solar_system, dt)
if i % report_freq == 0:
for index, trajectory in enumerate(trajectories):
trajectory["x"].append(solar_system[index].position.x)
trajectory["y"].append(solar_system[index].position.y)
trajectory["z"].append(solar_system[index].position.z)
return trajectories
可視化と解析
def plot_output(solar_system, outfile = None):
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')
max_range = 0
for current in solar_system:
temp = max(max(current["x"]), max(current["y"]), max(current["z"]))
if temp > max_range:
max_range = temp
ax.plot(current["x"], current["y"], current["z"], label = current["name"])
ax.set_xlim([-max_range, max_range])
ax.set_ylim([-max_range, max_range])
ax.set_zlim([-max_range, max_range])
ax.legend(loc='best')
if outfile:
plt.savefig(outfile)
else:
plt.show()
3Dでぐりぐり回したり、拡大縮小できたりすればいいんだろうけど、探すのがめんどくさいのでやめた。可視化は今回の目的ではないので。
def calculate_distance(src_trajectory, dst_trajectory):
distance = []
for t in range(len(src_trajectory[0])):
r = (src_trajectory[0][t] - dst_trajectory[0][t])**2 + (src_trajectory[1][t] - dst_trajectory[1][t])**2 + (src_trajectory[2][t] - dst_trajectory[2][t])**2
distance.append(np.sqrt(r))
return distance
def get_planet_index(solar_system, name):
for index, target in enumerate(solar_system):
if target.name==name:
return index
def plot_distance(trajectories, solar_system, src='earth', dst=None):
src_index = get_planet_index(solar_system, src)
src_trajectory = [ trajectories[src_index]['x'], trajectories[src_index]['y'], trajectories[src_index]['z'] ]
time = np.arange(0, len(src_trajectory[0]))
distance = []
labels = []
if dst==None:
for index, trajectory in enumerate(trajectories):
if index!=src_index:
name = trajectories[index]['name']
pos_trajectory = [ trajectories[index]['x'], trajectories[index]['y'], trajectories[index]['z'] ]
labels.append(src+'-'+name)
d = calculate_distance(pos_trajectory, src_trajectory)
distance.append(d)
else:
dst_index = get_planet_index(solar_system, dst)
pos_trajectory = [ trajectories[dst_index]['x'], trajectories[dst_index]['y'], trajectories[dst_index]['z'] ]
d = calculate_distance(src_trajectory, pos_trajectory)
distance.append(d)
labels.append(src+'-'+dst)
plt.rcParams['figure.figsize'] = 18,9
plt.rcParams["font.size"] = 16
for k in range(len(labels)):
plt.plot(time, distance[k], label=labels[k], color=c_cycle[k])
plt.xlabel('time [day]')
plt.ylabel('distance [km]')
plt.legend(loc='upper left')
return distance
距離の計算とグラフ化。
メイン関数
trajectories = simulation(solar_system, dt = 175.5, number_of_steps = 365*500, report_freq = 500)
d1 = plot_output(trajectories)
dtとreport_freqをうまく調整して、trajectoryの1step分が一日くらいに相当するようにした。20年分くらい計算しても地球がどこかへ飛んでいかなかったので、まあこれくらいでいいのだろう。
結果
動作確認を兼ねた、地球と他の惑星の距離。
Wikipedia見たときから理解していたけど、火星大接近って、0.5e11kmか0.7e11kmなのかという、宇宙スケールでは小さな違いだっていうのがよくわかる。
横軸は時間(単位は日),縦軸は距離(単位はkm)。地球と火星の距離をグラフにした。オレンジ色が公転面が同一、青色が公転面が傾いているもの。
極小値を拾ってグラフにすればいいのだろうが、手で拾う以外のアイデアが思いつかなかったので、これでご勘弁を。
オレンジの極小の変化よりも青の極小の変化が激しいように見えるので、公転面の傾き(僅か1.8497°)はかなり影響が大きいようである。
コメントお待ちしてます