提交 cf5948b7 编写于 作者: M march3

三体运行模拟器

上级
from bodies.body import Body
from bodies.earth import Earth
from bodies.sun import Sun
\ No newline at end of file
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
import json
import numpy as np
import math
"""
天文单位
"""
AU: float = 149597870.700
class Body:
def __init__(self, name, mass, init_position, init_velocity,
density=5e3, color=(125 / 255, 125 / 255, 125 / 255),
texture=None, size_scale=1.0, distance_scale=1.0):
"""
天体类
:param name: 天体名称
:param mass: 天体质量 (kg)
:param init_position: 初始位置 (km)
:param init_velocity: 初始速度 (km/s)
:param density: 平均密度 (kg/m³)
:param color: 天体颜色(纹理图片优先)
:param texture: 纹理图片
:param size_scale: 尺寸缩放
:param distance_scale: 距离缩放
"""
self.name = name
self.__mass = mass
self.init_position = np.array(init_position, dtype='float32')
self.init_velocity = np.array(init_velocity, dtype='float32')
self.position = self.init_position
self.velocity = self.init_velocity
self.__density = density
self.color = color
self.texture = texture
self.size_scale = size_scale
self.distance_scale = distance_scale
# 初始化后,加速度为0,只有多个天体的引力才会影响到加速度
self.acceleration = np.array([0, 0, 0], dtype='float32')
self.__his_pos = []
self.__his_vel = []
self.__his_acc = []
self.__his_reserved_num = 10
def append_history(self, his_list, data):
"""
:param his_list:
:param data:
:return:
"""
# 如果历史记录为0 或者 新增数据和最后的历史数据不相同,则添加
if len(his_list) == 0 or \
np.sum(data == his_list[-1]) < len(data):
his_list.append(data.copy())
def record_history(self):
"""
记录历史
:return:
"""
# 如果历史记录数超过了保留数量,则截断,只保留 __his_reserved_num 数量的历史
if len(self.__his_pos) > self.__his_reserved_num:
self.__his_pos = self.__his_pos[len(self.__his_pos)-self.__his_reserved_num:]
self.__his_vel = self.__his_vel[len(self.__his_vel)-self.__his_reserved_num:]
self.__his_acc = self.__his_acc[len(self.__his_acc)-self.__his_reserved_num:]
# 追加历史记录(位置、速度、加速度)
self.append_history(self.__his_pos, self.position)
self.append_history(self.__his_vel, self.velocity)
self.append_history(self.__his_acc, self.acceleration)
print(self.name, "his pos->", self.__his_pos)
def his_position(self):
"""
历史位置
:return:
"""
return self.__his_pos
def his_velocity(self):
"""
历史瞬时速度
:return:
"""
return self.__his_vel
def his_acceleration(self):
"""
历史瞬时加速度
:return:
"""
return self.__his_acc
@property
def mass(self):
"""
天体质量 (kg)
:return:
"""
return self.__mass
@property
def density(self):
"""
平均密度 (kg/m³)
:return:
"""
return self.__density
@property
def volume(self):
"""
天体的体积
"""
# v = m/ρ
# 体积(m³) = 质量(kg) / 密度(kg/m³)
v = self.mass / self.density
return v
@property
def raduis(self):
"""
天体的半径
:return:
"""
# V = ⁴⁄₃πr³ -> r = pow((3V)/(4π),1/3)
return pow(3 * self.volume / (4 * math.pi), 1 / 3)
@property
def diameter(self):
"""
天体的直径
:return:
"""
return self.raduis * 2
def __repr__(self):
return '<%s> m=%.3e(kg), r=%.3e(km), p=[%.3e,%.3e,%.3e](km), v=%s(km/s)' % \
(self.name, self.mass, self.raduis,
self.position[0], self.position[1], self.position[2], self.velocity)
def position_au(self):
pos = self.position
pos_au = pos / AU
return pos_au
def change_velocity(self, dv):
self.velocity += dv
def move(self, dt):
self.position += self.velocity * dt
def reset(self):
self.position = self.init_position
self.velocity = self.init_velocity
def kinetic_energy(self):
"""
计算动能(千焦耳)
表示动能,单位为焦耳j,m为质量,单位为千克,v为速度,单位为米/秒。
ek=(1/2).m.v^2
m(kg) v(m/s) -> j
m(kg) v(km/s) -> kj
"""
v = self.velocity
return 0.5 * self.mass * (v[0] ** 2 + v[1] ** 2 + v[2] ** 2)
@staticmethod
def build_bodies_from_json(json_file):
bodies = []
with open(json_file, "r") as read_content:
json_data = json.load(read_content)
for body_data in json_data["bodies"]:
# print(body_data)
body = Body(**body_data)
bodies.append(body)
# print(body.position_au())
return bodies
if __name__ == '__main__':
# build_bodies_from_json('../data/sun.json')
bodies = Body.build_bodies_from_json('../data/sun_earth.json')
# 太阳半径 / 地球半径
print(bodies[0].raduis / bodies[1].raduis)
for body in bodies:
print(body.kinetic_energy())
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
from bodies.body import Body
class Earth(Body):
"""
地球
------------------------
远日点距离: 152097701 km
近日点距离: 147098074 km
 逃逸速度: 11.186 km/s
 公转速度: 29.79 km/s
 天体质量: 5.97237✕10²⁴ kg
 平均密度: 5507.85 kg/m³
"""
def __init__(self, name="earth", mass=5.97237e24,
init_position=[149597870.700, 0, 0],
init_velocity=[29.79, 0, 0],
texture="", size_scale=1.0, distance_scale=1.0):
params = {
"name": name,
"mass": mass,
"init_position": init_position,
"init_velocity": init_velocity,
"density": 5507.85,
"color": [
125,
125,
125
],
"texture": texture,
"size_scale": size_scale,
"distance_scale": distance_scale
}
super().__init__(**params)
if __name__ == '__main__':
earth = Earth()
print(earth)
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
from bodies.body import Body
class Sun(Body):
def __init__(self, name="sun", mass=1.9891e30, init_position=[0, 0, 0], init_velocity=[0, 0, 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.408e3,
"color": [
125,
125,
125
],
"texture": texture,
"size_scale": size_scale,
"distance_scale": distance_scale
}
super().__init__(**params)
if __name__ == '__main__':
print(Sun())
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
import numpy as np
import math
class MechanicalCalculator:
"""
力学计算器:
牛顿引力常数就有万有引力常数:G=6.67x10^-11 (N·m^2 /kg^2)
计算万有引力时会用到
万有引力公式:F=G*m1*m2/r*
"""
@staticmethod
def getAcc(pos, mass, G, softening):
"""
Calculate the acceleration on each particle due to Newton's Law
pos is an N x 3 matrix of positions
mass is an N x 1 vector of masses
G is Newton's Gravitational constant
softening is the softening length
a is N x 3 matrix of accelerations
"""
# positions r = [x,y,z] for all particles
x = pos[:, 0:1]
y = pos[:, 1:2]
z = pos[:, 2:3]
# matrix that stores all pairwise particle separations: r_j - r_i
dx = x.T - x
dy = y.T - y
dz = z.T - z
# matrix that stores 1/r^3 for all particle pairwise particle separations
inv_r3 = (dx ** 2 + dy ** 2 + dz ** 2 + softening ** 2)
inv_r3[inv_r3 > 0] = inv_r3[inv_r3 > 0] ** (-1.5)
ax = G * (dx * inv_r3) @ mass
ay = G * (dy * inv_r3) @ mass
az = G * (dz * inv_r3) @ mass
# pack together the acceleration components
a = np.hstack((ax, ay, az))
return a
@staticmethod
def getEnergy(pos, vel, mass, G):
"""
Get kinetic energy (KE) and potential energy (PE) of simulation
pos is N x 3 matrix of positions
vel is N x 3 matrix of velocities
mass is an N x 1 vector of masses
G is Newton's Gravitational constant
KE is the kinetic energy of the system
PE is the potential energy of the system
"""
# Kinetic Energy:
KE = 0.5 * np.sum(np.sum(mass * vel ** 2))
# Potential Energy:
# positions r = [x,y,z] for all particles
x = pos[:, 0:1]
y = pos[:, 1:2]
z = pos[:, 2:3]
# matrix that stores all pairwise particle separations: r_j - r_i
dx = x.T - x
dy = y.T - y
dz = z.T - z
# matrix that stores 1/r for all particle pairwise particle separations
inv_r = np.sqrt(dx ** 2 + dy ** 2 + dz ** 2)
inv_r[inv_r > 0] = 1.0 / inv_r[inv_r > 0]
# sum over upper triangle, to count each interaction only once
PE = G * np.sum(np.sum(np.triu(-(mass * mass.T) * inv_r, 1)))
return KE, PE;
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
from common.system import System
from bodies import Body, Sun, Earth
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
X_MIN, X_MAX = -2e+14, 2e+14 # the x range of the bounding box (m)
Y_MIN, Y_MAX = -2e+14, 2e+14 # the y range of the bounding box (m)
Z_MIN, Z_MAX = -2e+14, 2e+14 # the z range of the bounding box (m)
def show(bodies, idx=0):
# Creating figures for the plot
fig = plt.figure(figsize=(16, 12))
ax = plt.axes(projection="3d")
ax.set_xlim(X_MIN, X_MAX)
ax.set_ylim(Y_MIN, Y_MAX)
ax.set_zlim(Z_MIN, Z_MAX)
# Creating a plot using the random datasets
colors = ['red', 'blue']
sizes = [800, 500]
for idx, body in enumerate(bodies):
pos = body.position
ax.scatter(pos[0], pos[1], pos[2], color=colors[idx], s=sizes[idx])
for _his_pos in body.his_position():
ax.scatter3D(_his_pos[0], _his_pos[1], _his_pos[2], color=colors[idx], alpha=0.5)
# ax.scatter(his_pos[0], his_pos[1], his_pos[2], color=colors[idx], s=10)
his_pos = body.his_position()
if len(his_pos) > 2:
_his_pos = list(zip(*his_pos))
ax.plot3D(_his_pos[0], _his_pos[1], _his_pos[2], color=colors[idx], alpha=0.5)
plt.title("3D scatter plot %s" % idx)
# display the plot
plt.show()
import time
class Simulator:
def __init__(self, system, dt):
self.system = system
self.dt = dt
def evolve(self, dt):
self.system.update_acceleration()
# next_p = p + dt * v + 0.5 * a * (dt ** 2)
# next_v = v + dt * a
# for body1 in bodies:
# for body2 in bodies:
# if body1 == body2: # 相等说明是同一个天体,跳过计算
# continue
# r = body2.position - body1.position
# # F = G x (m1 x m2) / r² 万有引力定律
# F = Configs.G * body1.mass * body2.mass * r / np.linalg.norm(r) / r.dot(r)
# body1.momentum = body1.momentum + F * dt # 动量定理
# body1.position = body1.position + (body1.momentum / body1.mass) * dt
# body1.update_source_data()
for body in self.system.bodies:
body.position += body.velocity * dt + 0.5 * body.acceleration * (dt ** 2)
# acceleration 加速度
body.velocity += body.acceleration * dt
print(body.name, body.position)
body.record_history()
show(self.system.bodies)
def run(self, t):
n = int(t / self.dt)
for i in range(n):
time.sleep(1)
self.evolve(self.dt)
def run_always(self):
dt = 1 # 时间变化1秒
while True:
time.sleep(1)
self.evolve(dt)
if __name__ == '__main__':
t = 60 * 60 * 24
_sys = System([Sun(), Earth()])
_sim = Simulator(_sys, t)
_sim.run(t * 100)
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
import numpy as np
class System(object):
def __init__(self, bodies):
self.bodies = bodies
def add(self, body):
self.bodies.append(body)
def total_mass(self):
total_mass = 0
for body in self.bodies:
total_mass = body.mass
return total_mass
def __repr__(self):
return 'System({})'.format(self.bodies)
def kinetic_energy(self):
"""
KE是系统的动能
:return:
"""
ke = 0
for body in self.bodies:
v = body.velocity # velocity
ke = 0.5 * body.mass * (v[0] ** 2 + v[1] ** 2) # ke = 0.5 * body.mass * (v[0]**2 v[1]**2)
return ke
def kinetic_energy(self, vx, vy, vz, mass):
"""
计算动能
"""
v = body.velocity # velocity
return 0.5 * mass * (vx ** 2 + vy ** 2 + vz ** 2)
def potential_energy(self, x, y, z, mass, G, M):
"""
计算势能
"""
return -G * M * mass / np.sqrt(x ** 2 + y ** 2 + z ** 2)
def potential_energy(self):
"""
PE是系统的势能
:return:
"""
pe = 0
for body1 in self.bodies:
for body2 in self.bodies:
if body1 is body2:
continue
r = body1.position - body2.position
# pe -= body1.mass * body2.mass / np.sqrt(r[0]**2 r[1]**2)
pe -= body1.mass * body2.mass / np.sqrt(r[0] ** 2 + r[1] ** 2)
return pe
def momentum(self):
"""
动量
:return:
"""
p = np.zeros(2)
for body in self.bodies:
p = body.mass * body.velocity
return p
def momentum(self, vx, vy, vz, mass):
"""
计算动量
"""
return mass * np.sqrt(vx ** 2 + vy ** 2 + vz ** 2)
def angular_momentum(self):
"""
角动量
:return:
"""
L = 0
for body in self.bodies:
r = body.position
v = body.velocity
L = body.mass * (r[0] * v[1] - r[1] * v[0])
return L
def angular_momentum(self, x, y, z, vx, vy, vz):
"""
计算角动量
"""
return self.mass * (x * vy - y * vx)
def center_of_mass(self):
r = np.zeros(2)
for body in self.bodies:
r = body.mass * body.position
return r / self.total_mass()
# def step(m: np.ndarray, v: np.ndarray, p: np.ndarray):
# ''' Calculate the next state of the N-star system.
# args:
# m: np.ndarray[(N,), np.float64], mass.
# v: np.ndarray[(N,3), np.float64], velocity.
# p: np.ndarray[(N,3), np.float64], position.
# returns:
# (next_v, next_p), the next state.
# '''
# N = m.shape[0]
# a = np.zeros_like(p)
# for i in range(N):
# for j in range(N):
# if i == j: continue
# r = np.sqrt(np.sum(np.square(p[j, :] - p[i, :])))
# aij = G * m[j] / (r * r)
# dir = (p[j, :] - p[i, :]) / r
# a[i, :] += aij * dir
# next_p = p + (1 / FS) * v + 0.5 * a * ((1 / FS) ** 2)
# next_v = v + (1 / FS) * a
# return (next_v, next_p)
# next_p = p + dt * v + 0.5 * a * (dt ** 2)
# next_v = v + dt * a
# for body1 in bodies:
# for body2 in bodies:
# if body1 == body2: # 相等说明是同一个天体,跳过计算
# continue
# r = body2.position - body1.position
# # F = G x (m1 x m2) / r² 万有引力定律
# F = Configs.G * body1.mass * body2.mass * r / np.linalg.norm(r) / r.dot(r)
# body1.momentum = body1.momentum + F * dt # 动量定理
# body1.position = body1.position + (body1.momentum / body1.mass) * dt
# body1.update_source_data()
def update_acceleration(self):
for body in self.bodies:
body.acceleration = np.zeros(3)
for body1 in self.bodies:
body1.record_history()
for body2 in self.bodies:
if body1 is body2:
continue
# r = np.sqrt(np.sum(np.square(body2.position - body1.position)))
# aij = 6.67e-11 * body2.mass / (r * r)
# dir = (body2.position - body1.position) / r
# body.acceleration1 += aij * dir
r = body2.position - body1.position
body.acceleration += -6.67e-11 * body2.mass * r / np.linalg.norm(r) ** 3
# body1.acceleration = body.acceleration2
# r = body2.position - body1.position
# # body1.acceleration = -body2.mass * r / np.linalg.norm(r) ** 3
# body1.acceleration += 6.67e-11 * body2.mass * r / np.linalg.norm(r) ** 3
# body1.acceleration = (G * body2.mass * (x / pow(x ** 2 + y ** 2 + z ** 2, 3 / 2)))
# pow(x ** 2 + y ** 2 + z ** 2, 3 / 2)
# G = 6.67e-11 # 万有引力常数
# m1 = 1 # 第一个天体的质量
# m2 = 1 # 第二个天体的质量
# r = 1 # 两个天体之间的距离
#
# a = G * m2 / math.pow(r, 2)
"""
# Python 3
def get_acceleration(mass1,mass2,position1,position2):
G = 6.67*10**-11
r = np.linalg.norm(position2-position1)
return G*mass2*(position2-position1)/r**3
# Python 2
# def get_acceleration(mass1,mass2,position1,position2):
# G = 6.67*10**-11
# r = np.linalg.norm(position2-position1)
# return G*mass2*(position2-position1)/r**3
根据两个天体body1 和 body2 的质量,距离 position[0],position[1],position[2],用python计算两个天体的加速度
# 假设两个天体的质量分别为m1和m2
m1 = 5.974 * 10 ** 24
m2 = 7.348 * 10 ** 22
# 计算两个天体的位置,假设分别为position1 和 position2
position1 = [2.21, 3.45, 4.67]
position2 = [1.78, 6.34, -1.23]
# 计算两个天体的距离
dist = ( (position1[0]-position2[0])**2 + (position1[1]-position2[1])**2 + (position1[2]-position2[2])**2 ) ** (1/2)
# 计算两个天体的加速度
G = 6.674 * 10 ** -11
acc_1 = G * m2 / (dist**2)
acc_2 = G * m1 / (dist**2)
# 输出结果
print("body1的加速度为:", acc_1)
print("body2的加速度为:", acc_2)
ax1 = (G * m2 * (x/pow(x**2 + y**2 + z**2,3/2)))
ay1 = (G * m2 * (y/pow(x**2 + y**2 + z**2,3/2)))
az1 = (G * m2 * (z/pow(x**2 + y**2 + z**2,3/2)))
ax2 = (G * m1 * (-x/pow(x**2 + y**2 + z**2,3/2)))
ay2 = (G * m1 * (-y/pow(x**2 + y**2 + z**2,3/2)))
az2 = (G * m1 * (-z/pow(x**2 + y**2 + z**2,3/2)))
加速度可以用万有引力定律计算:
a = G * m2 / r2
其中G为万有引力常数,m2和r2分别为两个天体的质量和距离平方。
因此,可以使用以下python代码来计算两个天体的加速度:
import math
G = 6.67e-11 # 万有引力常数
m1 = 1 # 第一个天体的质量
m2 = 1 # 第二个天体的质量
r = 1 # 两个天体之间的距离
a = G * m2 / math.pow(r, 2)
print(a)
"""
{
"bodies": [
{
"name": "sun",
"mass": 1.9891e30,
"init_position": [
0,
0,
0
],
"init_velocity": [
0,
0,
0
],
"density": 1.408e3,
"color": [
125,
125,
125
],
"texture": "",
"size_scale": 1.0,
"distance_scale": 1.0
}
]
}
\ No newline at end of file
{
"bodies": [
{
"name": "sun",
"mass": 1.9891e30,
"init_position": [
0,
0,
0
],
"init_velocity": [
0,
0,
0
],
"density": 1.408e3,
"color": [
125,
125,
125
],
"texture": "",
"size_scale": 1.0,
"distance_scale": 1.0
},
{
"name": "earth",
"mass": 5.97237e24,
"init_position": [
149597870.700,
0,
0
],
"init_velocity": [
29.79,
0,
0
],
"density": 5507.85,
"color": [
125,
125,
125
],
"texture": "",
"size_scale": 1.0,
"distance_scale": 1.0
}
]
}
\ No newline at end of file
# -*- coding:utf-8 -*-
# title :
# description :
# author :Python超人
# date :2023-01-22
# notes :
# python_version :3.8
# ==============================================================================
# from mayavi import mlab
from bodies import Sun, Earth
bodies = [
Sun(size_scale=1.5e-2),
Earth(size_scale=9e3) # 地球
]
print(bodies)
# 运行天体的动画
# run_anim(bodies)
# azimuth:
# 观测方位角,可选,float类型(以度为单位,0-360),用x轴投影到x-y平面上的球体上的位置矢量所对的角度。
# elevation:
# 观测天顶角,可选,float类型(以度为单位,0-180), 位置向量和z轴所对的角度。
# distance:
# 观测距离,可选,float类型 or 'auto',一个正浮点数,表示距放置相机的焦点的距离。
# Mayavi 3.4.0中的新功能:'auto' 使得距离为观察所有对象的最佳位置。
# focalpoint:
# 观测焦点,可选,类型为一个由3个浮点数组成的数组 or 'auto',,代表观测相机的焦点
# Mayavi 3.4.0中的新功能:'auto',则焦点位于场景中所有对象的中心。
# roll:
# 控制滚动,可选,float类型,即摄影机围绕其轴的旋转
# reset_roll:
# 布尔值,可选。如果为True,且未指定“滚动”,则重置相机的滚动方向。
# figure:
# 要操作的Mayavi图形。如果为 None,则使用当前图形。
# mlab.view(azimuth=-45, elevation=45, distance=50e8 * 2 * 2 * 4 * 4, focalpoint=[5e10, 0, 0])
# mlab.show()
Markdown is supported
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册