diff --git "a/docs/\346\213\211\346\240\274\346\234\227\346\227\245\347\202\271\350\256\241\347\256\227\345\205\254\345\274\217.md" "b/docs/\346\213\211\346\240\274\346\234\227\346\227\245\347\202\271\350\256\241\347\256\227\345\205\254\345\274\217.md" new file mode 100644 index 0000000000000000000000000000000000000000..baa42c331c78c762af21210ec403f0df63249942 --- /dev/null +++ "b/docs/\346\213\211\346\240\274\346\234\227\346\227\245\347\202\271\350\256\241\347\256\227\345\205\254\345\274\217.md" @@ -0,0 +1,39 @@ +### 拉格朗日点公式 + +$$ +a = \frac{m2}{m1+m2} +$$ +``` +a = m2 / (m1 + m2) +``` + +$$ +L1 = (R ⋅ (1 - \sqrt[3]{\frac{a}{3}} ) , 0) +$$ +``` +l1 = (r * (1 - pow(a / 3, 1 / 3)), 0) +``` +$$ +L2 = (R ⋅ (1 + \sqrt[3]{\frac{a}{3}} ) , 0) +$$ +``` +l2 = (r * (1 + pow(a / 3, 1 / 3)), 0) +``` +$$ +L3 = (-R ⋅ (1 + {\frac{5 ⋅ a}{12}} ) , 0) +$$ +``` +l3 = (-r * (1 + (5 * a) / 12), 0) +``` +$$ +L4 = (\frac{R}{2} ⋅ {\frac{m1-m2}{m1+m2}} , \frac{\sqrt{3}}{2} ⋅R) +$$ +``` +l4 = ((r / 2) * ((m1 - m2) / (m1 + m2)), pow(3, 1 / 2) / 2 * r) +``` +$$ +L5 = (\frac{R}{2} ⋅ {\frac{m1-m2}{m1+m2}} , -\frac{\sqrt{3}}{2} ⋅R) +$$ +``` +l5 = ((r / 2) * ((m1 - m2) / (m1 + m2)), -pow(3, 1 / 2) / 2 * r) +``` \ No newline at end of file diff --git a/objs/obj.py b/objs/obj.py index e41ca3c3a114d38560571ba664714ba966e0dec4..5a05dcf6dace2a60522ce35711b93a3fc680b116 100644 --- a/objs/obj.py +++ b/objs/obj.py @@ -23,7 +23,7 @@ class Obj(metaclass=ABCMeta): 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, - parent=None, ignore_mass=False, + parent=None, ignore_mass=False, trail_scale_factor=1.0, trail_color=None, show_name=False, rotation=None, gravity_only_for=[], model=None): @@ -95,6 +95,7 @@ class Obj(metaclass=ABCMeta): self.light_disable = False self.__has_rings = False + self.trail_scale_factor = trail_scale_factor def find_model(self, model: str): if not model.endswith(".obj"): @@ -317,7 +318,8 @@ class Obj(metaclass=ABCMeta): def __repr__(self): return '<%s(%s):%s(%s)> m=%.3e(kg), d=%.3e(kg/m³), p=[%.3e,%.3e,%.3e](km), v=%s(km/s), s=%.3e' % \ - (self.name, self.__class__.__name__, os.path.basename(self.model), os.path.basename(self.texture), self.mass, self.density, + (self.name, self.__class__.__name__, os.path.basename(self.model), os.path.basename(self.texture), + self.mass, self.density, self.position[0], self.position[1], self.position[2], self.velocity, self.size_scale) def ignore_gravity_with(self, body): diff --git a/objs/satellite.py b/objs/satellite.py index 80a3f810e9a67bff8aea88838acc1b6f581544f5..d3e81623c925e577ce41faff683fd4adbd941519 100644 --- a/objs/satellite.py +++ b/objs/satellite.py @@ -20,7 +20,7 @@ class Satellite(Obj): texture="satellite.png", size_scale=1.0, distance_scale=1.0, ignore_mass=False, density=1e3, color=(7, 0, 162), trail_color=(255, 255, 255), show_name=False, - model="satellite.obj", + trail_scale_factor=5.0, model="satellite.obj", parent=None, gravity_only_for=[]): params = { "name": name, @@ -34,6 +34,7 @@ class Satellite(Obj): "distance_scale": distance_scale, "ignore_mass": ignore_mass, "trail_color": trail_color, + "trail_scale_factor": trail_scale_factor, "show_name": show_name, "parent": parent, "gravity_only_for": gravity_only_for, diff --git a/sim_lab/lagrangian_points.py b/sim_lab/lagrangian_points.py index 16c4eee606320a4e887f8b805a8dc3d9e02fef9c..2d2259d9b4018e7fc6e0c00208cca0d23b9d0dc7 100644 --- a/sim_lab/lagrangian_points.py +++ b/sim_lab/lagrangian_points.py @@ -1,80 +1,211 @@ -# https://www.163.com/dy/article/G5F1016F053102ZV.html -# https://www.sciencedirect.com/topics/physics-and-astronomy/lagrangian-points - - -# 以下是太阳和地球的第一、二、三个拉格朗日点的真实坐标和速度数据： +""" +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 +""" # -# L1点： 坐标： x = 0.010205 AU， y = 0 AU， z = 0 AU 速度： vx = 0 m/s， vy = 246.593 m/s， vz = 0 m/s +# AU = 1.496e8 # -# 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 +# 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) # -# 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 +# 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 # -# 其中，AU表示“天文单位”，即地球与太阳之间的平均距离，约为1.496 x 10^8公里。速度以m/s为单位。 -# 这些数据是基于2021年3月的真实数据，并经过了卫星组织、NASA和欧洲航天局等机构的验证。 - -# AU = 1.496e8 -AU = 1 - -import matplotlib.pyplot as plt - -points = [(0.010205 * AU, 0), (-0.010205 * AU, 0), - # (0.990445 * AU, 0), - (0.500 * AU, 0.866025 * AU), - (0.500 * AU, -0.866025 * AU)] - -# plt.plot(AU, 0, "r.") -plt.plot(0, 0, "b.") -for x,y in points: - plt.plot(x, y, "g.") -# x = [1, 2, 3, 4] -# y = [2, 4, 2, 6] -# y1 = [e + 1 for e in y] -# y2 = [e + 2 for e in y] -# plt.plot(x, y, "b.") # b：蓝色，.：点 -# plt.plot(x, y1, "ro") # r：红色，o：圆圈 -# plt.plot(x, y2, "kx") # k：黑色，x：x字符(小叉) -plt.show() # 在窗口显示该图片 - +# @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 # -# import matplotlib.pyplot as plt # -# # Data from NASA -# L1 = [1.5e6, 0] -# L2 = [-1.5e6, 0] -# L3 = [-1.5e6, 0] -# L4 = [0.5e6, 0.87e6] -# L5 = [0.5e6, -0.87e6] -# sun = [0, 0] -# earth = [0, -1.5e8] +# def show_figure(points, p1_name, p2_name, unit, barycenter=None): +# import matplotlib.pyplot as plt # -# # Create plot and set axis limits -# fig, ax = plt.subplots() -# # ax.set_xlim(-2.5e8, 2.5e8) -# # ax.set_ylim(-2.5e8, 2.5e8) +# 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 # -# # Plot positions of Lagrange points, Sun and Earth -# ax.plot(sun[0], sun[1], 'o', markersize=10, color='yellow') -# # ax.plot(earth[0], earth[1], 'o', markersize=5, color='blue') -# ax.plot(L1[0], L1[1], 'x', markersize=10, color='red') -# ax.plot(L2[0], L2[1], 'x', markersize=10, color='green') -# ax.plot(L3[0], L3[1], 'x', markersize=10, color='purple') -# ax.plot(L4[0], L4[1], '+', markersize=10, color='red') -# ax.plot(L5[0], L5[1], '+', markersize=10, color='green') +# 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 # -# # Plot labels for Lagrange points, Sun and Earth -# # ax.annotate('Sun', (sun[0]+2e7, sun[1]+2e7)) -# # ax.annotate('Earth', (earth[0]+2e7, earth[1]+2e7)) -# ax.annotate('L1', (L1[0]+2e7, L1[1]+2e7)) -# ax.annotate('L2', (L2[0]+2e7, L2[1]+2e7)) -# ax.annotate('L3', (L3[0]+2e7, L3[1]+2e7)) -# ax.annotate('L4', (L4[0]+2e7, L4[1]+2e7)) -# ax.annotate('L5', (L5[0]+2e7, L5[1]+2e7)) +# if barycenter is not None: +# plt.plot(barycenter[0], barycenter[1], "gx", markersize=10, label=f"L{idx}") # -# # Set title and show plot -# plt.title('Positions of Lagrange Points L1 to L5, Sun and Earth') -# plt.show() \ No newline at end of file +# # 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 ursina_run, camera_look_at +from bodies.body import AU +from simulators.ursina.entities.body_timer import TimeData +from simulators.ursina.entities.entity_utils import create_directional_light +from simulators.ursina.ursina_event import UrsinaEvent +from simulators.ursina.ursina_mesh import create_line + + +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 + # TODO: 可以抵消月球带动地球的力，保持地球在原地 + # OFFSETTING = 0.01265 + 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.03, 0, 0], size_scale=1e1) # 月球放大 10 倍，距离保持不变 + ] + earth = bodies[0] + moon = bodies[1] + + points = get_lagrangian_points(earth.mass, moon.mass, 363104) + offset_points = [ + [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], + [-1.265, 0, 0], + [1.03, 0, 0], + [0, 0, 0], + [0, 0, 0], + ] + for i, point in enumerate(points): + satellite = Satellite(name=f'卫星{i + 1}', mass=1.4e10, size_scale=1e3, color=(255, 200, 0), + init_position=[point[0] + offset_points[i][0], + point[1] + offset_points[i][1], + point[2] + offset_points[i][2]], + init_velocity=velocities[i]) + bodies.append(satellite) + + + def on_ready(): + # 运行前触发 + # 运行开始前，将摄像机指向地球 + + # 摄像机看向地球 + camera_look_at(moon) + + + def on_timer_changed(time_data: TimeData): + from ursina import destroy + if hasattr(earth, "line"): + destroy(earth.line) + earth.line = create_line(from_pos=earth.planet.position, to_pos=moon.planet.main_entity.position) + + + # 订阅事件后，上面的函数功能才会起作用 + # 运行前会触发 on_ready + UrsinaEvent.on_ready_subscription(on_ready) + # 运行中，每时每刻都会触发 on_timer_changed + UrsinaEvent.on_timer_changed_subscription(on_timer_changed) + + # 使用 ursina 查看的运行效果 + # 常用快捷键： P：运行和暂停 O：重新开始 I：显示天体轨迹 + # position = 左-右+、上+下-、前+后- + ursina_run(bodies, SECONDS_PER_HOUR, + position=(-300000, 1500000, -100), + show_timer=True, + show_trail=True) diff --git a/simulators/ursina/entities/entity_utils.py b/simulators/ursina/entities/entity_utils.py index e9e7706224866cd5034bc6c5779cd4a8d1c96254..ccd6f4143a7a1dae78b2308c4c75da106758d8e5 100644 --- a/simulators/ursina/entities/entity_utils.py +++ b/simulators/ursina/entities/entity_utils.py @@ -57,8 +57,16 @@ def trail_init(parent, scale): trail_color = conv_to_vec4_color(parent.body_view.body.trail_color) trail_color = adjust_brightness(trail_color, 0.4) parent.trail_color = color.rgba(trail_color[0], trail_color[1], trail_color[2], 0.6) - # 拖尾球体的大小为该天体的 1/5 - parent.trail_scale = scale / 5 + + if hasattr(parent.body, "trail_scale_factor"): + if parent.body.trail_scale_factor is not None: + parent.trail_scale = parent.body.trail_scale_factor * scale + else: + # 拖尾球体的大小为该天体的 1/5 + parent.trail_scale = scale / 5 + else: + # 拖尾球体的大小为该天体的 1/5 + parent.trail_scale = scale / 5 if parent.trail_scale < 1: # 如果太小，则 pass