Finding Earth-Centered Coordinates of Periapsis

Question

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()

Answer 1

Currently the question as written is too general to answer specifically. Nonetheless I'll provide some helpful recommended resources as an extended comment.

You could start with this explanation of the longitude of periapsis for an example of a set of orbital elements that does not suffer from misbehavior when the inclination approaches zero.

And/or you could use formulae for relations between state vectors and orbital elements given in Bate, Muller and White (1971) Fundamentals of Astrodynamics, ...

or for example in A E Roy, Orbital Motion (or google for other sources of the same title).

You might also find some of the formulae useful in Preliminary orbit-determination method having no co-planar singularity by R M L Baker and N H Jacoby, Celestial Mechanics 15 (1977) 137-160.

I hope that helps.

Answer 2

In the situation where the orbit is not inclined, the convention generally is that the Argument of Periapsis is the angle from the Reference Direction (in your case, the x-axis) to the Periapsis.

The code you have assumes that the Argument of Periapsis is always zero if the inclination is zero, which isn't true. Here's how I'd rewrite your Argument of Periapsis section.

    #Argument of periapsis

    if i != 0 and e!= 0: # Orbit is inclined and eccentric
        w = np.arccos(np.dot(n_vec,e_vec)/(np.linalg.norm(n_vec)*np.linalg.norm(e_vec)))
        if e_vec[2] < 0:
            w = (2*np.pi) - w
    elif e != 0: # Orbit is not inclined, but is eccentric
        w = numpy.arctan2(e_vec[1], e_vec[0])
    else: #Orbit is circular.
        w = 0.0

As such, as mentioned above, you can then use your argument of periapsis value as rotation_deg=w to rotate your ellipse in the fashion you desire, if your orbit is equatorial and prograde (inclination = 0 ) and the viewpoint is from above the Z-axis.