Angle problem between solar system orbital using keplerian orbital elements in python

Question

I'm currently working on a project in python that aim to plot the trajectory of an asteroid in the solar system and see how a specific impact can change it's trajectory.

After searching the subject, I came accros the concept of orbital elements. So I created a code that take in parameter:

1. the Semimajor axis (AU)
2. the Orbital eccentricity
3. the Orbital inclination (deg)
4. the Longitude of ascending node (deg)
5. Longitude of perihelion (deg)

This information is stored in a dictionnary such as:

dic_planetary_datasheet = { 'Mercure': { 'a': 0.38709893, 'e': 0.20563069, 'i': 7.00487, 'w': 77.45645, 'omega': 48.33167, 'color': 'indianred' }, 'Venus': { 'a': 0.72333199, 'e': 0.00677323, 'i': 3.39471, 'w': 131.53298, 'omega': 76.68069, 'color': 'wheat' }, 'Terre': { 'a': 1.00000011, 'e': 0.01671022, 'i': 0.00005, 'w': 102.94719, 'omega': -11.26064, 'color': 'cornflowerblue' }, 'Mars': { 'a': 1.52366231, 'e': 0.09341233, 'i': 1.85061, 'w': 336.04084, 'omega': 49.57854, 'color': "sandybrown" } }``` 

I already checked on different website for the value and they seem correct. I then convert this information to cartesian-heliocentric coordinates and I plot the result:

 import numpy as np import matplotlib.pyplot as plt from mpl_toolkits import mplot3d AU = 149597870.7 # astronomical unit [km] G = 6.67428e-11 # Newtonian constant of gravitation [m^3*kg^-1*s^-2] M_soleil = 1.9884e30 # solar mass [kg] mu = G * M_soleil # standard gravitational parameter [m^3*s^-2] def from_DEG_to_RAD(deg_value): return ((deg_value * np.pi)/180) def convert_KOP_to_cartesian(a, e, i, omega, w): ''' Convert from keplerian orbital object param. to cartesian state vector for 1 revolution :param a: Semi-major axis [AU] -> [km] :param e: Eccentricity [1] :param i: Inclination [deg] -> [rad] :param omega: Longitude of ascending node (LAN) [deg] -> [rad] :param w: Argument of perihelion [deg] -> [rad] :return: x, y, z (Cartesian coordinates in km) ''' a_km = a*AU num_points = 10000 x = [] y = [] z = [] x_p = np.zeros(num_points) y_p = np.zeros(num_points) z_p = np.zeros(num_points) i_rad = from_DEG_to_RAD(i) w_rad = from_DEG_to_RAD(w) omega_rad = from_DEG_to_RAD(omega) n = np.sqrt(mu/a_km**3) for M in np.linspace(0, 2 * np.pi, num_points): E = M for _ in range(100): # Méthode itérative de Newton-Raphson E = E - (E - e * np.sin(E) - M) / (1.0 - e * np.cos(E)) nu = 2*np.arctan(np.sqrt((1.0 + e)/(1.0 - e))*np.tan(E/2.0)) h = np.sqrt(mu * a_km * (1 - e ** 2)) r = a_km * (1 - e * np.cos(E)) x.append( r*( np.cos(w_rad)*np.cos(nu+omega_rad) - np.sin(w_rad)*np.sin(nu+omega_rad)*np.cos(i_rad) ) ) y.append( r*( np.sin(w_rad)*np.cos(nu+omega_rad) + np.cos(w_rad)*np.sin(nu+omega_rad)*np.cos(i_rad) ) ) z.append( r*( np.sin(nu+omega_rad)*np.sin(i_rad) ) ) return x, y, z def plot_orbit_KOE_other_YR4(dic_planetary_datasheet): for planet, KOE in dic_planetary_datasheet.items(): x, y, z = convert_KOP_to_cartesian(KOE['a'], KOE['e'], KOE['i'], KOE['omega'], KOE['w']) ax.plot3D(x, y, z, label=f'Trajectoire {planet}', color=KOE['color']) #asteroid trajectory x_a, y_a, z_a = convert_KOP_to_cartesian(2.516, 0.6615, 3.408, 134.36, 271.37) fig = plt.figure() ax = plt.axes(projection='3d') ax.scatter(0,0,0, label='Soleil', color='yellow', marker='o') ax.plot3D(x_a, y_a, z_a, label=f'2024 YR4', color = 'black') plot_orbit_KOE_other_YR4(dic_planetary_datasheet) ax.set_title('3D trajectories (KOE method)') plt.show() 

The problem is I get very big angle between all the trajectory and after searching I can't find a solution to that. I compared it with JPL simulation and they only have a few degrees in differences... My result using this code: JPL simulation result

If someone can help me with that it would be lovely !

Ps: If anyone know, can we get the same result(with angles and elliptic orbit) using newton method to solve a system of equations based on Netwon Gravitational Law because I only get spherical and plane trajectory

Answer 1

...aim to plot the trajectory of an asteroid in the solar system and see how a specific impact can change it's trajectory.

This is a great project, and an excellent way to learn a lot along the way!

There are several ways to treat the planets. The primary effect on the asteroid's trajectory is the Sun. Unless the asteroid passes close to a planet, their effect on the asteroid is real and important, but usually small (unless you are including Earth impact, where Earth's effect is of course now critical). So you can either

1. EASIEST: Simply propagate the asteroid's trajectory around the Sun twice; once with an impactor, and once without. That will give you a pretty good idea of the effect of the impact by looking at the difference in those two trajectories. You can then include in your plots the trajectories of the planets from an ephemeris from the Python package Skyfield or from JPL's Horizons which has a Python API as well. I am lazy so I just download a table as text. Use heliocentric coordinates for simplicity. That means you have a one-body orbit (asteroid around the Sun at the origin) and from Horizons and Skyfield you ask for the heliocentric coordinates for the planet (not barycentric).
2. BETTER: Go back and include the effects gravity from the planets in your orbital calculation. Instead of just a simple central force from the Sun, you include forces from all, or most of the planets. You still use interpolated ephemeris positions for all the planets to calculate their forces on your asteroid. In other words, you don't need to consider the asteroid's effects on the orbits of the planets.
3. BEST, HARDEST, MOST CHALLENGING: Forget the ephemerides, and Do a complete n-body simulation yourself. Consider this a "holy grail" project because it's a little challenging. I did it here How to calculate the planets and moons beyond Newtons's gravitational force?) once. just to see if I could do it, and now that I did I'll never do it again because those JPL ephemerides that Horizons and Skyfield are based on are so good and so useful, it would be silly to compete.

Forget Kelperian orbital elements

There is only one really good, straightforward, and easy to implement a solution to this problem. Forget Kelperian orbital elements. Get some initial conditions. i.e. state vectors $(x, y, z, v_x, v_y, v_z)$ for your asteroid at some starting time. Then integrate the motion using some initial value problem solver like SciPy's solve_ivp

Write the equations of motion as two first-order equations

$$\mathbf{\dot{x}} = \mathbf{v}$$

$$\mathbf{\dot{v}} = \mathbf{a} = -GM_{\odot} \frac{\mathbf{\hat{r}}}{r^2}= -GM_{\odot} \frac{\mathbf{x}}{|x|^3}$$

and solve them with Python as shown below.

Note that we should always use published values for the standard gravitational parameter because for astronomical bodies, we can measure the product $GM$ much more accurately than the masses and $G$ separately.

Python script with some helpful links in comments:

import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy.integrate import solve_ivp def deriv(time, state, GM): x, v = state.reshape(2, -1) acc = -GM * x / (x**2).sum()**1.5 # 1/r^2, vector pointing towards the sun return np.hstack([v, acc]) # https://en.wikipedia.org/wiki/Standard_gravitational_parameter GM_Sun = 1.327124400189E+20 # m^3/s^2, # assume Sun remains at origin # initial state vector for asteroid x0 = np.array([200E+09, 0, 0]) # about 1.3 AU v0 = np.array([0, 27E+03, 10E+03]) # some reasonable speed and direction state_0 = np.hstack([x0, v0]) # t_eval are the times for which you are requesting interpolated output # the internal, variable timesteps are controlled by the integrator. t_eval = np.linspace(0, 7E+07, 201) # seconds t_span = t_eval[[0, -1]] # first and last time method = 'DOP853' # for non-stiff (no super-close approaches) https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html args = (GM_Sun, ) rtol = 1E-08 result = solve_ivp(fun=deriv, t_span=t_span, y0=state_0, method=method, t_eval=t_eval, args=args, rtol=rtol) xyz, vxyz = result.y.reshape(2, 3, -1) rmax = np.linalg.norm(xyz, axis=0).max() # furthest distance from Sun AU = 149597870700. # meters https://en.wikipedia.org/wiki/Astronomical_unit if True: fig = plt.figure() ax = fig.add_subplot(1, 1, 1, projection='3d', proj_type = 'ortho') origin = np.zeros(3) # the "*" unpacks ax.plot(*origin, 'o', ms=20, mfc='yellow', mec='red') ax.plot(*xyz[:, :1] / AU, 'ok') ax.plot(*xyz / AU, '-k') # https://stackoverflow.com/a/68242226/3904031 # make a nice cube hw = 1.1 * rmax / AU ax.set_xlim(-hw, hw) ax.set_ylim(-hw, hw) ax.set_zlim(-hw, hw) ax.set_xlabel('x (AU)') ax.set_ylabel('y (AU)') ax.set_zlabel('z (AU)') ax.set_aspect('equal') plt.show()