Finden der erdzentrierten Koordinaten der Periapsis

Ich schreibe gerade ein Programm, das einen Zustandsvektor akzeptiert, die entsprechende Keplersche Ellipse berechnet und dann die Ellipse zeichnet.

Aber ich hänge an folgendem Punkt fest. Ich brauche derzeit einen Wert, um dem Programm mitzuteilen, um wie viel die Orbitallipse gedreht werden soll. Dazu benötige ich einen Wert des Winkels von der x-Achse (ECI-Koordinaten) zur Periapsis.

Das Problem ist, dass das Argument der Periapsis mir nicht hilft, wenn die Umlaufbahn nicht geneigt ist.

Gibt es eine Möglichkeit, diesen Winkel zu berechnen oder die Koordinaten der Periapsis zu berechnen, damit ich dies erreichen kann?

Im Moment zeichne ich die Umlaufbahn mit der Funktion matplotlib add ellipse (Python), die einen Winkel erfordert, um den die Ellipse gedreht werden sollte, rotation_degist der Winkel, nach dem ich suche.

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

Im folgenden Beispiel habe ich sorgfältig einen Anfangszustandsvektor mit ausgewählt$x=\text{periapsis}, y=0, z=0$so dass der Winkel Null ist, aber ich muss für jeden beliebigen Zustandsvektor innerhalb einer bestimmten Umlaufbahn verallgemeinern.

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

Hier mein bisheriger Code:

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