Python超人-宇宙模拟器

sim_lab/halley_comet_research_calc.py

+# -*- coding:utf-8 -*-
+# title           :
+# description     :
+# author          :Python超人
+# date            :2023-11-16
+# link            :https://gitcode.net/pythoncr/
+# python_version  :3.9
+# ==============================================================================
+from bodies import Body, Sun, Earth, Moon
+from objs import Obj, Satellite, Satellite2
+from common.consts import SECONDS_PER_HOUR, SECONDS_PER_DAY, SECONDS_PER_WEEK, SECONDS_PER_MONTH
+from bodies.body import AU
+from ursina import camera, application
+
+from common.celestial_data_service import init_bodies_reality_pos_vels, conv_to_astropy_time, \
+    set_solar_system_celestial_position
+from common.consts import SECONDS_PER_YEAR, AU
+from common.func import calculate_distance
+from objs import HalleComet
+from sim_scenes.func import create_text_panel
+from sim_scenes.func import ursina_run, create_sphere_sky
+from sim_scenes.solar_system.halley_comet_lib import HalleyCometSimBase, HalleyCometParams, \
+    create_halley_comet, create_orbit_line
+from simulators.ursina.entities.body_timer import TimeData, BodyTimer
+from simulators.ursina.entities.entity_utils import get_value_direction_vectors
+from simulators.ursina.ui.control_ui import ControlUI
+from simulators.ursina.ursina_config import UrsinaConfig
+from simulators.ursina.ursina_event import UrsinaEvent
+from simulators.ursina.ursina_mesh import create_label
+from simulators.ursina_simulator import UrsinaSimulator
+
+
+class HalleyCometSim(HalleyCometSimBase):
+    """
+    哈雷彗星场景模拟
+    """
+
+    def __init__(self, params=None):
+        super(HalleyCometSim, self).__init__()
+        if params is None:
+            self.params = HalleyCometParams()
+        else:
+            self.params = params
+
+        if isinstance(params.start_time, str):
+            self.start_time = conv_to_astropy_time(params.start_time)
+        else:
+            self.start_time = params.start_time
+
+        # print("北京时间：", dt.to_datetime(timezone=pytz.timezone('Asia/Shanghai')))
+
+    def build(self):
+        dt = self.start_time
+        # dt = None
+        self.build_solar_system(ignore_gravity=True, start_time=dt)
+        # self.bodies = [
+        #     self.sun,  # 太阳
+        #     self.mars,  # 火星
+        #     self.neptune,  # 海王星
+        # ]
+
+        # self.bodies = self.bodies[:1]
+
+        # 创建哈雷彗星创建哈雷彗星
+        self.halley_comet = create_halley_comet(self.params.init_velocity, self.params.init_position)
+        self.bodies.append(self.halley_comet)
+
+    # def on_ready(self):
+    #     """
+    #     事件绑定后，模拟器运行前会触发
+    #     @return:
+    #     """
+    #     # 创建天空
+    #     # UrsinaConfig.trail_type = "line"
+    #     # UrsinaConfig.trail_thickness_factor = 3
+    #     # UrsinaConfig.trail_length = 91
+    #     UrsinaConfig.trail_length = 225
+    #     # camera.clip_plane_near = 0.1
+    #     camera.clip_plane_far = 1000000
+    #     create_sphere_sky(scale=200000)
+    #
+    #     # WorldGrid().draw_axises(10)
+    #
+    #     application.time_scale = 5
+    #     self.orbit_lines = []
+    #     for body in self.bodies[1:]:
+    #         if isinstance(body, HalleComet):
+    #             continue
+    #         orbit_line = create_orbit_line(self.sun, body, self.start_time)
+    #         self.orbit_lines.append(orbit_line)
+    #
+    #     self.text_panel = create_text_panel()
+
+    def show_clock(self, dt):
+        """
+        显示时钟
+        @param dt: 时间 datetime
+        @return:
+        """
+        print(dt.strftime('%Y-%m-%d %H:%M:%S'))
+
+    # def set_bodies_position(self, dt):
+    #     """
+    #     设置天体的位置（包含速度和加速度的信息）
+    #     @param time_data:
+    #     @return:
+    #     """
+    #     # t = self.start_time + time_data.total_days
+    #     set_solar_system_celestial_position(self.bodies, dt, False)
+
+
+def calc_simulator():
+    from sim_scenes.func import calc_run
+    from simulators.calc_simulator import CalcSimulator, CalcContext
+
+    CalcSimulator.init(True)
+
+    params = HalleyCometParams(
+        start_time='1982-09-24 00:00:00',
+        init_velocity=[-2.836, 4.705, 8.85],
+        init_position=[0, -5 * AU, -10 * AU]
+    )
+
+    sim = HalleyCometSim(params)
+    sim.total_times = 0
+    sim.build()
+
+    def on_ready(context: CalcContext):
+        # for body in sim.bodies:
+        #     body.reset()
+        # sim.on_ready()
+        context.init_param("years", 0)
+        pass
+
+    def evolve_next(context: CalcContext):
+        years = context.get_param("years")
+        return years < 42  # context.evolve_count < 1000
+
+    def before_evolve(dt, context: CalcContext):
+        sim.start_time = conv_to_astropy_time(params.start_time)
+        # td = TimeData(sim.total_times * dt, BodyTimer.MIN_UNIT_SECONDS)
+        time_data = TimeData(sim.total_times * dt, BodyTimer.MIN_UNIT_SECONDS)
+        sim.total_times += 1
+
+        t = sim.start_time + time_data.total_days
+        set_solar_system_celestial_position(sim.bodies, t, False)
+
+        _dt = time_data.get_datetime(str(sim.start_time))
+        sim.show_clock(_dt)
+        # sim.set_bodies_position(_dt)
+        # 距离太阳远日点: 35.018 AU(40y310y)
+        # 距离太阳近日点: 0.580 AU(3y138d)
+
+        # 距离太阳远日点: 34.733 AU(40y178y)
+        # 距离太阳近日点: 0.581 AU(3y138d)
+        context.put_param("time_data", time_data)
+
+    def after_evolve(dt, context: CalcContext):
+        d_sun = calculate_distance(sim.halley_comet.position, sim.sun.position)
+        d_earth = calculate_distance(sim.halley_comet.position, sim.earth.position)
+        time_data = context.get_param("time_data")
+        # td = TimeData(sim.total_times * dt * UrsinaSimulator.INTERVAL_FATOR, BodyTimer.MIN_UNIT_SECONDS)
+
+
+        # sim.set_bodies_position(td)
+
+        context.put_param("years", time_data.years)
+
+        # 哈雷彗星离太阳最远的点称为 "aphelion of Halley's Comet"（远日点）
+        if not hasattr(sim, "comet_aphel"):
+            sim.comet_aphel = d_sun
+            sim.comet_aphel_dt = f'{time_data.years}y{time_data.days}d'
+        elif d_sun > sim.comet_aphel:
+            sim.comet_aphel = d_sun
+            sim.comet_aphel_dt = f'{time_data.years}y{time_data.days}d'
+            print("year days:", time_data.years, time_data.days)
+            # sim.comet_aphel_dt = dt.strftime("%Y-%m-%d")
+            # sim.set_bodies_position(td)
+            print("距离太阳远日点: %.3f AU(%s)" % (sim.comet_aphel / AU, sim.comet_aphel_dt))
+            print("距离太阳近日点: %.3f AU(%s)" % (sim.comet_peri / AU, sim.comet_peri_dt))
+            # print("距离地球: %.3f AU" % (d_earth / AU))
+
+        # 哈雷彗星离太阳最近的点称为 "perihelion of Halley's Comet"（近日点：comet_peri），
+        if not hasattr(sim, "comet_peri"):
+            sim.comet_peri = d_sun
+            sim.comet_peri_dt = f'{time_data.years}y{time_data.days}d'
+        elif d_sun < sim.comet_peri:
+            sim.comet_peri = d_sun
+            sim.comet_peri_dt = f'{time_data.years}y{time_data.days}d'
+            print("距离太阳远日点: %.3f AU(%s)" % (sim.comet_aphel / AU,sim.comet_aphel_dt))
+            print("距离太阳近日点: %.3f AU(%s)" % (sim.comet_peri / AU,sim.comet_peri_dt))
+            # print("距离地球: %.3f AU" % (d_earth / AU))
+
+            # sim.comet_peri_dt = dt.strftime("%Y-%m-%d")
+
+        # sim.set_bodies_position(td)
+        # print("year days:", td.years, td.days)
+        # print("距离太阳: %.3f AU" % (d_sun / AU))
+        # print("距离地球: %.3f AU" % (d_earth / AU))
+
+    def on_finished(context: CalcContext):
+        pass
+        print("on_finished", context)
+        # import matplotlib.pyplot as plt
+        # # 解决中文显示问题
+        # plt.rcParams['font.sans-serif'] = ['SimHei']
+        # plt.rcParams['axes.unicode_minus'] = False
+        # acc_values = context.params["acc_values"]
+        # vel_values = context.params["vel_values"]
+        # max_value = max(acc_values)
+        # min_value = min(acc_values)
+        # fig = plt.figure(figsize=(10, 6))
+        # ax1 = fig.add_subplot(111)
+        # ax1.plot(acc_values, 'blue', label='加速度')
+        # ax1.set_ylabel('加速度', fontdict={'weight': 'normal', 'size': 15, 'color': 'blue'})
+        # # ax1.set_title("加速度/速度", fontdict={'weight': 'normal', 'size': 15})
+        # plt.title("%s(%.4f) max:%.4f min:%.4f diff:%.4f" % (target.name,
+        #                                                     context.get_param("init_vel"), max_value * 1e6,
+        #                                                     min_value * 1e6, (max_value - min_value) * 1e6))
+        # l1 = ax1.legend(loc='lower left', labels=['加速度'])
+        #
+        # ax2 = ax1.twinx()  # this is the important function
+        # ax2.plot(vel_values, 'red', label='速度')
+        # ax2.set_ylabel('速度', fontdict={'weight': 'normal', 'size': 15, 'color': 'red'})
+        # ax2.set_xlabel('Same')
+        #
+        # l2 = ax2.legend(loc='upper right', labels=['速度'])
+        # plt.show()
+
+    CalcSimulator.on_ready_subscription(on_ready)
+    CalcSimulator.on_after_evolve_subscription(after_evolve)
+    CalcSimulator.on_before_evolve_subscription(before_evolve)
+    CalcSimulator.on_finished_subscription(on_finished)
+
+    calc_run(sim.bodies, SECONDS_PER_YEAR, evolve_next=evolve_next)
+
+
+def ursina_simulator():
+    from sim_scenes.func import ursina_run, camera_look_at
+    from simulators.ursina.entities.body_timer import TimeData
+    from simulators.ursina.ursina_event import UrsinaEvent
+
+    params = HalleyCometParams(
+        start_time='1982-09-24 00:00:00',
+        init_velocity=[-2.836, 4.705, 8.85],
+        init_position=[0, -5 * AU, -10 * AU]
+    )
+
+    sim = HalleyCometSim(params)
+    sim.build()
+
+    def on_ready():
+        pass
+
+    def on_timer_changed(time_data: TimeData):
+        pass
+
+    # 订阅事件后，上面的函数功能才会起作用
+    # 运行前会触发 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(sim.bodies, SECONDS_PER_YEAR,
+               # position=(-300000, 1500000, -100),
+               position=(0, 0.5 * AU, -5 * AU),
+               show_timer=True,
+               show_trail=True
+               )
+
+
+if __name__ == '__main__':
+    calc_simulator()
+    # calc_simulator(moon)
+    # calc_simulator(satellite)
+    # ursina_simulator()

simulators/calc_simulator.py

@@ -9,6 +9,7 @@
 # ==============================================================================
 from simulators.simulator import Simulator
 from common.system import System
+from simulators.ursina_simulator import UrsinaSimulator
 from simulators.views.body_view import BodyView
 from common.singleton import Singleton
 
@@ -170,6 +171,7 @@ class CalcSimulator(Simulator):
         super().__init__(bodies_sys, CalcView)
 
     def run(self, dt, **kwargs):
+        dt = dt * UrsinaSimulator.INTERVAL_FATOR
         evolve_next = None
         if "evolve_next" in kwargs:
             evolve_next = kwargs["evolve_next"]

simulators/ursina_simulator.py

@@ -34,6 +34,7 @@ class UrsinaSimulator(Simulator):
     """
     Ursina官网： https://www.ursinaengine.org/
     """
+    INTERVAL_FATOR = 0.01
 
     def __init__(self, bodies_sys: System):
         # window.borderless = False
@@ -230,7 +231,7 @@ class UrsinaSimulator(Simulator):
 
             UrsinaEvent.on_evolving(evolve_dt)
             # interval_fator 能让更新天体运行状态（位置、速度）更精确
-            evolve_dt = evolve_dt * self.interval_fator
+            evolve_dt = evolve_dt * self.INTERVAL_FATOR
             if run_speed_factor < 3:
                 evolve_dt *= 1.666  # 人为加入一个针对秒级计算的误差，保证模拟器的1秒和现实同步（不要求精确可以注释掉）
             evolve_args = {"evolve_dt": evolve_dt}
@@ -323,11 +324,10 @@ class UrsinaSimulator(Simulator):
 
         # interval_fator 能让更新天体运行状态（位置、速度）更精确
         # 设定时间间隔为0.01秒
-        self.interval_fator = 0.01
         self.evolve_dt = dt
 
         # interval 和 last_time 用于检查时间间隔是否已过期
-        self.interval = datetime.timedelta(seconds=self.interval_fator)
+        self.interval = datetime.timedelta(seconds=self.INTERVAL_FATOR)
         self.last_time = datetime.datetime.now() - datetime.timedelta(seconds=2)
 
         if "show_grid" in kwargs: