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() 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$)$(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
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy.integrate import solve_ivp # https://en.wikipedia.org/wiki/Standard_gravitational_parameterGM_Sun = 1.327124400189E+20 # m^3/s^2, # assume Sun remains at origin def deriv(time, state, GM): #print('time:',time) # print('state: ', state) # print('GM: ', 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, 25E+0327E+03, 10E+03]) # some reasonable speed and direction state_0 = np.hstack([x0, v0]) # t_eval areare the times for which you are requesting interpolated output # the internal, variable timesteps are controlled by the integrator. t_eval = np.linspace(0, 5E+077E+07, 201) 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) # hwrmax = (np.abs(xyz)**2)linalg.sumnorm(xyz, axis=0).max()**0.5 # uglyfurthestwaydistance from Sunhw AU = np.linalg149597870700.norm(xyz,axis=0)# meters https://en.max()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 (mAU)') ax.set_ylabel('y (mAU)') ax.set_zlabel('z (mAU)') ax.set_aspect('equal') plt.show() 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
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy.integrate import solve_ivp # https://en.wikipedia.org/wiki/Standard_gravitational_parameterGM_Sun = 1.327124400189E+20 # m^3/s^2, # assume Sun remains at origin def deriv(time, state, GM): #print('time:',time) # print('state: ', state) # print('GM: ', 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]) x0 = np.array([200E+09, 0, 0]) # about 1.3 AU v0 = np.array([0, 25E+03, 10E+03]) # some reasonable speed and direction state_0 = np.hstack([x0, v0]) # t_eval areare the times for which you are requesting interpolated output # the internal, variable timesteps are controlled by the integrator. t_eval = np.linspace(0, 5E+07, 201) 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) # hw = (np.abs(xyz)**2).sum(axis=0).max()**0.5 # uglywayhw = np.linalg.norm(xyz,axis=0).max() 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], 'ok') ax.plot(*xyz, '-k') # https://stackoverflow.com/a/68242226/3904031 # make a nice cube ax.set_xlim(-hw, hw) ax.set_ylim(-hw, hw) ax.set_zlim(-hw, hw) ax.set_xlabel('x (m)') ax.set_ylabel('y (m)') ax.set_zlabel('z (m)') ax.set_aspect('equal') plt.show() 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
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) 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() # furthestdistance 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() ...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
- 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).
- 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.
- 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 # https://en.wikipedia.org/wiki/Standard_gravitational_parameter GM_Sun = 1.327124400189E+20 # m^3/s^2, # assume Sun remains at origin def deriv(time, state, GM): # print('time: ', time) # print('state: ', state) # print('GM: ', 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]) x0 = np.array([200E+09, 0, 0]) # about 1.3 AU v0 = np.array([0, 25E+03, 10E+03]) # some reasonable speed and direction state_0 = np.hstack([x0, v0]) # t_eval are are the times for which you are requesting interpolated output # the internal, variable timesteps are controlled by the integrator. t_eval = np.linspace(0, 5E+07, 201) 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) # hw = (np.abs(xyz)**2).sum(axis=0).max()**0.5 # ugly way hw = np.linalg.norm(xyz, axis=0).max() 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], 'ok') ax.plot(*xyz, '-k') # https://stackoverflow.com/a/68242226/3904031 # make a nice cube ax.set_xlim(-hw, hw) ax.set_ylim(-hw, hw) ax.set_zlim(-hw, hw) ax.set_xlabel('x (m)') ax.set_ylabel('y (m)') ax.set_zlabel('z (m)') ax.set_aspect('equal') plt.show() 