Python超人-宇宙模拟器

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画图，并写上详细的注释。
+# """

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

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)

sim_lab/lagrangian_points_3.py

+"""
+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()

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',

simulators/calc_simulator.py

+# -*- 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)