Нахождение геоцентрических координат перицентра

В настоящее время я пишу программу, которая принимает вектор состояния, вычисляет соответствующий кеплеровский эллипс, а затем рисует эллипс.

Но я застрял в следующем моменте. В настоящее время мне нужно какое-то значение, чтобы сообщить программе, насколько нужно повернуть орбитальный эллипс. Для этого мне нужно значение угла от оси x (координаты ECI) до перицентра.

Проблема в том, что аргумент перицентра мне не помогает, если орбита не наклонена.

Есть ли способ вычислить этот угол или вычислить координаты перицентра, чтобы я мог этого добиться?

Прямо сейчас я рисую орбиту, используя функцию добавления эллипса matplotlib (python), которая требует некоторого угла, на который должен быть повернут эллипс, rotation_deg- это угол, который я ищу.

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

В приведенном ниже примере я тщательно выбрал вектор начального состояния с$x=\text{periapsis}, y=0, z=0$так что угол равен нулю, но мне нужно обобщить для любого произвольного вектора состояния в пределах заданной орбиты.

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

Вот мой код:

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