2D Static Orbits

You will need several libraries to help plot these orbits:

import numpy as np
from PyAstronomy import pyasl
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import mpl_toolkits.mplot3d.axes3d as p3

The beauty of using a third-party library (that has been properly vetted) is that the library can handle the hard work, and you only need one line to create the elliptical orbit object KeplerEllipse. you will start with a basic ellipse:

• semi-major axis of 1.0 units
• a time period of 1.0 units
• an eccentricity of 0.50
• an inclination of 10 degrees
• everything else 0.

You will also use Numpy to create time intervals to be used in determining position later.

orbit = pyasl.KeplerEllipse(a=1.0, per=1.0, e=0.5, Omega=0.0, i=10.0, w=0.0)

t = np.linspace(0, 4, 200)

The linspace() function accepts a starting point, ending point, and a number of samples that you want in between. In this case, it calculates 200 samples between 0 and 4 (inclusive). It means that you do not have to calculate the actual step size; Numpy will do that for you.

All that is left to do is call the xyzpos() method, which accepts a time array argument and returns a 3 column array (basically a list) of the same length as the input time array. It is important to note that the returned array is in the format Y coordinates, X coordinates, Z coordinates. For a 2D plot, you will ignore the Z coordinates, so our viewpoint is looking down at the North Pole of Earth from above Earth.

pos = orbit.xyzPos(t)

Congratulations, that’s all of the orbital mechanics that is required! The rest of the program will deal with Matplotlib and setting the proper parameters. Because of the convention used within PyAstronomy, the first point in pos corresponds to perigee, the closest point to the Earth. You will plot perigee separately, so you have your own reference point when you tweak the orbital parameters. The array-like notation contained within pos means that you will have to use a peculiar double-colon ∷ to get all of the elements within that array for that coordinate system. Earth will be at (0, 0, 0) in this coordinate system, and you will plot it with a big marker size so that you can see it on the plot.

plt.plot(0, 0, 'bo', markersize=9, label="Earth")

plt.plot(pos[::, 1], pos[::, 0], 'k-', label="Satellite Trajectory")

plt.plot(pos[0, 1], pos[0, 0], 'r*', label="Periapsis")

Since the pos object’s array coordinates are in the system (Y, X, Z), to get the X coordinates, you need to call the “first” (in reality, the second) element in the array. Add a legend and a title, and let’s show the plot:

plt.legend(loc="upper right")
plt.title('Orbital Simulation')
plt.show()

Now, you can see the plot: