I am currently writing a program that accepts a state vector, calculates the corresponding Keplerian ellipse, and then draws the ellipse.
But I'm stuck at the following point. I currently need some value to tell the program how much to rotate the orbital ellipse. For this I need a value of the angle from the x-axis (ECI coordinates) to the periapsis.
The problem is that argument of periapsis does not help me if the orbit is not inclined.
Is there some way of calculating this angle or calculating the coordinates of the periapsis so I can achieve this?
Right now I am plotting the orbit by using the matplotlib add ellipse function (python) which requires some angle that the ellipse should be rotated, rotation_deg
is the angle I am looking for.
orbit = patches.Ellipse((Cx,Cy), major, minor, rotation_deg,
facecolor="none", edgecolor="k", linestyle="--")
In the example below, I've carefully chosen an initial state vector with $x=\text{periapsis}, y=0, z=0$ such that the angle is zero, but I need to generalize for any arbitrary state vector within a given orbit.
r = np.array([r_earth+500000, 0, 0])
v = np.array([ 0, 9000, 0])
Here is my code so far:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
G = 6.67e-11
M_earth = 5.972e24
mu_earth = G*M_earth
r_earth = 6.3781e6
def plot_setup_earth():
global f, ax
f, ax = plt.subplots()
circle = plt.Circle((0,0),r_earth,color = "b")
ax.set_aspect("equal", "box")
ax.add_artist(circle)
ax.set_title("Vehicle Orbit (ECI Coordinates)")
ax.set_xlabel("X (meters)")
ax.set_ylabel("Y (meters)")
def elements_from_vectors(r,v,radius,mu):
global apoapsis, periapsis, a, b, e
K = np.array([0,0,1])
h_vec = np.cross(r,v) #momentum vector
n_vec = np.cross(K,h_vec) #line of nodes vector
e_vec = (1/mu)* np.cross(v,h_vec) - (r/np.linalg.norm(r)) #eccentricity vector
E = ((np.linalg.norm(v)**2)/2) - (mu/np.linalg.norm(r)) #orbital energy
e = np.linalg.norm(e_vec) #eccentricity
a = (-1)*(mu/(2*E)) #semimajor axis
b = a*np.sqrt(1-(e**2))
p = ((np.linalg.norm(h_vec)**2)/mu) #????
i = np.arccos(np.dot((h_vec/np.linalg.norm(h_vec)),K)) #Inclination (rad)
i_deg = np.degrees(i) #Inclination (deg)
#Right ascension of the ascending node
if i == 0:
raan = 0
elif n_vec[1] >= 0:
raan = np.arccos(n_vec[0]/np.linalg.norm(n_vec))
elif n_vec[1] <0:
raan = (2*np.pi) - np.arccos(n_vec[0]/np.linalg.norm(n_vec))
#Argument of periapsis
if i == 0:
w = 0
elif e_vec[2] >= 0:
w = np.arccos(np.dot(n_vec,e_vec)/(np.linalg.norm(n_vec)*np.linalg.norm(e_vec)))
elif e_vec[2] < 0:
w = (2*np.pi)-np.arccos(np.dot(n_vec,e_vec)/(np.linalg.norm(n_vec)*np.linalg.norm(e_vec)))
#True anomoly
if np.dot(r,v) >= 0:
f = np.arccos(np.dot(r,e_vec)/(np.linalg.norm(r)*np.linalg.norm(e_vec)))
if np.dot(r,v) < 0:
f = (2*pi)-np.arccos(np.dot(r,e_vec)/(np.linalg.norm(r)*np.linalg.norm(e_vec)))
periapsis = a*(1-e)
apoapsis = a*(1+e)
def plot_orbit(apoapsis, periapsis, semimajor, semiminor, eccentricity):
major = 2*semimajor #Aka major axis
minor = 2*semiminor #Aka minor axis
ae = semimajor*eccentricity #Focus to center distance (where focus 1 is center of body being orbited)
rotation_deg = 0
rotation_rad = np.radians(rotation_deg)
Px = periapsis*np.cos(rotation_rad)
Py = periapsis*np.sin(rotation_rad)
Ax = (-1)*apoapsis*np.cos((-1)*rotation_rad)
Ay = apoapsis*np.sin((-1)*rotation_rad)
Cx = 0 - ae*np.cos(-1*rotation_rad)
Cy = 0 + ae*np.sin(-1*rotation_rad)
orbit = patches.Ellipse((Cx,Cy), major, minor, rotation_deg, facecolor = "none", edgecolor = "k", linestyle = "--")
ax.add_artist(orbit)
plt.scatter(Px,Py, color = "r",marker = "+")
plt.annotate("Periapsis %f km" %round((periapsis-r_earth)/1000,1), (Px,Py), fontsize = 9)
plt.scatter(Ax,Ay,color = "b", marker = "+")
plt.annotate("Apoapsis %f km" %round((apoapsis-r_earth)/1000,1),(Ax,Ay), fontsize = 9)
r = np.array([r_earth+500000,0,0])
v = np.array([0,9000,0])
plot_setup_earth()
elements_from_vectors(r,v,r_earth,mu_earth)
plot_orbit(apoapsis, periapsis, a, b, e)
plt.xlim(-1.5*apoapsis,4*periapsis)
plt.ylim(-1*(b*2),1*((b*2)))
plt.show()