提交 7cd7d202 编写于 作者: M march3

三体运行模拟器

上级 0d59b92d
from bodies.body import Body from bodies.body import Body
from bodies.earth import Earth from bodies.earth import Earth
from bodies.sun import Sun from bodies.sun import Sun
from bodies.jupiter import Jupiter from bodies.jupiter import Jupiter
\ No newline at end of file from bodies.mars import Mars
from bodies.mercury import Mercury
from bodies.neptune import Neptune
from bodies.pluto import Pluto
from bodies.saturn import Saturn
from bodies.uranus import Uranus
from bodies.venus import Venus
...@@ -9,11 +9,7 @@ ...@@ -9,11 +9,7 @@
import json import json
import numpy as np import numpy as np
import math import math
from common.consts import AU
"""
天文单位
"""
AU: float = 149597870.700
class Body: class Body:
...@@ -59,8 +55,6 @@ class Body: ...@@ -59,8 +55,6 @@ class Body:
self.__acceleration = np.array([0, 0, 0], dtype='float32') self.__acceleration = np.array([0, 0, 0], dtype='float32')
self.__record_history() self.__record_history()
@property @property
def position(self): def position(self):
return self.__position return self.__position
......
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
from bodies.body import Body, AU
class Mars(Body):
"""
火星
------------------------
远日点距离: 1.666 天文单位
近日点距离: 1.382 天文单位
逃逸速度: 5.027 km/s
 公转速度: 24.13 km/s
 天体质量: 6.4171✕10²³
 平均密度: 3.9335 g/cm³ -> 3.9335✕10³ kg/m³
"""
def __init__(self, name="Mars", mass=6.4171e23,
init_position=[1.5 * AU, 0, 0],
init_velocity=[0, 24.13, 0],
texture="", size_scale=1.0, distance_scale=1.0):
params = {
"name": name,
"mass": mass,
"init_position": init_position,
"init_velocity": init_velocity,
"density": 3.9335e3,
"color": [
125,
125,
125
],
"texture": texture,
"size_scale": size_scale,
"distance_scale": distance_scale
}
super().__init__(**params)
if __name__ == '__main__':
mars = Mars()
print(mars)
\ No newline at end of file
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
from bodies.body import Body, AU
class Mercury(Body):
"""
水星
------------------------
远日点距离: 0.466697 天文单位
近日点距离: 0.307499 天文单位
逃逸速度: 4.25 km/s
 公转速度: 47.87 km/s
 天体质量: 3.3011✕10²³ kg
 平均密度: 5.427 g/cm³ -> 5.427×10³ kg/m³
"""
def __init__(self, name="Mercury", mass=3.3011e23,
init_position=[0.4 * AU, 0, 0],
init_velocity=[0, 47.87, 0],
texture="", size_scale=1.0, distance_scale=1.0):
params = {
"name": name,
"mass": mass,
"init_position": init_position,
"init_velocity": init_velocity,
"density": 5.427e3,
"color": [
125,
125,
125
],
"texture": texture,
"size_scale": size_scale,
"distance_scale": distance_scale
}
super().__init__(**params)
if __name__ == '__main__':
mercury = Mercury()
print(mercury)
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
from bodies.body import Body, AU
class Neptune(Body):
"""
海王星
------------------------
远日点距离: 30.33 天文单位
近日点距离: 29.81 天文单位
逃逸速度: 23.5 km/s
 公转速度: 5.43 km/s
 天体质量: 1.0241✕10²⁶ kg
 平均密度: 1.638 g/cm³ -> 1.638×10³ kg/m³
"""
def __init__(self, name="Neptune", mass=1.0241e26,
init_position=[30 * AU, 0, 0],
init_velocity=[0, 5.43, 0],
texture="", size_scale=1.0, distance_scale=1.0):
params = {
"name": name,
"mass": mass,
"init_position": init_position,
"init_velocity": init_velocity,
"density": 1.638e3,
"color": [
125,
125,
125
],
"texture": texture,
"size_scale": size_scale,
"distance_scale": distance_scale
}
super().__init__(**params)
if __name__ == '__main__':
neptune = Neptune()
print(neptune)
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
from bodies.body import Body, AU
class Pluto(Body):
"""
冥王星
------------------------
远日点距离: 49.305 天文单位(73.760 亿千米)
近日点距离: 29.658 天文单位(44.368 亿千米)
逃逸速度: 1.212 km/s
 公转速度: 4.7 km/s
 天体质量: 1.303✕10²² kg(±0.003)
 平均密度: 1.854 g/cm³(±0.006) -> 1.854×10³ kg/m³
"""
def __init__(self, name="Pluto", mass=1.303e22,
init_position=[40 * AU, 0, 0],
init_velocity=[0, 4.7, 0],
texture="", size_scale=1.0, distance_scale=1.0):
params = {
"name": name,
"mass": mass,
"init_position": init_position,
"init_velocity": init_velocity,
"density": 1.854e3,
"color": [
125,
125,
125
],
"texture": texture,
"size_scale": size_scale,
"distance_scale": distance_scale
}
super().__init__(**params)
if __name__ == '__main__':
pluto = Pluto()
print(pluto)
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
from bodies.body import Body, AU
class Saturn(Body):
"""
土星
------------------------
远日点距离: 10.1238 天文单位
近日点距离: 9.0412 天文单位
逃逸速度: 35.49 km/s
 公转速度: 9.64 km/s
 天体质量: 5.6834✕10²⁶ kg
 平均密度: 0.687 g/cm³ -> 0.687×10³ kg/m³
"""
def __init__(self, name="Saturn", mass=5.6834e26,
init_position=[10 * AU, 0, 0],
init_velocity=[0, 9.64, 0],
texture="", size_scale=1.0, distance_scale=1.0):
params = {
"name": name,
"mass": mass,
"init_position": init_position,
"init_velocity": init_velocity,
"density": 0.687e3,
"color": [
125,
125,
125
],
"texture": texture,
"size_scale": size_scale,
"distance_scale": distance_scale
}
super().__init__(**params)
if __name__ == '__main__':
saturn = Saturn()
print(saturn)
...@@ -10,6 +10,13 @@ from bodies.body import Body ...@@ -10,6 +10,13 @@ from bodies.body import Body
class Sun(Body): class Sun(Body):
"""
太阳
------------------------
天体质量: 1.9891×10³⁰ kg
平均密度: 1.408×10³ kg/m³
"""
def __init__(self, name="Sun", mass=1.9891e30, def __init__(self, name="Sun", mass=1.9891e30,
init_position=[0, 0, 0], init_position=[0, 0, 0],
init_velocity=[0, 0, 0], init_velocity=[0, 0, 0],
......
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
from bodies.body import Body, AU
class Uranus(Body):
"""
天王星
------------------------
远日点距离: 20.11 天文单位
近日点距离: 18.33 天文单位
 逃逸速度: 21.3 km/s
 公转速度: 6.81 km/s
 天体质量: 8.681✕10²⁵ kg(±0.0013)
 平均密度: 1.27 g/cm³ -> 1.27×10³ kg/m³
"""
def __init__(self, name="Uranus", mass=8.681e25,
init_position=[19 * AU, 0, 0],
init_velocity=[0, 6.81, 0],
texture="", size_scale=1.0, distance_scale=1.0):
params = {
"name": name,
"mass": mass,
"init_position": init_position,
"init_velocity": init_velocity,
"density": 1.27e3,
"color": [
125,
125,
125
],
"texture": texture,
"size_scale": size_scale,
"distance_scale": distance_scale
}
super().__init__(**params)
if __name__ == '__main__':
uranus = Uranus()
print(uranus)
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
from bodies.body import Body, AU
class Venus(Body):
"""
金星
------------------------
远日点距离: 0.728213 天文单位
近日点距离: 0.718440天文单位
逃逸速度: 10.36 km/s
 公转速度: 35 km/s
 天体质量: 4.8675✕10²⁴ kg
 平均密度: 5.24g/cm3 -> 5.24×10³ kg/m³
"""
def __init__(self, name="Venus", mass=4.8675e24,
init_position=[0.72 * AU, 0, 0],
init_velocity=[0, 35, 0],
texture="", size_scale=1.0, distance_scale=1.0):
params = {
"name": name,
"mass": mass,
"init_position": init_position,
"init_velocity": init_velocity,
"density": 5.24e3,
"color": [
125,
125,
125
],
"texture": texture,
"size_scale": size_scale,
"distance_scale": distance_scale
}
super().__init__(**params)
if __name__ == '__main__':
venus = Venus()
print(venus)
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
"""
天文单位
"""
AU: float = 149597870.700
"""
万有引力常数
"""
G: float = 6.67e-11
"""
一小时多少秒
"""
SECONDS_PER_HOUR = 60 * 60
"""
一天多少秒
"""
SECONDS_PER_DAY = SECONDS_PER_HOUR * 24
"""
一周多少秒
"""
SECONDS_PER_WEEK = SECONDS_PER_DAY * 7
"""
一月多少秒(按照30天)
"""
SECONDS_PER_MONTH = SECONDS_PER_DAY * 30
...@@ -6,8 +6,9 @@ ...@@ -6,8 +6,9 @@
# notes : # notes :
# python_version :3.8 # python_version :3.8
# ============================================================================== # ==============================================================================
from common.consts import SECONDS_PER_DAY, SECONDS_PER_HOUR, SECONDS_PER_MONTH, SECONDS_PER_WEEK
from common.system import System from common.system import System
from bodies import Body, Sun, Earth, Jupiter from bodies import Body, Sun, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
...@@ -26,9 +27,9 @@ Y_MIN, Y_MAX = -1e+9, 1e+9 # the y range of the bounding box (m) ...@@ -26,9 +27,9 @@ Y_MIN, Y_MAX = -1e+9, 1e+9 # the y range of the bounding box (m)
Z_MIN, Z_MAX = -1e+9, 1e+9 # the z range of the bounding box (m) Z_MIN, Z_MAX = -1e+9, 1e+9 # the z range of the bounding box (m)
# X_MIN, X_MAX = -8e+8, 8e+8 # the x range of the bounding box (m) X_MIN, X_MAX = -8e+8, 8e+8 # the x range of the bounding box (m)
# Y_MIN, Y_MAX = -8e+8, 8e+8 # the y range of the bounding box (m) Y_MIN, Y_MAX = -8e+8, 8e+8 # the y range of the bounding box (m)
# Z_MIN, Z_MAX = -8e+8, 8e+8 # the z range of the bounding box (m) Z_MIN, Z_MAX = -8e+8, 8e+8 # the z range of the bounding box (m)
def show(bodies, idx=0): def show(bodies, idx=0):
...@@ -59,7 +60,7 @@ def show(bodies, idx=0): ...@@ -59,7 +60,7 @@ def show(bodies, idx=0):
if len(his_pos) > 1: if len(his_pos) > 1:
_his_pos = list(zip(*his_pos)) _his_pos = list(zip(*his_pos))
ax.plot3D(_his_pos[0], _his_pos[1], _his_pos[2], color=color, alpha=0.8) ax.plot3D(_his_pos[0], _his_pos[1], _his_pos[2], color=color, alpha=0.8)
plt.title("3D scatter plot %s" % idx) plt.title("3D %s" % idx)
# display the plot # display the plot
plt.show() plt.show()
...@@ -80,7 +81,7 @@ class Simulator: ...@@ -80,7 +81,7 @@ class Simulator:
# acceleration 加速度 # acceleration 加速度
body.velocity += body.acceleration * dt body.velocity += body.acceleration * dt
# body.position += 0.5 * body.acceleration * (dt ** 2) # body.position += 0.5 * body.acceleration * (dt ** 2)
body.position += 0.5 * body.velocity * dt body.position += body.velocity * dt
print(body) print(body)
# print(body.name, body.position) # print(body.name, body.position)
...@@ -97,9 +98,28 @@ class Simulator: ...@@ -97,9 +98,28 @@ class Simulator:
self.evolve(dt) self.evolve(dt)
if __name__ == '__main__': def solar_system():
t = 60 * 60 * 24 * 100 t = SECONDS_PER_WEEK
_sys = System([
Sun(), # 太阳
Mercury(), # 水星
Venus(), # 金星
Earth(), # 地球
Mars(), # 火星
Jupiter(), # 木星
Saturn(), # 土星
Uranus(), # 天王星
Neptune(), # 海王星
Pluto() # 冥王星(从太阳系的行星中排除)
])
_sim = Simulator(_sys, t)
_sim.run(t * 100)
if __name__ == '__main__':
# t = 60 * 60 * 24 * 10
solar_system()
# _sys = System([Sun(init_position=[0, 0, 149597870.700]), Earth()]) # _sys = System([Sun(init_position=[0, 0, 149597870.700]), Earth()])
# _sys = System([Sun(init_position=[849597870.700, 0, 0], init_velocity=[0, 9.79, 0]), # _sys = System([Sun(init_position=[849597870.700, 0, 0], init_velocity=[0, 9.79, 0]),
...@@ -112,7 +132,7 @@ if __name__ == '__main__': ...@@ -112,7 +132,7 @@ if __name__ == '__main__':
# Sun(name="sun3", init_position=[0, -849597870.700, 0], init_velocity=[18.0, 0, 0]), # Sun(name="sun3", init_position=[0, -849597870.700, 0], init_velocity=[18.0, 0, 0]),
# Earth(init_position=[0, -349597870.700, 0], init_velocity=[15.50, 0, 0])]) # Earth(init_position=[0, -349597870.700, 0], init_velocity=[15.50, 0, 0])])
_sys = System([Sun(), Earth(), Jupiter()]) # _sys = System([Sun(), Earth(), Jupiter()])
#
_sim = Simulator(_sys, t) # _sim = Simulator(_sys, t)
_sim.run(t * 100) # _sim.run(t * 100)
...@@ -7,6 +7,7 @@ ...@@ -7,6 +7,7 @@
# python_version :3.8 # python_version :3.8
# ============================================================================== # ==============================================================================
import numpy as np import numpy as np
from common.consts import AU, G
class System(object): class System(object):
...@@ -36,14 +37,14 @@ class System(object): ...@@ -36,14 +37,14 @@ class System(object):
ke = 0.5 * body.mass * (v[0] ** 2 + v[1] ** 2) # ke = 0.5 * body.mass * (v[0]**2 v[1]**2) ke = 0.5 * body.mass * (v[0] ** 2 + v[1] ** 2) # ke = 0.5 * body.mass * (v[0]**2 v[1]**2)
return ke return ke
def kinetic_energy(self, vx, vy, vz, mass): # def kinetic_energy(self, vx, vy, vz, mass):
""" # """
计算动能 # 计算动能
""" # """
v = body.velocity # velocity # v = body.velocity # velocity
return 0.5 * mass * (vx ** 2 + vy ** 2 + vz ** 2) # return 0.5 * mass * (vx ** 2 + vy ** 2 + vz ** 2)
def potential_energy(self, x, y, z, mass, G, M): def potential_energy(self, x, y, z, mass,M):
""" """
计算势能 计算势能
""" """
...@@ -139,8 +140,7 @@ class System(object): ...@@ -139,8 +140,7 @@ class System(object):
# body1.position = body1.position + (body1.momentum / body1.mass) * dt # body1.position = body1.position + (body1.momentum / body1.mass) * dt
# body1.update_source_data() # body1.update_source_data()
def acceleration(self, body1, body2):
def acceleration(self, body1, body2, G=6.67e-11):
""" """
计算两个天体之间的加速度 计算两个天体之间的加速度
:param body1: 天体1 :param body1: 天体1
...@@ -155,7 +155,7 @@ class System(object): ...@@ -155,7 +155,7 @@ class System(object):
distance = math.sqrt(dx ** 2 + dy ** 2 + dz ** 2) distance = math.sqrt(dx ** 2 + dy ** 2 + dz ** 2)
mag = G * body2.mass / (distance ** 2) mag = G * body2.mass / (distance ** 2)
acc = [dx / distance * mag, dy / distance * mag, dz / distance * mag] acc = [dx / distance * mag, dy / distance * mag, dz / distance * mag]
return np.array(acc)/1000/1000/1000 return np.array(acc) / 1000 / 1000 / 1000
def calc_acceleration(self): def calc_acceleration(self):
for body1 in self.bodies: for body1 in self.bodies:
...@@ -170,7 +170,7 @@ class System(object): ...@@ -170,7 +170,7 @@ class System(object):
# m/s² = kg * m / m**3 # m/s² = kg * m / m**3
# km/s² = kg * m / m**3 / 1e9 # km/s² = kg * m / m**3 / 1e9
# acceleration = G * body2.mass * dx / (d ** 3) # acceleration = G * body2.mass * dx / (d ** 3)
acceleration += (6.67e-11 * body2.mass * r / np.linalg.norm(r) ** 3) / 1e9 acceleration += (G * body2.mass * r / np.linalg.norm(r) ** 3) / 1e9
body1.acceleration = acceleration body1.acceleration = acceleration
......
Markdown is supported
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册