1
$\begingroup$

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

example plot of elliptical orbit

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()
$\endgroup$
  • $\begingroup$ If you have proper state vectors, you just plot them and that's your orbit. Do you really mean you have a table of $x, y, z, v_x, v_y, v_z$ values? Maybe you could show an example of what you have to plot so it's easier to understand why you still need more information. $\endgroup$ – uhoh Apr 1 '19 at 23:25
  • $\begingroup$ The way I am going about it is converting the vectors into orbital elements and plotting an ellipse from those values. I am using the matplotlib library for python to do this plotting. I am not sure how I would simply plot the vector and see an orbit. $\endgroup$ – Griffin J Apr 2 '19 at 15:51
  • $\begingroup$ I don't think we can discuss further unless you show a sample of what you have. A plot of an orbit is really a series of line segments connecting a sequence of $x, y, z$ points. See this answer for example. That orbit is a sequence of 201 state vectors, with the $x, y, z$ part of the state vectors plotted. $\endgroup$ – uhoh Apr 2 '19 at 22:30
  • $\begingroup$ I updated my answer, if you want you can take a look at my code to see what you think. $\endgroup$ – Griffin J Apr 4 '19 at 17:45
  • $\begingroup$ Thanks very much for the update and script! I've adjusted the wording a bit and added some formatting. Now I see that each orbit is from one initial state vector and not from state vectors (plural). Can you clarify if you need this to work for any arbitrary initial state vector, including all RAAN and inclinations? Right now there's an answer which seems helpful but I'm not sure it covers your case. It would be great if you could leave a comment for that author directly under their answer. Thanks! $\endgroup$ – uhoh Apr 5 '19 at 1:09
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ This is the start of a great answer, but the answer is in the links, not here, making this a link-only answer mostly, which is generally discouraged in Stack Exchange. If/when the links rot, there is no answer here for future readers to see. Can you either grab a block-quite of a small passage or an equation or two and summarize the answer(s) here in you post? Thanks! $\endgroup$ – uhoh Apr 1 '19 at 22:57
  • $\begingroup$ You've certainly written quite thorough answers before! $\endgroup$ – uhoh Apr 1 '19 at 23:13
  • $\begingroup$ @uhoh : I understand a 'link-only answer' to be one that gives a URL without identifying what should be addressed by it. Here I have given citations to the publications by title and author &c., so that even if the links in themselves 'rot', the cited references can still be identified and searched for. So I believe this is not in fact a link-only answer. The difficulty about selecting formulae for the questioner is that so little has been said about the problem that it's not clear which computational operations would be most relevant. $\endgroup$ – terry-s Apr 1 '19 at 23:14
  • $\begingroup$ If the question is too general to answer, then add a comment to the question asking for clarification. I'm pretty sure the term link-only includes citation-only as well, the point being that the answer is not here in the answer but is somewhere else outside of Stack Exchange. I've modified your preface, then added a comment to the question asking for clarification. $\endgroup$ – uhoh Apr 1 '19 at 23:22
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I've added some formatting here. Have a look at the question, it's been updated with a full script and with some clarifications in wording. There are additional comments there as well. I wonder if you can expand on this to include inclination also? Thanks! $\endgroup$ – uhoh Apr 5 '19 at 1:11
  • $\begingroup$ Thanks for the clarification. Doesn't using the argument of periapsis as the rotation value assume that the line of nodes is along the x-axis? I'm just struggling to understand why the argument of periapsis works as the rotation angle. $\endgroup$ – Griffin J Apr 5 '19 at 16:56
  • $\begingroup$ Argument of Periapsis is the angle, in the plane of the orbit, measured from the Right Ascension, of the Ascending Node, through the center of the body being orbited to the Periapsis, in the direction of travel around the orbit. So if your orbit is equatorial and prograde (i =0), then raan=0, and thus the Argument of Periapsis is the angle you want. $\endgroup$ – notovny Apr 5 '19 at 17:07
  • $\begingroup$ Awesome! Thank you for the explanation. And this even works with the ECI coordinates that I am using for the graph? I don't have to be graphing the orbit from a view perpendicular to the plane of the orbit? $\endgroup$ – Griffin J Apr 5 '19 at 17:41
  • $\begingroup$ @GriffinJ unfortunately, using argument of periapsis alone only works for the non-inclined, orbit viewed from the perpendicular case, which is what I thought you were asking. Apologies. I'm not personally familiar with matplotlib, so I can't really make a good recommendation on how to deal with making it draw an elliptical orbit from an arbitrary viewpoint. $\endgroup$ – notovny Apr 5 '19 at 23:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.