Question: 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) 

Answer: This problem has mismatched units. In the following code:  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]])  The obj.vectors() query does not return output in meters and meters per second. As the documentation describes, it returns output in AU and AU per day.     Column Name Definition     ... ...   x x-component of position vector (float, au, X)   vx x-component of velocity vector (float, au/d, VX)   ... ...     Source. Later on, you're using a gravitational constant which is only appropriate for units of kilogram, meter, and second. G = 6.67430e-11  The simplest way to fix this is to convert the units that Horizons uses into meters and meters per second. Alternatively, you could change G to a value suitable for AU and days. 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 }  You should get the following result.