为什么我的行星轨道模拟产生的是直线而不是椭圆轨道？

Question

我正在尝试使用

compute_orbit

方法模拟行星的轨道，但是当我绘制结果位置时，我得到的是一条直线而不是预期的椭圆轨道。以下是我的代码的相关部分。

代码片段

def get_initial_conditions(self, planet_name):
    planet_id = self.PLANETS[planet_name]
    obj = Horizons(id=planet_id, location='@sun', epochs=2000.0)
    eph = obj.vectors()
    position = np.array([eph['x'][0], eph['y'][0], eph['z'][0]]) 
    velocity = np.array([eph['vx'][0], eph['vy'][0], eph['vz'][0]])   
    scale_factor_position = 1  
    scale_factor_velocity = 1  
    return {
        "position": np.array(position) / scale_factor_position,
        "velocity": np.array(velocity) / scale_factor_velocity
    }

def compute_orbit(self, central_mass=1.989e30, dt=10000, total_time=31536000):
    initial_conditions_data = self.get_initial_conditions(self.planet_name)
    initial_conditions = [
        initial_conditions_data['position'][0], initial_conditions_data['velocity'][0],
        initial_conditions_data['position'][1], initial_conditions_data['velocity'][1],
        initial_conditions_data['position'][2], initial_conditions_data['velocity'][2]
    ]

    def f(t, state):
        x, vx, y, vy, z, vz = state
        r = np.sqrt(x**2 + y**2 + z**2) + 1e-5  
        G = 6.67430e-11
        Fx = -G * central_mass * self.mass * x / r**3
        Fy = -G * central_mass * self.mass * y / r**3
        Fz = -G * central_mass * self.mass * z / r**3
        return [vx, Fx / self.mass, vy, Fy / self.mass, vz, Fz / self.mass]

    t_span = (0, total_time)
    t_eval = np.arange(0, total_time, dt)  
    sol = scipy.integrate.solve_ivp(
        f, 
        t_span, 
        initial_conditions, 
        t_eval=t_eval, 
        rtol=1e-3, 
        atol=1e-6
    )
    
    x = sol.y[0]  
    y = sol.y[2]  
    z = sol.y[4]  
    positions = np.column_stack((x, y, z))
    return positions

我尝试过使用不同的初始条件进行实验，但我总是得到绘制数据的直线，下面是一个最小的可重现示例

代码片段

from astroquery.jplhorizons import Horizons
import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def get_initial_conditions(planet_id):
        obj = Horizons(id=planet_id, location='@sun', epochs=2000.0)
        eph = obj.vectors()
        position = np.array([eph['x'][0], eph['y'][0], eph['z'][0]]) 
        velocity = np.array([eph['vx'][0], eph['vy'][0], eph['vz'][0]])   
        scale_factor_position = 1  
        scale_factor_velocity = 1  
        return {
            "position": np.array(position) / scale_factor_position,
            "velocity": np.array(velocity) / scale_factor_velocity
        }

def compute_orbit(central_mass=1.989e30,rotating_mass=5.972e24, dt=10000, total_time=31536000):
        initial_conditions_data = get_initial_conditions(399)
        initial_conditions = [
            initial_conditions_data['position'][0], initial_conditions_data['velocity'][0],
            initial_conditions_data['position'][1], initial_conditions_data['velocity'][1],
            initial_conditions_data['position'][2], initial_conditions_data['velocity'][2]
                            ]

        def f(t, state):
            x, vx, y, vy, z, vz = state
            r = np.sqrt(x**2 + y**2 + z**2) + 1e-5  
            G = 6.67430e-11
            Fx = -G * central_mass * rotating_mass * x / r**3
            Fy = -G * central_mass * rotating_mass * y / r**3
            Fz = -G * central_mass * rotating_mass * z / r**3
            return [vx, Fx / rotating_mass, vy, Fy / rotating_mass, vz, Fz / rotating_mass]


        t_span = (0, total_time)
        t_eval = np.arange(0, total_time, dt)  
        sol = scipy.integrate.solve_ivp(
        f, 
        t_span, 
        initial_conditions, 
        t_eval=t_eval, 
        rtol=1e-3, 
        atol=1e-6)
        x = sol.y[0]  
        y = sol.y[2]  
        z = sol.y[4]  
        positions = np.column_stack((x, y, z))
        return positions

def plot_orbit(positions):
    x = positions[:, 0]
    y = positions[:, 1]
    z = positions[:, 2]
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(x, y, z, label='Orbit Path', color='blue')
    ax.set_xlabel('X Position (m)')
    ax.set_ylabel('Y Position (m)')
    ax.set_zlabel('Z Position (m)')
    ax.set_title('Planet Orbit Simulation')
    ax.scatter(0, 0, 0, color='yellow', s=100, label='Central Mass (Sun)')
    ax.legend()
    ax.set_box_aspect([1, 1, 1]) 
    plt.show()


positions = compute_orbit()
plot_orbit(positions)

Answer 1

此问题的单位不匹配。

在以下代码中：

    obj = Horizons(id=planet_id, location='@sun', epochs=2000.0)
    eph = obj.vectors()
    position = np.array([eph['x'][0], eph['y'][0], eph['z'][0]]) 
    velocity = np.array([eph['vx'][0], eph['vy'][0], eph['vz'][0]])   

obj.vectors()

查询不返回以米和米每秒为单位的输出。正如文档所述，它返回 AU 和每天 AU 的输出。

列名称	定义
...	...
x	位置向量的 x 分量（float、au、X）
vx	速度矢量的 x 分量（float、au/d、VX）
...	...

来源。

稍后，您将使用仅适用于千克、米和秒单位的引力常数。

G = 6.67430e-11

解决此问题的最简单方法是将 Horizons 使用的单位转换为米和米每秒。或者，您可以将

G

更改为适合 AU 和天数的值。

def get_initial_conditions(planet_id):
        obj = Horizons(id=planet_id, location='@sun', epochs=2000.0)
        eph = obj.vectors()
        meters_per_au = 1.496e+11
        seconds_per_day = 86400
        position = np.array([eph['x'][0], eph['y'][0], eph['z'][0]]) * meters_per_au
        velocity = np.array([eph['vx'][0], eph['vy'][0], eph['vz'][0]]) * meters_per_au / seconds_per_day
        scale_factor_position = 1  
        scale_factor_velocity = 1  
        return {
            "position": np.array(position) / scale_factor_position,
            "velocity": np.array(velocity) / scale_factor_velocity
        }

您应该得到以下结果。