行星在太阳系模拟中不遵循轨道

Question

我正在尝试模拟一个太阳系，但是行星似乎并不遵循轨道，而是行星与太阳之间的距离增加了。我不知道这段代码有什么问题。这是一个 MRE.

import numpy as np
import matplotlib.pyplot as plt


GRAVITATIONAL_CONSTANT = 6.674e-11
EARTH_MASS = 5.972 * (10**24)
SUN_MASS = 332954.355179 * EARTH_MASS  # Relative to earth
MERCURY_MASS = 0.06 * EARTH_MASS
AU2M = 1.495979e11
AUD2MS = AU2M / 86400


class SolarSystemBody:
    def __init__(self, name, mass, position, velocity):
        self.name = name
        self.mass = mass
        self.position = np.array(position, dtype=float) * AU2M
        self.velocity = np.array(velocity, dtype=float) * AUD2MS
        self.acceleration = np.zeros(3, dtype=float)

    def update_position(self, dt):
        self.position += self.velocity * dt + 0.5 * self.acceleration * dt**2

    def update_velocity(self, dt):
        self.velocity += self.acceleration * dt

    def gravitational_force(self, sun):
        r = sun.position - self.position
        distance = np.linalg.norm(r)
        direction = r / distance
        force_magnitude = GRAVITATIONAL_CONSTANT * self.mass * sun.mass / distance**2
        return force_magnitude * direction

    def calculate_acceleration(self, sun):
        force = self.gravitational_force(sun)
        self.acceleration = force / self.mass


mercury_position = [0.1269730114807624, 0.281031132701101, 0.01131924496141172]
mercury_velocity = [-0.03126433724097832, 0.01267637703164289, 0.00390363008183905]

sun = SolarSystemBody("Sun", SUN_MASS, [0, 0, 0], [0, 0, 0])
mercury = SolarSystemBody("Mercury", MERCURY_MASS, mercury_position, mercury_velocity)

dt = 3600 * 24
total_time = 365 * dt

pos = np.zeros((365, 3), dtype=float)
i = 0
for t in np.arange(0, total_time, dt):
    print(np.linalg.norm(sun.position - mercury.position))
    pos[i, :] = mercury.position
    mercury.calculate_acceleration(sun)
    mercury.update_velocity(dt)
    mercury.update_position(dt)
    i += 1


fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot(pos[:, 0], pos[:, 1], pos[:, 2])
plt.show()

使用

astroquery

的HorizonsClass

获取初始数据

from astropy.time import Time
from astroquery.jplhorizons import Horizons
import numpy as np

sim_start_date = "2023-01-01"

data = []
planet_id = 199  # Mercury

obj = Horizons(id=planet_id, location="@sun", epochs=Time(sim_start_date).jd, id_type='id').vectors()
name = obj["targetname"].data[0].split('(')[0].strip()
r = [np.double(obj[xi]) for xi in ['x', 'y', 'z']]
v = [np.double(obj[vxi]) for vxi in ['vx', 'vy', 'vz']]

Answer 1

我认为这是因为水星的位置和速度是以天文单位和天文单位给出的，而引力常数和天体的质量是以国际单位给出的。它可能搞乱了计算。尝试将它们全部转换为 SI 单位。将 AU 转换为米，将速度从 AU/day 转换为 m/s。