[Python] Why does my planetary orbit simulation produce a straight line instead of an...

I'm trying to simulate the orbit of a planet using the compute_orbit method, but when I plot the resulting positions, I get a straight line instead of an expected elliptical orbit. Below are the relevant parts of my code.

Code Snippet


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


i have tried experimenting with different initial conditions but i always get straight lines plotting the data , below is a minimal reproduceable example

Code Snippet


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)

Continue reading...