Basic 2D orbit plotting using matplotlib

basicorbit.py

# Orbit plotting example in 2D
# www.chryseus.co.uk

import numpy as np
import matplotlib.pyplot as plt

a = 1 # Semi-major axis
ε = 0.4 # Eccentricity
points = 100 # Number of points to use in plot
step = (2*np.pi)/points # Step size
# Arrays for X and Y data
xdata = np.array([])
ydata = np.array([])

for i in range(0,points+1):
    θ = i * step # Theta will increase by step with each loop until it reaches 2*pi
    r = (a * (1-ε**2)) / (1+ε*np.cos(θ)) # Polar orbit equation
    # Convert to rectangular cartesian coordinates
    x = r * np.cos(θ)
    y = r * np.sin(θ)
    # Add X and Y position to arrays
    xdata = np.append(xdata, x)
    ydata = np.append(ydata, y)

# Calculate the minimum (periapsis) and maximum (apoapsis) distance
rmin = a * (1-ε**2) / (1+ε)
rmax = a * (1-ε**2) / (1-ε)
# Calculate other parameters
b = a * np.sqrt(1-ε**2) # Semi-minor axis
c = np.sqrt(a**2 - b**2) # Linear eccentricity
l = b**2 / a # Semi-latus rectum

# Plotting
fig,ax = plt.subplots() # Prepare a new plot
ax.set_xlim(-2,2) # Axis x limit
ax.set_ylim(-2,2) # Axis y limit
ax.set_aspect(1) # Ensure a 1:1 aspect ratio to avoid visual distortion
ax.grid(True) # Turn on grid lines
ax.plot(xdata, ydata, 'b') # Plot the orbit in blue
ax.plot(rmin, 0, 'ro') # Periapsis point
ax.plot(-rmax, 0, 'ro') # Apoapsis point
ax.plot([-rmax,rmin],[0,0],'g', linestyle='solid') # Draw green major axis
ax.plot([-c,-c], [b,-b], 'r') # Draw red minor axis
ax.plot([0,0], [l,-l],color='purple') # Draw latus rectum
ax.plot(0, 0, 'yo') # Yellow focal point
plt.show()