Encontrar coordenadas de periapsis centradas en la Tierra

Actualmente estoy escribiendo un programa que acepta un vector de estado, calcula la elipse Kepleriana correspondiente y luego dibuja la elipse.

Pero estoy atascado en el siguiente punto. Actualmente necesito algún valor para decirle al programa cuánto girar la elipse orbital. Para esto necesito un valor del ángulo desde el eje x (coordenadas ECI) hasta el periapsis.

El problema es que el argumento del periapsis no me sirve si la órbita no está inclinada.

¿Hay alguna forma de calcular este ángulo o calcular las coordenadas del periapsis para poder lograr esto?

En este momento estoy trazando la órbita usando la función matplotlib add ellipse (python) que requiere algún ángulo en el que se debe rotar la elipse, rotation_deges el ángulo que estoy buscando.

orbit = patches.Ellipse((Cx,Cy), major, minor, rotation_deg, 
                        facecolor="none", edgecolor="k", linestyle="--")

En el siguiente ejemplo, elegí cuidadosamente un vector de estado inicial con$x=\text{periapsis}, y=0, z=0$tal que el ángulo es cero, pero necesito generalizar para cualquier vector de estado arbitrario dentro de una órbita dada.

r = np.array([r_earth+500000,      0,   0])
v = np.array([             0,   9000,   0])

Aquí está mi código hasta ahora:

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