# -*- 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): 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, evolve_next=evolve_next) return sim.result def target_function(x, y, z): 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) if __name__ == '__main__': pass # 2023年12月9日 - 1986年2月9日 = 13817 # 35.1AU 0.586AU # 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 # a, b, c = target_function(-2.47, 5.13, 8.73) # [-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()