From 4bafb14aab9414920f3ace04cc5c32ed8c73f989 Mon Sep 17 00:00:00 2001 From: march3 Date: Mon, 24 Apr 2023 15:44:38 +0800 Subject: [PATCH] =?UTF-8?q?Python=E8=B6=85=E4=BA=BA-=E5=AE=87=E5=AE=99?= =?UTF-8?q?=E6=A8=A1=E6=8B=9F=E5=99=A8?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- common/func.py | 300 ++++++++++++++++++--------------- common/system.py | 2 +- sim_lab/lagrangian_points.py | 6 +- sim_lab/lagrangian_points_3.py | 254 ++++++++++++++++++++++++++++ sim_scenes/func.py | 21 ++- simulators/calc_simulator.py | 139 +++++++++++++++ 6 files changed, 577 insertions(+), 145 deletions(-) create mode 100644 sim_lab/lagrangian_points_3.py create mode 100644 simulators/calc_simulator.py diff --git a/common/func.py b/common/func.py index 9800bde..f34c543 100644 --- a/common/func.py +++ b/common/func.py @@ -108,142 +108,164 @@ def calculate_distance(pos1, pos2=[0, 0, 0]): pow(np.array(pos1[2]) - np.array(pos2[2]), 2), 1 / 2) return d - -def calculate_velocity(mass, semimajor_axis, eccentricity): - """ - 计算天体在椭圆轨道上的速度。 - - 参数： - - mass: 天体质量，单位 kg - - semimajor_axis: 轨道半长轴，单位 m - - eccentricity: 轨道离心率 - - 返回值： - 天体在轨道上的速度，单位 m/s。 - """ - # 计算轨道的半短轴和半焦距 - semiminor_axis = semimajor_axis * math.sqrt(1 - eccentricity ** 2) - focus_distance = semimajor_axis * eccentricity - - # 计算轨道的第一和第二离心率角 - theta = math.atan2(focus_distance, semiminor_axis) - psi = math.atan2(math.sqrt(1 - eccentricity ** 2) * math.sin(theta), eccentricity + math.cos(theta)) - - # 计算轨道的速率 - v = math.sqrt(G * mass * (2 / semimajor_axis - 1 / semiminor_axis)) - - # 计算天体在轨道上的速度，注意要考虑太阳的质量 - return v * math.sqrt(1 + (2 * mass) / (v ** 2 * AU)) * math.sin(psi) - - -def get_lagrange_points(m1, m2, r, R, omega): - """ - 函数需要5个输入参数： - - m1: 大质点的质量（单位：kg） - m2: 小质点的质量（单位：kg） - r: 小质点到拉格朗日点的距离（单位：km） - R: 大质点到拉格朗日点的距离（单位：km） - omega: 两个星体之间的夹角（单位：rad） - 函数返回一个元组，其中包括五个元素（Tuple），每个元素又是由七个数据（Tuple）组成： - - @param m1: - @param m2: - @param r: - @param R: - @param omega: - @return: - x: 拉格朗日点的x坐标（单位：km） - y: 拉格朗日点的y坐标（单位：km） - z: 拉格朗日点的z坐标（单位：km） - vx: 拉格朗日点物体的x方向速度（单位：km/s） - vy: 拉格朗日点物体的y方向速度（单位：km/s） - vz: 拉格朗日点物体的z方向速度（单位：km/s） - """ - G = 6.673e-20 # gravitational constant in km^3/kg*s^2 - M = m1 + m2 - - # Calculate mass ratio - mu = m2 / M - - # Calculate distance between primary and secondary bodies - d = np.sqrt((R * mu) ** 2 + r ** 2 - 2 * R * r * mu * np.cos(omega)) - - # Calculate L1 point - a = (d - R) / (2 - mu) - x1 = R - a - y1 = 0 - z1 = 0 - v1 = np.sqrt(G * m2 * a / d) * (1 - mu) - vx1 = 0 - vy1 = v1 * (R - x1) / d - vz1 = v1 * y1 / d - - # Calculate L2 point - a = (d - R) / (2 - mu) - x2 = R + a - y2 = 0 - z2 = 0 - v2 = np.sqrt(G * m2 * a / d) * (1 - mu) - vx2 = 0 - vy2 = v1 * (R - x2) / d - vz2 = v1 * y2 / d - - # Calculate L3 point - a = (d + R) / (2 + mu) - x3 = -a - y3 = 0 - z3 = 0 - v3 = np.sqrt(G * m1 * a / d) * (1 + mu) - vx3 = 0 - vy3 = v3 * (R - x3) / d - vz3 = v3 * -y3 / d - - # Calculate L4 and L5 points - x4, y4, z4, v4, vx4, vy4, vz4 = get_l4_l5_points(m1, m2, R, omega, 1) - x5, y5, z5, v5, vx5, vy5, vz5 = get_l4_l5_points(m1, m2, R, omega, -1) - - return ((x1, y1, z1, vx1, vy1, vz1), - (x2, y2, z2, vx2, vy2, vz2), - (x3, y3, z3, vx3, vy3, vz3), - (x4, y4, z4, vx4, vy4, vz4), - (x5, y5, z5, vx5, vy5, vz5)) - - -def get_l4_l5_points(m1, m2, R, omega, sign): - # G = 6.673e-20 # gravitational constant in km^3/kg*s^2 - M = m1 + m2 - - # Calculate mass ratio - mu = m2 / M - - # Calculate position of L4 or L5 point - x = R * np.cos(omega + sign * 60 * np.pi / 180) - y = R * np.sin(omega + sign * 60 * np.pi / 180) - z = 0 - - # Calculate velocity of L4 or L5 point - v = np.sqrt(G * M / (3 * R)) - vx = -v * y / R - vy = v * x / R - vz = 0 - - return x, y, z, v, vx, vy, vz - - -if __name__ == '__main__': - # print(calculate_distance([6, 8, 0], [3, 4, 0])) - # print(find_file("common/func.py")) - - # 使用地球数据测试 - mass_earth = 5.972e24 - semimajor_axis_earth = AU - eccentricity_earth = 0.0167 - - velocity_earth = calculate_velocity(mass_earth, semimajor_axis_earth, eccentricity_earth) - - print("地球在轨道上的速度是：{:.2f} km/s".format(velocity_earth / 1000)) - -""" -你现在是一位资深的天体学家、请使用python完成一个获取拉格朗日点的函数，获取拉格朗日点的坐标(px,py 单位：km)同时，也要返回拉格朗日点物体的速度（vx,vy 单位：km/s）,需要返回L1-L5所有的点和速度。返回的L1、L2、L3、L4、L5 格式为：(x1, y1, vx1, vy1),(x2, y2, vx2, vy2), (x3, y3, vx3, vy3), (x4, y4, vx4, vy4) (x5, y5, vx5, vy5)，并针对太阳、地球拉格朗日点以及 地球、月球拉格朗日点举例。并使用matlibplot画图，并写上详细的注释。 -""" +# +# def calculate_velocity(mass, semimajor_axis, eccentricity): +# """ +# 计算天体在椭圆轨道上的速度。 +# +# 参数： +# - mass: 天体质量，单位 kg +# - semimajor_axis: 轨道半长轴，单位 m +# - eccentricity: 轨道离心率 +# +# 返回值： +# 天体在轨道上的速度，单位 m/s。 +# """ +# # 计算轨道的半短轴和半焦距 +# semiminor_axis = semimajor_axis * math.sqrt(1 - eccentricity ** 2) +# focus_distance = semimajor_axis * eccentricity +# +# # 计算轨道的第一和第二离心率角 +# theta = math.atan2(focus_distance, semiminor_axis) +# psi = math.atan2(math.sqrt(1 - eccentricity ** 2) * math.sin(theta), eccentricity + math.cos(theta)) +# +# # 计算轨道的速率 +# v = math.sqrt(G * mass * (2 / semimajor_axis - 1 / semiminor_axis)) +# +# # 计算天体在轨道上的速度，注意要考虑太阳的质量 +# return v * math.sqrt(1 + (2 * mass) / (v ** 2 * AU)) * math.sin(psi) + +# +# def get_lagrange_points(m1, m2, r, R, omega): +# """ +# 函数需要5个输入参数： +# +# m1: 大质点的质量（单位：kg） +# m2: 小质点的质量（单位：kg） +# r: 小质点到拉格朗日点的距离（单位：km） +# R: 大质点到拉格朗日点的距离（单位：km） +# omega: 两个星体之间的夹角（单位：rad） +# 函数返回一个元组，其中包括五个元素（Tuple），每个元素又是由七个数据（Tuple）组成： +# +# @param m1: +# @param m2: +# @param r: +# @param R: +# @param omega: +# @return: +# x: 拉格朗日点的x坐标（单位：km） +# y: 拉格朗日点的y坐标（单位：km） +# z: 拉格朗日点的z坐标（单位：km） +# vx: 拉格朗日点物体的x方向速度（单位：km/s） +# vy: 拉格朗日点物体的y方向速度（单位：km/s） +# vz: 拉格朗日点物体的z方向速度（单位：km/s） +# """ +# G = 6.673e-20 # gravitational constant in km^3/kg*s^2 +# M = m1 + m2 +# +# # Calculate mass ratio +# mu = m2 / M +# +# # Calculate distance between primary and secondary bodies +# d = np.sqrt((R * mu) ** 2 + r ** 2 - 2 * R * r * mu * np.cos(omega)) +# +# # Calculate L1 point +# a = (d - R) / (2 - mu) +# x1 = R - a +# y1 = 0 +# z1 = 0 +# v1 = np.sqrt(G * m2 * a / d) * (1 - mu) +# vx1 = 0 +# vy1 = v1 * (R - x1) / d +# vz1 = v1 * y1 / d +# +# # Calculate L2 point +# a = (d - R) / (2 - mu) +# x2 = R + a +# y2 = 0 +# z2 = 0 +# v2 = np.sqrt(G * m2 * a / d) * (1 - mu) +# vx2 = 0 +# vy2 = v1 * (R - x2) / d +# vz2 = v1 * y2 / d +# +# # Calculate L3 point +# a = (d + R) / (2 + mu) +# x3 = -a +# y3 = 0 +# z3 = 0 +# v3 = np.sqrt(G * m1 * a / d) * (1 + mu) +# vx3 = 0 +# vy3 = v3 * (R - x3) / d +# vz3 = v3 * -y3 / d +# +# # Calculate L4 and L5 points +# x4, y4, z4, v4, vx4, vy4, vz4 = get_l4_l5_points(m1, m2, R, omega, 1) +# x5, y5, z5, v5, vx5, vy5, vz5 = get_l4_l5_points(m1, m2, R, omega, -1) +# +# return ((x1, y1, z1, vx1, vy1, vz1), +# (x2, y2, z2, vx2, vy2, vz2), +# (x3, y3, z3, vx3, vy3, vz3), +# (x4, y4, z4, vx4, vy4, vz4), +# (x5, y5, z5, vx5, vy5, vz5)) +# +# +# def get_l4_l5_points(m1, m2, R, omega, sign): +# # G = 6.673e-20 # gravitational constant in km^3/kg*s^2 +# M = m1 + m2 +# +# # Calculate mass ratio +# mu = m2 / M +# +# # Calculate position of L4 or L5 point +# x = R * np.cos(omega + sign * 60 * np.pi / 180) +# y = R * np.sin(omega + sign * 60 * np.pi / 180) +# z = 0 +# +# # Calculate velocity of L4 or L5 point +# v = np.sqrt(G * M / (3 * R)) +# vx = -v * y / R +# vy = v * x / R +# vz = 0 +# +# return x, y, z, v, vx, vy, vz +# +# +# def get_v(M, m, r): +# import math +# +# G = 6.67e-11 # 引力常数，单位为N·m²/kg² +# v = math.sqrt(G * M / r) # 计算地球的线速度，单位为米/秒 +# return v +# +# +# if __name__ == '__main__': +# M = 5.97237e24 +# r = 363104*1000 +# m = 5.97e24 +# print(get_v(M, m, r)) + # import math + # + # G = 6.67e-11 # 引力常数，单位为N·m²/kg² + # M = 1.99e30 # 太阳质量，单位为kg + # # m = 5.97e24 # 地球质量，单位为kg + # r = 1.5e11 # 地球到太阳的距离，单位为米 + # + # v = math.sqrt(G * M / r) # 计算地球的线速度，单位为米/秒 + # + # print("天体的线速度为：%.2f 公里/秒" % (v / 1000)) +# # print(calculate_distance([6, 8, 0], [3, 4, 0])) +# # print(find_file("common/func.py")) +# +# # 使用地球数据测试 +# mass_earth = 5.972e24 +# semimajor_axis_earth = AU +# eccentricity_earth = 0.0167 +# +# velocity_earth = calculate_velocity(mass_earth, semimajor_axis_earth, eccentricity_earth) +# +# print("地球在轨道上的速度是：{:.2f} km/s".format(velocity_earth / 1000)) +# +# """ +# 你现在是一位资深的天体学家、请使用python完成一个获取拉格朗日点的函数，获取拉格朗日点的坐标(px,py 单位：km)同时，也要返回拉格朗日点物体的速度（vx,vy 单位：km/s）,需要返回L1-L5所有的点和速度。返回的L1、L2、L3、L4、L5 格式为：(x1, y1, vx1, vy1),(x2, y2, vx2, vy2), (x3, y3, vx3, vy3), (x4, y4, vx4, vy4) (x5, y5, vx5, vy5)，并针对太阳、地球拉格朗日点以及 地球、月球拉格朗日点举例。并使用matlibplot画图，并写上详细的注释。 +# """ diff --git a/common/system.py b/common/system.py index 41d5452..315126d 100644 --- a/common/system.py +++ b/common/system.py @@ -226,7 +226,7 @@ class System(object): if body1.mass / body2.mass < 1e10: # 10亿倍的差距 self.fast_calc_list[body1].append(body2) else: - print(f"{body2.name}相对{body1.name}质量太小，加速度可以忽略不计！") + # print(f"{body2.name}相对{body1.name}质量太小，加速度可以忽略不计！") continue r = body2.position - body1.position diff --git a/sim_lab/lagrangian_points.py b/sim_lab/lagrangian_points.py index 2d2259d..deb372c 100644 --- a/sim_lab/lagrangian_points.py +++ b/sim_lab/lagrangian_points.py @@ -159,14 +159,14 @@ if __name__ == '__main__': points = get_lagrangian_points(earth.mass, moon.mass, 363104) offset_points = [ - [0, 0, 21590], # 调整加速度为0 + [0, 0, 0], # [0, 0, 21590], # 调整加速度为0 [0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0], ] velocities = [ - [-0.7, -0.1, 0], # [-0.859, 0, 0], + [-0.872, 0, 0], [-1.265, 0, 0], [1.03, 0, 0], [0, 0, 0], @@ -205,7 +205,7 @@ if __name__ == '__main__': # 使用 ursina 查看的运行效果 # 常用快捷键： P：运行和暂停 O：重新开始 I：显示天体轨迹 # position = 左-右+、上+下-、前+后- - ursina_run(bodies, SECONDS_PER_HOUR, + ursina_run(bodies, SECONDS_PER_HOUR*10, position=(-300000, 1500000, -100), show_timer=True, show_trail=True) diff --git a/sim_lab/lagrangian_points_3.py b/sim_lab/lagrangian_points_3.py new file mode 100644 index 0000000..fad36ef --- /dev/null +++ b/sim_lab/lagrangian_points_3.py @@ -0,0 +1,254 @@ +""" +https://www.163.com/dy/article/G5F1016F053102ZV.html +https://www.sciencedirect.com/topics/physics-and-astronomy/lagrangian-points +以下是太阳和地球的第一、二、三个拉格朗日点的真实坐标和速度数据： +L1点： 坐标： x = 0.010205 AU， y = 0 AU， z = 0 AU 速度： vx = 0 m/s， vy = 246.593 m/s， vz = 0 m/s +L2点： 坐标： x = -0.010205 AU， y = 0 AU， z = 0 AU 速度： vx = 0 m/s， vy = -246.593 m/s， vz = 0 m/s +L3点： 坐标： x = 0.990445 AU， y = 0 AU， z = 0 AU 速度： vx = 0 m/s， vy = 11.168 m/s， vz = 0 m/s +L4点： 坐标： x = 0.500 AU， y = 0.866025 AU， z = 0 AU 速度： vx = -2446.292 m/s， vy = -1412.901 m/s， vz = 0 m/s +L5点： 坐标： x = 0.500 AU， y = -0.866025 AU， z = 0 AU 速度： vx = -2446.292 m/s， vy = 1412.901 m/s， vz = 0 m/s +https://baike.baidu.com/pic/%E6%8B%89%E6%A0%BC%E6%9C%97%E6%97%A5%E7%82%B9/731078/0/dbb44aed2e738bd4510fa07aa98b87d6277ff94b?fr=lemma&fromModule=lemma_content-image&ct=single#aid=0&pic=dbb44aed2e738bd4510fa07aa98b87d6277ff94b +""" +# +# AU = 1.496e8 +# +# +# def compute_barycenter(masses, positions): +# """ +# Compute the barycenter position of celestial objects in 3D space +# masses: a list of masses of celestial objects +# positions: a list of positions of celestial objects, each position is a tuple (x, y, z) +# """ +# m_sum = sum(masses) +# x_sum = 0 +# y_sum = 0 +# z_sum = 0 +# for i in range(len(masses)): +# x_sum += masses[i] * positions[i][0] / m_sum +# y_sum += masses[i] * positions[i][1] / m_sum +# z_sum += masses[i] * positions[i][2] / m_sum +# return (x_sum, y_sum, z_sum) +# +# +# def get_lagrangian_points(m1, m2, r): +# """ +# https://baike.baidu.com/item/%E6%8B%89%E6%A0%BC%E6%9C%97%E6%97%A5%E7%82%B9/731078 +# +# @param m1: 大质量 +# @param m2: 小质量 +# @param r: 半径 +# @return: +# """ +# a = m2 / (m1 + m2) +# l1 = (r * (1 - pow(a / 3, 1 / 3)), 0) +# l2 = (r * (1 + pow(a / 3, 1 / 3)), 0) +# l3 = (-r * (1 + (5 * a) / 12), 0) +# l4 = ((r / 2) * ((m1 - m2) / (m1 + m2)), pow(3, 1 / 2) / 2 * r) +# l5 = ((r / 2) * ((m1 - m2) / (m1 + m2)), -pow(3, 1 / 2) / 2 * r) +# +# # print(l1[0]/AU, l2[0]/AU, l3[0]/AU, l4[0]/AU, l5[0]/AU) +# return l1, l2, l3, l4, l5 +# +# +# def show_figure(points, p1_name, p2_name, unit, barycenter=None): +# import matplotlib.pyplot as plt +# +# plt.figure(figsize=(16, 12)) +# plt.plot(0, 0, "ro", markersize=20, label=p1_name) +# plt.text(-unit / 20, -unit / 10, p1_name, fontsize=30, color="r") +# plt.plot(unit, 0, "b.", markersize=4, label=p2_name) +# plt.text(unit - unit / 20, -unit / 20, p2_name, fontsize=20, color="b") +# idx = 1 +# +# for x, y in points: +# plt.plot(x, y, "gx", markersize=3, label=f"L{idx}") +# if idx == 1: +# x_offset = -unit / 22 +# else: +# x_offset = unit / 300 +# plt.text(x + x_offset, y + unit / 300, f"L{idx}", fontsize=18, color="g") +# idx += 1 +# +# if barycenter is not None: +# plt.plot(barycenter[0], barycenter[1], "gx", markersize=10, label=f"L{idx}") +# +# # plt.plot(x, y, "b.") # b：蓝色，.：点 +# # plt.plot(x, y1, "ro") # r：红色，o：圆圈 +# # plt.plot(x, y2, "kx") # k：黑色，x：x字符(小叉) +# plt.show() # 在窗口显示该图片 +# +# +# barycenter = compute_barycenter([5.97237e24, 7.342e22], [[0, 0, 0], [363104, 0, 0]]) +# print(barycenter) +# # show_figure(get_lagrangian_points(1.9891e30, 5.97237e24, AU), "Sun", "Earth", AU) +# show_figure(get_lagrangian_points(5.97237e24, 7.342e22, 363104), "Earth", "Moon", 363104, barycenter) + +# ************************************************************************************************************** +# ************************************************************************************************************** + +# -*- coding:utf-8 -*- +# title :地月场景模拟 +# description :地月场景模拟 +# author :Python超人 +# date :2023-02-11 +# link :https://gitcode.net/pythoncr/ +# python_version :3.8 +# ============================================================================== +from bodies import Sun, Earth, Moon +from objs import Satellite, Satellite2 +from common.consts import SECONDS_PER_HOUR, SECONDS_PER_HALF_DAY, SECONDS_PER_DAY, SECONDS_PER_WEEK, SECONDS_PER_MONTH +from sim_scenes.func import calc_run +from bodies.body import AU +from simulators.calc_simulator import CalcSimulator + + +def compute_barycenter(masses, positions): + """ + Compute the barycenter position of celestial objects in 3D space + masses: a list of masses of celestial objects + positions: a list of positions of celestial objects, each position is a tuple (x, y, z) + """ + m_sum = sum(masses) + x_sum = 0 + y_sum = 0 + z_sum = 0 + for i in range(len(masses)): + x_sum += masses[i] * positions[i][0] / m_sum + y_sum += masses[i] * positions[i][1] / m_sum + z_sum += masses[i] * positions[i][2] / m_sum + return (x_sum, y_sum, z_sum) + + +def get_lagrangian_points(m1, m2, r): + """ + https://baike.baidu.com/item/%E6%8B%89%E6%A0%BC%E6%9C%97%E6%97%A5%E7%82%B9/731078 + + @param m1: 大质量 + @param m2: 小质量 + @param r: 半径 + @return: + """ + a = m2 / (m1 + m2) + l1 = (0, 0, r * (1 - pow(a / 3, 1 / 3))) + l2 = (0, 0, r * (1 + pow(a / 3, 1 / 3))) + l3 = (0, 0, -r * (1 + (5 * a) / 12)) + l4 = (pow(3, 1 / 2) / 2 * r, 0, (r / 2) * ((m1 - m2) / (m1 + m2))) + l5 = (-pow(3, 1 / 2) / 2 * r, 0, (r / 2) * ((m1 - m2) / (m1 + m2))) + + return l1, l2, l3, l4, l5 + + +if __name__ == '__main__': + """ + 地球、月球 + """ + OFFSETTING = 0 + bodies = [ + Earth(init_position=[0, 0, 0], texture="earth_hd.jpg", + init_velocity=[OFFSETTING, 0, 0], size_scale=0.5e1), # 地球放大 5 倍，距离保持不变 + Moon(init_position=[0, 0, 363104], # 距地距离约: 363104 至 405696 km + init_velocity=[-1.054152222, 0, 0], size_scale=1e1) # 月球放大 10 倍，距离保持不变 + ] # -1.0543 < -1.05435 < -1.0545 ？-1.4935 OK：-1.05918 + # -1.05415< <-1.05425 + # 1047.4197364163992 + earth = bodies[0] + moon = bodies[1] + # -1.0648500185012806 -1.05435 + # 1.0544061918995067 3.02326384371554e-06 <月球(Moon)> m=7.342e+22(kg), r|d=1.737e+03|3.474e+03(km), v=2.196e+10(km³), d=3.344e+03(kg/m³), p=[-3.796e+03,0.000e+00,3.631e+05](km), v=[-1.05435 0. -0.01088375](km/s) + # 1.0544608857544382 3.023593683297842e + points = get_lagrangian_points(earth.mass, moon.mass, 363104) + velocities = [] + for i in range(30): + v = round(-0.846 - (i / 10000), 4) + velocities.append([v, 0, 0]) + + satellites = [] + point = points[0] + for j, velocitie in enumerate(velocities): + satellite = Satellite(name=f'卫星{j + 1}', mass=1.4e10, size_scale=1e3, color=(255, 200, 0), + init_position=[point[0], + point[1], + point[2]], + init_velocity=velocities[j]) + # bodies.append(satellite) + satellites.append(satellite) + + CalcSimulator.init_velocity_x = moon.init_velocity[0] + CalcSimulator.offset_rate_x = 0.0001 + CalcSimulator.accelerations = [] + CalcSimulator.velocities = [] + + + def on_init(bodies): + return bodies + + + def on_ready(simulator): + pass + + + def evolve_next(simulator): + if not hasattr(simulator, "loop"): + simulator.loop = 1000 + else: + simulator.loop -= 1 + # print(simulator.bodies_sys.bodies) + moon = simulator.bodies_sys.bodies[1] + + from simulators.ursina.entities.entity_utils import get_value_direction_vectors + + velocity = get_value_direction_vectors(moon.velocity) + acceleration = get_value_direction_vectors(moon.acceleration) + vel_value = velocity[0] # km/s + acc_value = acceleration[0] # km/s² + print(vel_value, acc_value, moon) + + CalcSimulator.accelerations.append(acc_value) + CalcSimulator.velocities.append(vel_value) + + # if not hasattr(simulator, "init_acc_value") and acc_value > 0: + # simulator.init_acc_value = acc_value + # elif hasattr(simulator, "init_acc_value"): + # if simulator.loop == 1: + # diff = acc_value - simulator.init_acc_value + # if CalcSimulator.init_velocity_x > 0: + # d = 1 + # else: + # d = -1 + # if abs(diff) < 1e-10: + # print("完成", CalcSimulator.init_velocity_x) + # exit(0) + # elif diff > 0: + # CalcSimulator.init_velocity_x += CalcSimulator.offset_rate_x * d + # print("慢慢靠近", CalcSimulator.init_velocity_x) + # else: + # CalcSimulator.init_velocity_x -= CalcSimulator.offset_rate_x * d + # print("慢慢远离", CalcSimulator.init_velocity_x) + + return simulator.loop > 0 + + + loop = 1 + while loop > 0: + moon.init_velocity = [CalcSimulator.init_velocity_x, moon.init_velocity[1], moon.init_velocity[2]] + calc_run(bodies, SECONDS_PER_HOUR, + on_init=on_init, + on_ready=on_ready, + evolve_next=evolve_next) + loop -= 1 + + import matplotlib.pyplot as plt + + plt.figure(figsize=(8, 6)) + max_value = max(CalcSimulator.accelerations[1:]) + min_value = min(CalcSimulator.accelerations[1:]) + x = [] + for i in range(len(CalcSimulator.accelerations)): + x.append(i) + plt.title("%s max:%.4f mix:%.4f diff:%.4f" % (moon.init_velocity[0], max_value*1e6, min_value*1e6, (max_value - min_value)*1e6)) + plt.ylim(2.95e-6, 3.06e-6) + plt.plot(x[1:], CalcSimulator.accelerations[1:], label="accelerations") + # plt.plot(x, CalcSimulator.velocities, label="velocities") + + plt.legend() + plt.show() diff --git a/sim_scenes/func.py b/sim_scenes/func.py index c623e13..87eea48 100644 --- a/sim_scenes/func.py +++ b/sim_scenes/func.py @@ -11,13 +11,27 @@ from common.consts import SECONDS_PER_WEEK, SECONDS_PER_MINUTE, SECONDS_PER_HALF from common.func import calculate_distance from common.system import System from bodies import Body -from simulators.ursina.ursina_config import UrsinaConfig -from simulators.ursina.ursina_event import UrsinaEvent from common.consts import LIGHT_SPEED import math import numpy as np +def calc_run(bodies, dt=SECONDS_PER_WEEK, on_init=None, **kwargs): + from simulators.calc_simulator import CalcSimulator + import copy + if on_init is not None: + _bodies = copy.deepcopy(bodies) + _bodies = on_init(_bodies) + if _bodies is None: + _bodies = bodies + else: + _bodies = bodies + + body_sys = System(_bodies) + simulator = CalcSimulator(body_sys) + simulator.run(dt, **kwargs) + + def mayavi_run(bodies, dt=SECONDS_PER_WEEK, view_azimuth=0, view_distance='auto', view_focalpoint='auto', bgcolor=(1 / 255, 1 / 255, 30 / 255)): @@ -147,6 +161,8 @@ def ursina_run(bodies, ursina_view.update() import sys + from simulators.ursina.ursina_config import UrsinaConfig + from simulators.ursina.ursina_event import UrsinaEvent sys.modules["__main__"].update = callback_update if show_trail: UrsinaConfig.show_trail = show_trail @@ -254,6 +270,7 @@ def create_light_ship(size_scale, init_position, speed=LIGHT_SPEED): def create_text_panel(width=0.35, height=.5): # 创建一个 Panel 组件 from ursina import Text, Panel, color, camera, Vec3 + from simulators.ursina.ursina_config import UrsinaConfig panel = Panel( parent=None, model='quad', diff --git a/simulators/calc_simulator.py b/simulators/calc_simulator.py new file mode 100644 index 0000000..3075968 --- /dev/null +++ b/simulators/calc_simulator.py @@ -0,0 +1,139 @@ +# -*- coding:utf-8 -*- +# title :计算运行模拟器（无界面） +# description :计算运行模拟器（无界面） +# author :Python超人 +# date :2023-02-11 +# link :https://gitcode.net/pythoncr/ +# python_version :3.8 +# mayavi version :4.8.1 +# ============================================================================== +from mayavi import mlab +from simulators.simulator import Simulator +from common.system import System +from simulators.views.body_view import BodyView +from common.singleton import Singleton + + +class CalcView(BodyView): + """ + 无界面 + """ + + def update(self): + pass + + def appear(self): + if hasattr(self.body, "torus_stars"): + # 暂不支持环状小行星群 + return + + def disappear(self): + pass + + +class CalcSimulator(Simulator): + """ + 计算运行模拟器（无界面） + 主要用于天体测试数据计算 + """ + + def __init__(self, bodies_sys: System): + super().__init__(bodies_sys, CalcView) + + def run(self, dt, **kwargs): + on_finished = None + if "on_finished" in kwargs: + on_finished = kwargs["on_finished"] + + on_ready = None + if "on_ready" in kwargs: + on_ready = kwargs["on_ready"] + + evolve_next = None + if "evolve_next" in kwargs: + evolve_next = kwargs["evolve_next"] + + after_evolve = None + if "after_evolve" in kwargs: + after_evolve = kwargs["after_evolve"] + + before_evolve = None + if "before_evolve" in kwargs: + before_evolve = kwargs["before_evolve"] + + if on_ready is not None: + on_ready(self) + + if evolve_next is None: + if before_evolve is not None: + before_evolve(self) + + self.evolve(dt) + + if after_evolve is not None: + after_evolve(self) + else: + while evolve_next(self): + if before_evolve is not None: + before_evolve(self) + + self.evolve(dt) + + if after_evolve is not None: + after_evolve(self) + + if on_finished is not None: + on_finished(self) + + +if __name__ == '__main__': + from sim_scenes.func import calc_run + from bodies import Sun, Earth + from common.consts import SECONDS_PER_WEEK + + """ + 3个太阳、1个地球 + """ + bodies = [ + Sun(name='太阳1', mass=1.5e30, init_position=[849597870.700, 0, 0], init_velocity=[0, 7.0, 0], + size_scale=5e1, texture="sun1.jpg"), # 太阳放大 100 倍 + Sun(name='太阳2', mass=2e30, init_position=[0, 0, 0], init_velocity=[0, -8.0, 0], + size_scale=5e1, texture="sun2.jpg"), # 太阳放大 100 倍 + Sun(name='太阳3', mass=2.5e30, init_position=[0, -849597870.700, 0], init_velocity=[18.0, 0, 0], + size_scale=5e1, texture="sun2.jpg"), # 太阳放大 100 倍 + Earth(name='地球', init_position=[0, -349597870.700, 0], init_velocity=[15.50, 0, 0], + size_scale=4e3, distance_scale=1), # 地球放大 4000 倍，距离保持不变 + ] + + + def on_finished(simulator): + print(simulator) + + + def on_ready(simulator): + print(simulator) + + + def after_evolve(simulator): + print(simulator) + + + def before_evolve(simulator): + print(simulator) + + + def evolve_next(simulator): + if not hasattr(simulator, "loop"): + simulator.loop = 10 + else: + simulator.loop -= 1 + print(simulator.bodies_sys.bodies) + return simulator.loop > 0 + + + calc_run(bodies, SECONDS_PER_WEEK, + on_ready=on_ready, + on_finished=on_finished, + after_evolve=after_evolve, + before_evolve=before_evolve, + evolve_next=evolve_next) -- GitLab