# -*- coding:utf-8 -*- # title :天体数据服务类型 # description :天体数据服务类型 # author :Python超人 # date :2023-08-05 # link :https://gitcode.net/pythoncr/ # python_version :3.8 # ============================================================================== import math import numpy as np import pytz from astropy.coordinates import get_body_barycentric_posvel from astropy.time import Time from bodies import Body, Sun, Asteroids, Moon, Earth from common.consts import G, AU, SECONDS_PER_DAY from simulators.ursina.ursina_config import UrsinaConfig def calc_solar_acceleration(body_or_pos, big_body): """ 计算天体的加速度 @param body_or_pos: 需要计算的天体 @param big_body: 大天体（太阳 或者 地球） @return: """ if isinstance(body_or_pos, Body): body_pos = body_or_pos.position else: body_pos = body_or_pos x, y, z = body_pos[0] * 1000 - big_body.position[0] * 1000, \ body_pos[1] * 1000 - big_body.position[1] * 1000, \ body_pos[2] * 1000 - big_body.position[2] * 1000 r = math.sqrt(x ** 2 + y ** 2 + z ** 2) a = G * big_body.mass / r ** 2 ax = -a * x / r ay = -a * y / r az = -a * z / r return [ax / 1000, ay / 1000, az / 1000] # 设置天体的加速度（单位：km/s²） def get_body_posvel(body, time=None): """ 获取太阳系天体指定时间的位置和矢量速度 pip install -i http://pypi.douban.com/simple/ --trusted-host=pypi.douban.com de423 @param body: 天体（天体名称） @param time: 时间 @return: """ from astropy.coordinates import get_body_barycentric_posvel from astropy.time import Time if time is None: time = Time.now() if not isinstance(body, str): body = body.__class__.__name__ posvel = get_body_barycentric_posvel(body, time) return posvel def get_bodies_posvels(planet_names="sun,mercury,venus,earth,moon,mars,jupiter,saturn,uranus,neptune", time=None): from astropy.coordinates import get_body_barycentric_posvel from astropy.time import Time if time is None: time = Time.now() planets = planet_names.split(",") posvels = {} for planet in planets: try: position, velocity = get_body_barycentric_posvel(planet, time) posvels[planet] = position, velocity # print(planet, position) except Exception as e: print(planet, str(e)) return posvels def recalc_moon_position(moon_posvel, earth_pos): """ 重新计算月球的位置（由于月球和地球的距离在整个太阳系尺度下非常近，为了较好的显示效果，需要放大月球和地球的距离，但不要改变月球相对地球的位置） @param moon_posvel: 月球的三维坐标位置和三维矢量速度 @param earth_pos: 地球的三维坐标位置 @return: """ moon_pos, moon_vel = moon_posvel[0], moon_posvel[1] moon_pos_to_earth = moon_pos - earth_pos moon_pos_to_earth = moon_pos_to_earth * 50 return moon_pos_to_earth + earth_pos, moon_vel def get_celestial_body_data(body_name): # pip install ephem import ephem # 创建一个Observer对象，用于指定观测者的位置 observer = ephem.Observer() observer.lat = '0' # 观测者的纬度，这里使用0度作为示例 observer.lon = '0' # 观测者的经度，这里使用0度作为示例 # 创建一个Date对象，表示当前时间 current_time = ephem.now() # 根据天体名称创建一个CelestialBody对象 body = getattr(ephem, body_name)() # 计算天体的位置和速度 body.compute(observer) # 获取天体的实时位置和速度信息 position = (body.ra, body.dec) # 天体的赤经和赤纬 velocity = (body.ra_velocity, body.dec_velocity) # 天体的赤经速度和赤纬速度 return position, velocity def set_solar_system_celestial_position(bodies, dt, recalc_moon_pos, set_velocity=True, set_acceleration=True): """ 根据日期时间 dt 设置太阳系中天体的真实位置 @param bodies: 太阳系中天体 @param dt: 时间 @param recalc_moon_pos: 是否对月球的位置进行重新计算。 为了更好的展示效果，需要对月球的位置重新计算（使得地月距离放大，月球相对地球方向不变）， 重新计算位置后，地球和月球可以放大1000倍以上 @param set_velocity: 是否设置速度 @param set_acceleration: 是否设置加速度 @return: """ earth_pos = None sun_pos = None earth = None sun = None moon = None for body in bodies: if isinstance(body, Sun): sun = body elif isinstance(body, Earth): earth = body elif isinstance(body, Moon): moon = body for body in bodies: if isinstance(body, Asteroids): # 小行星带是模拟，不是正常的天体 posvel = None else: try: # 获取天体的三维位置和矢量速度 posvel = get_body_posvel(body, dt) except Exception: continue if isinstance(body, Moon): # 如果是月球，为了更好的展示效果，需要对月球的位置重新计算 moon_real_pos = [posvel[0].x.value * AU, posvel[0].z.value * AU, posvel[0].y.value * AU] # TODO:注释下行，月球就会在真实的位置 if recalc_moon_pos: posvel = recalc_moon_position(posvel, earth_pos) if posvel is None: # posvel 为空，则使用太阳的坐标 position, velocity = [sun_pos.x.value * AU, sun_pos.z.value * AU, sun_pos.y.value * AU], [0, 0, 0] else: # 坐标单位：千米 速度单位：千米/秒 position, velocity = [posvel[0].x.value * AU, posvel[0].z.value * AU, posvel[0].y.value * AU], \ [posvel[1].x.value * AU / SECONDS_PER_DAY, posvel[1].z.value * AU / SECONDS_PER_DAY, posvel[1].y.value * AU / SECONDS_PER_DAY] # 实时调整天体的位置和速度 body.position = np.array(position) if set_velocity: body.velocity = np.array(velocity) if isinstance(body, Asteroids): pass elif isinstance(body, Sun): # 记录太阳的位置 sun_pos = posvel[0] elif isinstance(body, Moon): if set_acceleration: # 月球受到2个影响比较大的天体引力（地球和太阳），计算引力引起的加速度和 acc_earth = calc_solar_acceleration(moon_real_pos, earth) acc_sun = calc_solar_acceleration(moon_real_pos, sun) body.acceleration = [acc_earth[0] + acc_sun[0], acc_earth[1] + acc_sun[1], acc_earth[2] + acc_sun[2]] # elif isinstance(body, Earth): # # 月球受到2个影响比较大的天体引力（地球和太阳），计算引力引起的加速度和 # acc_earth = calc_solar_acceleration(earth, moon) # acc_sun = calc_solar_acceleration(earth, sun) # body.acceleration = [acc_earth[0] + acc_sun[0], # acc_earth[1] + acc_sun[1], # acc_earth[2] + acc_sun[2]] else: if set_acceleration: # 其他天体受到太阳引力 body.acceleration = calc_solar_acceleration(body, sun) if isinstance(body, Earth): # 记录地球的位置 earth_pos = posvel[0] def set_earth_rotation(earth, dt): """ 根据指定的时间控制地球的旋转角度（保证地球的自转和北京时间同步） @param dt: 时间 datetime @return: """ # timetuple 可以获取当天的小时数、分数钟、秒数 timetuple = dt.timetuple() # 当年的第几天 day_of_year = timetuple.tm_yday # 根据当年的第几天计算出该日期当天的偏转角度：360度 / 365天 = 当天的偏转角度 angle_of_day = day_of_year * (360 / 365) # 计算出精确的小时数 total_hours = timetuple.tm_hour + timetuple.tm_min / 60 + timetuple.tm_sec / 60 / 60 # -total_hours： 负号控制地球的旋转方向、1天24小时，360度/24=15 # total_hours * 15：1天24小时，360度/24小时=1小时15度 # angle_of_day： 1年第几天的角度 earth.planet.rotation_y = -total_hours * 15 - angle_of_day + 15 # 精确调整 def conv_to_astropy_time(str_time): # 安装 astropy 包 # pip install -i http://pypi.douban.com/simple/ --trusted-host=pypi.douban.com astropy from astropy.time import Time # 时间 from datetime import datetime return Time(datetime.strptime(str_time + '+0800', '%Y-%m-%d %H:%M:%S%z'), format='datetime') def get_earth_pos_vel(dt): """ 根据日期时间获取地球的 @param dt: @return: """ import pytz # 安装 astropy 包 # pip install -i http://pypi.douban.com/simple/ --trusted-host=pypi.douban.com astropy # 获取天体的坐标pos和速度vel from astropy.coordinates import get_body_barycentric_posvel import astropy.units as u # 单位 # from astropy.coordinates import solar_system_ephemeris # print("天体名：", solar_system_ephemeris.bodies) # ('earth', 'sun', 'moon', 'mercury', 'venus', 'earth-moon-barycenter', # 'mars', 'jupiter', 'saturn', 'uranus', 'neptune') print('----------------------------------------') print("北京时间：", dt.to_datetime(timezone=pytz.timezone('Asia/Shanghai'))) # 获取地球的坐标和速度 earth_pos_vel = get_body_barycentric_posvel('earth', dt) print("地球坐标(公里)：", [earth_pos_vel[0].x.to(u.km), earth_pos_vel[0].y.to(u.km), earth_pos_vel[0].z.to(u.km)]) print("地球速度(公里/秒)：", [earth_pos_vel[1].x.to(u.km / u.second), earth_pos_vel[1].y.to(u.km / u.second), earth_pos_vel[1].z.to(u.km / u.second)]) print('----------------------------------------') # print("速度(公里/秒)：", posvel[1] * AU / SECONDS_PER_DAY) def show_bodies_posvels(dt): import pytz # 安装 astropy 包 # pip install -i http://pypi.douban.com/simple/ --trusted-host=pypi.douban.com astropy # 获取天体的坐标pos和速度vel from astropy.coordinates import get_body_barycentric_posvel import astropy.units as u # 单位 from astropy.coordinates import solar_system_ephemeris as sse print("北京时间：", dt.to_datetime(timezone=pytz.timezone('Asia/Shanghai'))) for body_name in sse.bodies: earth_pos_vel = get_body_barycentric_posvel(body_name, dt) print(f"{body_name}-坐标(公里)：", [earth_pos_vel[0].x.to(u.km), earth_pos_vel[0].y.to(u.km), earth_pos_vel[0].z.to(u.km)]) print(f"{body_name}-速度(公里/秒)：", [earth_pos_vel[1].x.to(u.km / u.second), earth_pos_vel[1].y.to(u.km / u.second), earth_pos_vel[1].z.to(u.km / u.second)]) def gen_real_pos_vel_bodies(body_names='sun,mercury,venus,earth,mars,jupiter,saturn,uranus,neptune', dt=None): import pytz # 安装 astropy 包 # pip install -i http://pypi.douban.com/simple/ --trusted-host=pypi.douban.com astropy # 获取天体的坐标pos和速度vel from astropy.coordinates import get_body_barycentric_posvel import astropy.units as u # 单位 if dt is None: dt = Time.now() # from astropy.coordinates import solar_system_ephemeris as sse print("北京时间：", dt.to_datetime(timezone=pytz.timezone('Asia/Shanghai'))) gen_body_str = " init_pos_vels = {}\n" for body_name in str(body_names).split(","): earth_pos_vel = get_body_barycentric_posvel(body_name, dt) pos, vel = earth_pos_vel[0], earth_pos_vel[1] gen_body_str += f" init_pos_vels['{body_name.lower()}']={{}}\n" gen_body_str += f" init_pos_vels['{body_name.lower()}']['pos']=[%.4f,%.4f,%.4f]\n" % (pos.x.to(u.km).value, pos.y.to(u.km).value, pos.z.to(u.km).value) gen_body_str += f" init_pos_vels['{body_name.lower()}']['vel']=[%.4f,%.4f,%.4f]\n" % ( vel.x.to(u.km / u.second).value, vel.y.to(u.km / u.second).value, vel.z.to(u.km / u.second).value) return gen_body_str # pos, vel = earth_pos_vel[0], earth_pos_vel[1] # print(f"{body_name}-坐标(公里)：", [pos.x.to(u.km).value, # pos.y.to(u.km).value, # pos.z.to(u.km)].value) # print(f"{body_name}-速度(公里/秒)：", # [vel.x.to(u.km / u.second).value, # vel.y.to(u.km / u.second).value, # vel.z.to(u.km / u.second).value]) def get_init_pos_vels(): """ 获取初始位置和速度（方便模拟） @return: """ init_pos_vels = {} init_pos_vels['sun'] = {} init_pos_vels['sun']['pos'] = [-4638653.0000, 0.0000, 1314332.6250] init_pos_vels['sun']['vel'] = [-0.0165, 0.0000, -0.0097] init_pos_vels['mercury'] = {} init_pos_vels['mercury']['pos'] = [-54412616.0000, 0.0000, 31280100.0000] init_pos_vels['mercury']['vel'] = [-25.2667, 0.0000, -40.0561] init_pos_vels['venus'] = {} init_pos_vels['venus']['pos'] = [-32341672.0000, 0.0000, -99861352.0000] init_pos_vels['venus']['vel'] = [34.3770, 0.0000, -10.6676] init_pos_vels['earth'] = {} init_pos_vels['earth']['pos'] = [140427392.0000, 0.0000, 45724572.0000] init_pos_vels['earth']['vel'] = [-8.1280, 0.0000, 28.2277] init_pos_vels['mars'] = {} init_pos_vels['mars']['pos'] = [218751648.0000, 0.0000, -45263688.0000] init_pos_vels['mars']['vel'] = [5.2205, 0.0000, 23.4601] init_pos_vels['asteroids'] = {} init_pos_vels['asteroids']['pos'] = [0.0000, 0.0000, 0.0000] init_pos_vels['asteroids']['vel'] = [0.0000, 0.0000, 0.0000] init_pos_vels['jupiter'] = {} init_pos_vels['jupiter']['pos'] = [769445760.0000, 0.0000, 2017638.5000] init_pos_vels['jupiter']['vel'] = [0.0389, 0.0000, 13.0879] init_pos_vels['saturn'] = {} init_pos_vels['saturn']['pos'] = [-1330536704.0000, 0.0000, -454165280.0000] init_pos_vels['saturn']['vel'] = [3.2189, 0.0000, -9.1974] init_pos_vels['uranus'] = {} init_pos_vels['uranus']['pos'] = [-1757163648.0000, 0.0000, 2272503808.0000] init_pos_vels['uranus']['vel'] = [-5.3918, 0.0000, -4.1634] init_pos_vels['neptune'] = {} init_pos_vels['neptune']['pos'] = [-1480893824.0000, 0.0000, 4351772160.0000] init_pos_vels['neptune']['vel'] = [-5.1466, 0.0000, -1.7189] init_pos_vels['pluto'] = {} init_pos_vels['pluto']['pos'] = [-1293909504.0000, 0.0000, 5770579968.0000] init_pos_vels['pluto']['vel'] = [-4.5847, 0.0000, -1.0475] return init_pos_vels def get_init_pos_vels_2(): init_pos_vels = {} init_pos_vels['sun'] = {} init_pos_vels['sun']['pos'] = [-8885067.0000, 0.0000, 2311787.7500] init_pos_vels['sun']['vel'] = [-0.0304, 0.0000, -0.0016] init_pos_vels['mercury'] = {} init_pos_vels['mercury']['pos'] = [-62266848.0000, 0.0000, -24917024.0000] init_pos_vels['mercury']['vel'] = [16.2605, 0.0000, -43.2059] init_pos_vels['venus'] = {} init_pos_vels['venus']['pos'] = [72525272.0000, 0.0000, 83115816.0000] init_pos_vels['venus']['vel'] = [-23.0332, 0.0000, 23.5375] init_pos_vels['earth'] = {} init_pos_vels['earth']['pos'] = [93683544.0000, 0.0000, -103641144.0000] init_pos_vels['earth']['vel'] = [22.8265, 0.0000, 19.8364] init_pos_vels['mars'] = {} init_pos_vels['mars']['pos'] = [-74898024.0000, 0.0000, -211774608.0000] init_pos_vels['mars']['vel'] = [23.1957, 0.0000, -7.7939] init_pos_vels['asteroids'] = {} init_pos_vels['asteroids']['pos'] = [0.0000, 0.0000, 0.0000] init_pos_vels['asteroids']['vel'] = [0.0000, 0.0000, 0.0000] init_pos_vels['jupiter'] = {} init_pos_vels['jupiter']['pos'] = [-2776470.5000, 0.0000, -768231616.0000] init_pos_vels['jupiter']['vel'] = [13.1263, 0.0000, 0.0694] init_pos_vels['saturn'] = {} init_pos_vels['saturn']['pos'] = [907136000.0000, 0.0000, -1055790848.0000] init_pos_vels['saturn']['vel'] = [7.4217, 0.0000, 6.3457] init_pos_vels['uranus'] = {} init_pos_vels['uranus']['pos'] = [-2782867968.0000, 0.0000, 714746048.0000] init_pos_vels['uranus']['vel'] = [-1.6953, 0.0000, -6.5996] init_pos_vels['neptune'] = {} init_pos_vels['neptune']['pos'] = [-2812751104.0000, 0.0000, 3651367424.0000] init_pos_vels['neptune']['vel'] = [-4.3237, 0.0000, -3.2580] init_pos_vels['pluto'] = {} init_pos_vels['pluto']['pos'] = [-2529432576.0000, 0.0000, 5336828928.0000] init_pos_vels['pluto']['vel'] = [-4.2408, 0.0000, -2.0510] return init_pos_vels def init_bodies_pos_vels(bodies, init_pos_vels_fun=None): if init_pos_vels_fun is None: init_pos_vels_fun = get_init_pos_vels # 获取模拟的初始位置和速度 init_pos_vels = init_pos_vels_fun() for body in bodies: pos_vels = init_pos_vels.get(type(body).__name__.lower(), None) if pos_vels is None: continue body.init_position = pos_vels['pos'] body.init_velocity = pos_vels['vel'] def get_reality_orbit_points(body_name, start_time=None, days=365, segments=100, scale_factor=1): if start_time is None: start_time = Time.now() if days < segments: s = 1 else: s = int(round(days / segments)) points = [] for d in range(0, days, s): dt = d + start_time pos, vel = get_body_barycentric_posvel(body_name, dt) x, y, z = pos.x.value * AU * scale_factor, pos.z.value * AU * scale_factor, pos.y.value * AU * scale_factor points.append((x, y, z)) return points def init_bodies_reality_pos_vels(bodies, dt=None): from astropy.coordinates import solar_system_ephemeris as sse if dt is None: dt = Time.now() # from astropy.coordinates import solar_system_ephemeris as sse print("北京时间：", dt.to_datetime(timezone=pytz.timezone('Asia/Shanghai'))) gen_body_str = " init_pos_vels = {}\n" # for body_name in str(body_names).split(","): # earth_pos_vel = get_body_barycentric_posvel(body_name, dt) # pos, vel = earth_pos_vel[0], earth_pos_vel[1] # # # 获取模拟的初始位置和速度 # init_pos_vels = gen_real_pos_vel_bodies(dt=dt) for body in bodies: body_name = type(body).__name__.lower() if body_name not in sse.bodies: continue earth_pos_vel = get_body_barycentric_posvel(body_name, dt) if earth_pos_vel is None: continue pos, vel = earth_pos_vel[0], earth_pos_vel[1] # pos_vels = init_pos_vels.get(type(body).__name__.lower(), None) # if pos_vels is None: # continue pos, vel = [pos.x.value * AU, pos.z.value * AU, pos.y.value * AU], \ [vel.x.value * AU / SECONDS_PER_DAY, vel.z.value * AU / SECONDS_PER_DAY, vel.y.value * AU / SECONDS_PER_DAY] body.init_position = pos body.init_velocity = vel if __name__ == '__main__': from astropy.time import Time # 时间 # bodies = [ # sun, # 太阳放大 80 倍 # Mercury(size_scale=4e3, distance_scale=1.3), # 水星放大 4000 倍，距离放大 1.3 倍 # Venus(size_scale=4e3, distance_scale=1.3), # 金星放大 4000 倍，距离放大 1.3 倍 # Earth(size_scale=4e3, distance_scale=1.3), # 地球放大 4000 倍，距离放大 1.3 倍 # Mars(size_scale=4e3, distance_scale=1.2), # 火星放大 4000 倍，距离放大 1.2 倍 # Asteroids(size_scale=1e2, parent=sun), # 小行星模拟(仅 ursina 模拟器支持) # Jupiter(size_scale=0.68e3, distance_scale=0.72), # 木星放大 680 倍，距离缩小到真实距离的 0.72 # Saturn(size_scale=0.68e3, distance_scale=0.52), # 土星放大 680 倍，距离缩小到真实距离的 0.52 # Uranus(size_scale=0.8e3, distance_scale=0.36), # 天王星放大 800 倍，距离缩小到真实距离的 0.36 # Neptune(size_scale=1e3, distance_scale=0.27), # 海王星放大 1000 倍，距离缩小到真实距离的 0.27 # Pluto(size_scale=10e3, distance_scale=0.23), # 冥王星放大 10000 倍，距离缩小到真实距离的 0.23(从太阳系的行星中排除) # ] # 获取当前时间 # dt = Time.now() # # 指定未来的日期时间 # dt = conv_to_astropy_time('2050-01-01 12:00:00') # # show_bodies_posvels(dt) # bs = gen_real_pos_vel_bodies('sun,Mercury,Venus,Earth,Mars,Jupiter,Saturn,Uranus,Neptune',dt) # print(bs) print(get_reality_orbit_points('earth')) # print(get_init_pos_vels()['sun'])