# -*- coding:utf-8 -*- # title : # description : # author :Python超人 # date :2023-11-16 # link :https://gitcode.net/pythoncr/ # python_version :3.9 # ChatGPT # ============================================================================== """ 我把我的详细需求写给你，你帮我分析用什么方法。我要计算哈雷彗星在太阳系中某个位置上面的初始速度（vx, vy, vz） ， 通过万有引力的计算，哈雷彗星会在一个轨道上有2个点，近日点和远日点。哈雷彗星从近日点和远日点花了多长时间。 我现在定义了一个函数，f(vx,vy,vz) 参数就是初始速度（vx, vy, vz）， 这个函数f返回3个值，分别是 近日点距离（单位：AU）、远日点距离（单位：AU）、以及哈雷彗星近日点到远日点花了多长时间（单位：天）。 我现在就是希望如果函数f返回值是 0.586AU, 35.1AU, 13817天，那么 vx,vy,vz的值是多少？ 以上的需求，我用什么算法可以实现。 """ import threading 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 calc_simulator(params, factor=1): from sim_scenes.func import calc_run from simulators.calc_simulator import CalcSimulator, CalcContext CalcSimulator.init(True) sim = HalleyCometSim(params) sim.total_times = 0 sim.build() sim.result = None 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 = current_dt(time_data) elif d_sun > sim.comet_aphel: sim.comet_aphel = d_sun sim.comet_aphel_dt = current_dt(time_data) # log(time_data) # 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 = current_dt(time_data) elif d_sun < sim.comet_peri: sim.comet_peri = d_sun sim.comet_peri_dt = current_dt(time_data) # log(time_data) # print("距离地球: %.3f AU" % (d_earth / AU)) def current_dt(time_data: TimeData): # f'{time_data.years}y{time_data.days}d' return time_data.years, time_data.days, time_data.total_days def log(time_data): print(params) comet_aphel = sim.comet_aphel / AU comet_peri = sim.comet_peri / AU diff_days = sim.comet_aphel_dt[2] - sim.comet_peri_dt[2] print("远日点:%.3f AU, 近日点:%.3f AU, 天数: %i" % (comet_aphel, comet_peri, diff_days)) """ 哈雷彗星 根据最新信息，哈雷彗星的平均运行速度约为每小时 70,000 英里或每小时 126,000 公里 。 该速度表示彗星穿过太阳系的速度。 值得注意的是，截至 2023 年，所提供的信息都是准确的。 质量: 10^15kg 远日点: 35.1 AU（2023年12月9日） 近日点: 0.586 AU 上次通过近日点：1986年2月9日 下次通过近日点：2061年7月28日 平均密度: 0.6（估计的范围在0.2至1.5g/cm） 均密度: 1 g/cm³ -> 1✕10³ kg/m³ 来源：https://www.cgmodel.com/model/500318.html 2023年12月9日 - 1986年2月9日 = 13817 35.1AU 0.586AU """ if abs(comet_peri - 0.586) < 0.01 and abs(comet_aphel - 35.1) < 0.1: # and abs(diff_days - 13817) < 2: print("OK ......", comet_peri, comet_aphel, diff_days) # print("year days:", time_data.years, time_data.days) # print("远日点: %.3f AU(%s)" % (sim.comet_aphel / AU, sim.comet_aphel_dt)) # print("近日点: %.3f AU(%s)" % (sim.comet_peri / AU, sim.comet_peri_dt)) # # 2023年12月9日 - 1986年2月9日 = 13817 # print("远\近日天数: %i 天 " % (sim.comet_aphel_dt[2] - sim.comet_peri_dt[2])) def on_finished(context: CalcContext): global result # print("on_finished", context) time_data = context.get_param("time_data") log(time_data) comet_aphel = sim.comet_aphel / AU comet_peri = sim.comet_peri / AU diff_days = sim.comet_aphel_dt[2] - sim.comet_peri_dt[2] sim.result = (comet_aphel, comet_peri, float(diff_days)) 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 * factor, evolve_next=evolve_next) return sim.result def target_function(x, y, z, factor): params = HalleyCometParams( start_time='1982-09-24 00:00:00', # init_velocity=[-2.836 + (x / 1000), 4.705 + (y / 1000), 8.85 + (z / 1000)], init_velocity=[x, y, z], init_position=[0, -5 * AU, -10 * AU] ) return calc_simulator(params, factor) if __name__ == '__main__': pass # 近日点 0.586 AU # 上次通过近日点： 1986年2月9日 # 下次通过近日点： 2061年7月28日 # 远日点 35.1 AU 2023年12月9日 # [3.34, 0, 10.7] 2060-4- # [3.34, 0, 10.712] 2061-5 # [3.34, 0, 10.715] 2061-6-24 # [3.34, 0, 10.718] 2061-8 # init_velocity=[3.34, 0, 10.718], # init_position=[0, 0.5 * AU, -10 * AU] # init_velocity=[-3.34, 3, 10.718], # init_position=[0, -2 * AU, -10 * AU]) \ # start_time='1982-09-24 00:00:00', # init_velocity=[-2.836, 4.705, 8.85], # init_position=[0, -5 * AU, -10 * AU] # 34.801 AU(2019-06-10) # 35.086 AU(2023-09-01) # 0.586 AU(1986-02-09) # start_time='1982-09-24 00:00:00', # init_velocity=[-2.837, 4.71, 8.852], # init_position=[0, -5 * AU, -10 * AU] # 34.840AU34.840 AU(2019-05-10) # 35.149 AU(2023-10-08) # 0.586 AU(1986-02-09) # 35.198 AU(2023-11-07) # 0.588 AU(1986-02-09) # start_time='1982-09-24 00:00:00', # init_velocity=[-2.841, 4.7, 8.86], # init_position=[0, -5 * AU, -10 * AU] # target_function(-2.836, 4.705, 8.85) # target_function(-2.816, 4.725, 8.86) # 远日点:35.276 AU, 近日点:0.576 AU, 天数: 13821 # target_function(-2.816, 4.725, 8.855) # 远日点:35.212 AU, 近日点:0.576 AU, 天数: 13785 # target_function(-2.816, 4.725, 8.857) # 远日点:35.237 AU, 近日点:0.576 AU, 天数: 13797 # target_function(-2.816, 4.725, 8.858) # 远日点:35.251 AU, 近日点:0.576 AU, 天数: 13809 # target_function(-2.81, 4.725, 8.858) # 远日点:35.230 AU, 近日点:0.573 AU, 天数: 13791 # target_function(-2.81, 4.735, 8.858) # 远日点:35.297 AU, 近日点:0.573 AU, 天数: 13834 # target_function(-2.81, 4.73, 8.858) # 远日点:35.264 AU, 近日点:0.573 AU, 天数: 13815 # target_function(-2.820, 4.73, 8.858) # 远日点:35.298 AU, 近日点:0.578 AU, 天数: 13834 # target_function(-2.820, 4.73, 8.85) # 远日点:35.196 AU, 近日点:0.578 AU, 天数: 13779 # target_function(-2.825, 4.73, 8.85) # 远日点:35.215 AU, 近日点:0.581 AU, 天数: 13791 # target_function(-2.825, 4.73, 8.848) # 远日点:35.190 AU, 近日点:0.581 AU, 天数: 13773 # x, y, z => a, b, c # x=a,b,c # y=a,c # z=a,c # target_function(-2.835, 4.72, 8.847) # 远日点:35.144 AU, 近日点:0.586 AU, 天数: 13748 # target_function(-2.835, 4.71, 8.86) # 远日点:35.243 AU, 近日点:0.584 AU, 天数: 13809 # target_function(-2.835, 4.73, 8.85) # OK 远日点:35.251 AU, 近日点:0.585 AU, 天数: 13815 # 2023年12月9日 - 1986年2月9日 = 13817 # 35.1AU 0.586AU a, b, c = target_function(-2.774, 5.126, 8.65, 1 / 8) # (-2.774, 5.126, 8.65, 1 / 8) 远日点:35.143 AU, 近日点:0.586 AU, 天数: 13753 # (-2.775, 5.126, 8.65, 1 / 8) 远日点:35.146 AU, 近日点:0.587 AU, 天数: 13755 # (-2.775, 5.126, 8.648, 1 / 8) 远日点:35.122 AU, 近日点:0.587 AU, 天数: 13743 # (-2.775, 5.125, 8.65, 1 / 8) OK 远日点:35.137 AU, 近日点:0.587 AU, 天数: 13751 # (-2.775, 5.12, 8.65, 1 / 8) OK 远日点:35.101 AU, 近日点:0.586 AU, 天数: 13728 # (-2.774, 5.12, 8.66, 1 / 8) 远日点:35.222(124) AU, 近日点:0.585 AU, 天数: 13803(77) # (-2.774, 5.13, 8.65, 1 / 8) 远日点:35.169(71) AU, 近日点:0.587 AU, 天数: 13770(51) # (-2.774, 5.12, 8.65, 1 / 8) OK 远日点:35.098 AU, 近日点:0.586 AU, 天数: 13726 # (-2.775, 5.13, 8.65, 1 / 8) 远日点:35.175 AU, 近日点:0.587 AU, 天数: 13773 # (-2.778, 5.13, 8.65, 1 / 8) OK 远日点:35.184 AU, 近日点:0.588 AU, 天数: 13779 # (-2.77, 5.12, 8.65, 1 / 8) OK 远日点:35.083 AU, 近日点:0.584 AU, 天数: 13714 # (-2.78, 5.11, 8.65, 1 / 8) OK 远日点:35.046 AU, 近日点:0.587 AU, 天数: 13696 # (-2.80, 5.10, 8.65, 1 / 8) OK 远日点:35.044 AU, 近日点:0.595 AU, 天数: 13698 # (-2.47, 5.13, 8.73, 1 / 12) 远日点:35.155 AU, 近日点:0.468 AU, 天数: 13690 # [-2.5, 5.02, 8.763] 远日点:35.032 AU, 近日点:0.477 AU, 天数: 13615 # [-2.835, 4.81, 8.8] 远日点:35.161 AU, 近日点:0.591 AU, 天数: 13754 # [-2.835, 4.8, 8.81] 远日点:35.220 AU, 近日点:0.590 AU, 天数: 13797 # [-2.835, 4.8, 8.8], 远日点:35.092 AU, 近日点:0.590 AU, 天数: 13718 # [-2.835, 4.83, 8.8] 远日点:35.299 AU, 近日点:0.592 AU, 天数: 13846 # [-2.835, 4.75, 8.83] 远日点:35.132 AU, 近日点:0.588 AU, 天数: 13742 # [-2.835, 4.74, 8.84] 远日点:35.191 AU, 近日点:0.587 AU, 天数: 13779 # [-2.835, 4.79, 8.8] 远日点:35.025 AU, 近日点:0.589 AU, 天数: 13682 # 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] # ) # year days: 42 6 # 远日点: 35.086 AU((40, 346, 14946.8522)) # 近日点: 0.586 AU((3, 133, 1228.3418)) # 远\近日天数: 13718 天 # R = 2 # # X, Y, Z = -2.835, 4.73, 8.85 # X, Y, Z = -2.835, 4.72, 8.847 # X, Y, Z = -2.826, 4.695, 8.86 # idx = 0 # for x in range(-R, R + 1): # for y in range(-R, R + 1): # for z in range(-R, R + 1): # x, y, z = X + (x / 100), Y + (y / 100), Z + (z / 100) # idx += 1 # print(f"Index:{idx} ---------------------------") # a, b, c = target_function(x, y, z) # def t(): # global idx # params = HalleyCometParams( # start_time='1982-09-24 00:00:00', # init_velocity=[-2.836 + (x / 100), 4.705 + (y / 100), 8.85 + (z / 100)], # init_position=[0, -5 * AU, -10 * AU] # ) # idx += 1 # print(f"Index:{idx} ---------------------------") # calc_simulator(params) # # threading.Thread(target=t).start() # t()