What are orbital state vectors?
The orbital state vectors are cartesian vectors of position(commonly referred as r, it describes the position of the body in the chosen frame of reference) and velocity(commonly referred as v, it describes its velocity in the same frame at the same time). Together with their time(commonly referred as epoch, which is a moment in time used as a reference point for some time-varying astronomical quantity) they can determine the trajectory of an orbiting body in space.
What are the Keplerian elements that we need in order to convert them to cartesian state vectors?
- the semi-major axis a (in meters)
- the eccentricity e
- the argument of periapsis ω (in radians)
- the longitude of the ascending node Ω (in radians)
- the inclination i (in radians)
- the mean anomaly M0=Mt(0) (in radians) at epoch t0 (in Julian dates)
- considered epoch t (in Julian dates) different from t0
- the mu which is given by the gravitational constant G and the Mass of the central body (mu = G*M) If the central body is earth for example, we have:
M = 5.972 × 10^24, G = 6.67408 × 10^-11, mu = 3.9860044188 × 10^14
Let’s start !
Now, that we know what we need we can start. Firstly we need to make a function and pass our parameters:
local function ConvertToCartesian( e, a, i, O, w, t, t0, M0)
end
As you have probably have guessed, e is the eccentricity, a is the semi-major axis, O is the longitude of the ascending node , w is the argument of periapsis ,M0 the mean anomaly at epoch t0, t is the considered epoch and t0 the different from t epoch, as we have seen.
Now, lets continue:
local function ConvertToCartesian( e, a, i, O, w, t, t0, M0)
G = 6.67408 * 10^-11
M = 5.972 * 10^24 -- We will calculate the mu of earth
local mu = G * M
end
We said that t must be different from t0. Lets make an if loop then.
local function ConvertToCartesian( e, a, i, O, w, t, t0, M0)
G = 6.67408 * 10^-11
M = 5.972 * 10^24 -- We will calculate the mu of earth
local mu = G * M
if(t==t0) then
t = t0
Mt = 0
else
-- We will continue this
end
end
Now, there are some things we need to do. Firstly, we said that t and t0 are in Julian Days. Now, in order to calculate the Δt which is the difference t - t0 , we need to convert Julian dates to seconds by using this:
So, the Δt will be:
Knowing the Δt will help us calculate the mean anomaly applied to the argument t
Now, this might seem difficult, but in reality it is really easy. We already know the M0 because it is a parameter of our function so it is given to us, we have calculated the Δt and the mu and the a are also given to us. So lets continue our program :
local function ConvertToCartesian( e, a, i, O, w, t, t0, M0)
G = 6.67408 * 10^-11
M = 5.972 * 10^24 -- We will calculate the mu of earth
local mu = G * M
if(t==t0) then
t = t0
Mt = 0
else
local dt = 86400 * (t-t0)
local Mt = M0 + dt * math.sqrt(mu/a^3)
-- We will continue this
end
end
Now, we need to solve the Kepler’s equation 
for the eccentric anomaly
. The problem with solving this equation is that it is a transcendental equation. In order to solve this numerically we need to use an appropriate method, in our case the Newton-Raphson method, which is a way to find a good approximation for the root of a real-valued function f(x) = 0 . You can see more here.
So, lets find the root of this function:
By applying the Newton’ method we have:
local function ConvertToCartesian( e, a, i, O, w, t, t0, M0)
G = 6.67408 * 10^-11
M = 5.972 * 10^24 -- We will calculate the mu of earth
local mu = G * M
if(t==t0) then
t = t0
Mt = 0
else
local dt = 86400 * (t-t0)
local Mt = M0 + dt * math.sqrt(mu/a^3)
local maxIter = 30 -- The desired accuracy
local E = Mt
local F = E - e * math.sin(E) - Mt -- Just a variable F
local i = 0
while ((math.abs(F)>delta) and (i<maxIter)) do
E = E - F/(1 -e* math.cos(E))
F = E - e * math.sin(E) - Mt
i = i + 1
end
-- Next
end
end
Then, we will need to calculate the true anomaly nu. Note, that it can take two values:
Let’s implement this in our code:
local function ConvertToCartesian( e, a, i, O, w, t, t0, M0)
G = 6.67408 * 10^-11
M = 5.972 * 10^24 -- We will calculate the mu of earth
local mu = G * M
if(t==t0) then
t = t0
Mt = 0
else
local dt = 86400 * (t-t0)
local Mt = M0 + dt * math.sqrt(mu/a^3)
local maxIter = 30 -- The desired accuracy
local E = Mt
local F = E - e * math.sin(E) - Mt -- Just a variable F
local i = 0
while ((math.abs(F)>delta) and (i<maxIter)) do
E = E - F/(1. -e* math.cos(E))
F = E - e * math.sin(E) - m
i = i + 1
end
nu = 2 * math.atan2(math.sqrt(1+e) * math.sin(E/2), math.sqrt (1 - e) * math.cos(E / 2))
-- Next
end
end
We have calculated the mean anomaly, yet we need to calculate and the distance to the central body:
local function ConvertToCartesian( e, a, i, O, w, t, t0, M0)
G = 6.67408 * 10^-11
M = 5.972 * 10^24 -- We will calculate the mu of earth
local mu = G * M
if(t==t0) then
t = t0
Mt = 0
else
local dt = 86400 * (t-t0)
local Mt = M0 + dt * math.sqrt(mu/a^3)
local maxIter = 30 -- The desired accuracy
local E = Mt
local F = E - e * math.sin(E) - Mt -- Just a variable F
local i = 0
while ((math.abs(F)>delta) and (i<maxIter)) do
E = E - F/(1. -e* math.cos(E))
F = E - e * math.sin(E) - m
i = i + 1
end
nu = 2 * math.atan2(math.sqrt(1+e) * math.sin(E/2), math.sqrt (1 - e) * math.cos(E / 2))
rc = a * (1 - e * math.cos(E))
-- next
end
end
Finally, we can start calculating the position and velocity vectots. First, remember that in three dimensions there are 3 axes, the x,y and z ones. In order to obtain the position vector o and the velocity vector odot we need the following:
Which simple means that :
- o.x = rc * cos(nu)
- o.y = rc * sin(nu)
- o.z = 0 (because rc*0=0)
- odot.x = (sqrt (mu * a) / rc) * sin(E)
- odot.y = (sqrt (mu * a) / rc) * sqrt(1-e^2)*cos(E)
- odot.z = 0
Let’s do this in our code:
local function ConvertToCartesian( e, a, i, O, w, t, t0, M0)
G = 6.67408 * 10^-11
M = 5.972 * 10^24 -- We will calculate the mu of earth
local mu = G * M
if(t==t0) then
t = t0
Mt = 0
else
local dt = 86400 * (t-t0)
local Mt = M0 + dt * math.sqrt(mu/a^3)
local maxIter = 30 -- The desired accuracy
local E = Mt
local F = E - e * math.sin(E) - Mt -- Just a variable F
local i = 0
while ((math.abs(F)>delta) and (i<maxIter)) do
E = E - F/(1. -e* math.cos(E))
F = E - e * math.sin(E) - m
i = i + 1
end
nu = 2 * math.atan2(math.sqrt(1+e) * math.sin(E/2), math.sqrt (1 - e) * math.cos(E / 2))
rc = a * (1 - e * math.cos(E))
o = Vector3.new(rc * math.os(nu), rc * math.sin(nu), 0)
odot = Vector3.new(math.sin (E), math.sqrt(1 - e * e) * math.cos(E), 0)
odot = (math.sqrt (mu * a) / rc) * odot
-- next
end
end
The only thing that is left to do, is to transform the o and the odot in bodycentric rectangular coordinates r and rdot. This can be done using the following( which is basically a transformation sequence of the rotation matrices Rz(φ), Rx(φ))
This just means that:
-
rx = ( odot.x * (cos (w) * cos (O) - sin (w) * cos (i) * sin (O)) -
odot.y * (sin (w) * cos (O) + cos (w) * cos (i) * sin (O))) -
ry = (odot.x * (cos (w) * sin (O) + sin (w) * cos (i) * cos (O)) +
odot.y * (cos (w) * cos (i) * cos (O) - sin (w) * sin (O))) -
rz = (odot.x * (sin (w) * sin (i)) + o.y * (cos (w) * sin (i)))
-
rdotx = ( odot.x * (cos (w) * cos (O) - sin (w) * cos (i) * sin (O)) -
odot.y * (sin (w) * cos (O) + cos (w) * cos (i) * sin (O))) -
rdoty = (odot.x * (cos (w) * sin (O) + sin (w) * cos (i) * cos (O)) +
odot.y * (cos (w) * cos (i) * cos (O) - sin (w) * sin (O))) -
rdotz = (odot.x * (sin (w) * sin (i)) + o.y * (cos (w) * sin (i)))
Adding this in our code, we have finally :
local function ConvertToCartesian( e, a, i, O, w, t, t0, M0)
G = 6.67408 * 10^-11
M = 5.972 * 10^24 -- We will calculate the mu of earth
local mu = G * M
if(t==t0) then
t = t0
Mt = 0
else
local dt = 86400 * (t-t0)
local Mt = M0 + dt * math.sqrt(mu/a^3)
local maxIter = 30 -- The desired accuracy
local E = Mt
local F = E - e * math.sin(E) - Mt -- Just a variable F
local i = 0
while ((math.abs(F)>delta) and (i<maxIter)) do
E = E - F/(1. -e* math.cos(E))
F = E - e * math.sin(E) - m
i = i + 1
end
nu = 2 * math.atan2(math.sqrt(1+e) * math.sin(E/2), math.sqrt (1 - e) * math.cos(E / 2))
rc = a * (1 - e * math.cos(E))
o = Vector3.new(rc * math.os(nu), rc * math.sin(nu), 0)
odot = Vector3.new(math.sin (E), math.sqrt(1 - e * e) * math.cos(E), 0)
odot = (math.sqrt (mu * a) / rc) * odot
rx = ( o.x * (math.cos(w) * math.cos(O) - math.sin (w) * math.cos(i) * math.sin(O)) -
o.y * (math.sin(w) * math.cos(O) + math.cos(w) * math.cos(i) * math.sin(O)))
ry = (o.x * (math.cos(w) * math.sin(O) + math.sin(w) * math.cos(i) * math.cos(O)) +
o.y * (math.cos(w) * math.cos(i) * math.cos(O) - math.sin(w) * math.sin(O)))
rz = (o.x * (math.sin(w) * math.sin(i)) + o.y * (math.cos(w) * math.sin(i)))
rdotx = ( odot.x * (math.cos(w) * math.cos(O) - math.sin(w) * math.cos(i) * math.sin(O)) -
odot.y * (math.sin(w) * math.cos(O) + math.cos(w) * math.cos(i) * math.sin(O)))
rdoty = (odot.x * (math.cos(w) * math.sin(O) + math.sin(w) * math.cos(i) * math.cos(O)) +
odot.y * (math.cos(w) * math.cos(i) * math.cos(O) - math.sin(w) * math.sin(O)))
rdotz = (odot.x * (math.sin(w) * math.sin(i)) + odot.y * (math.cos(w) * math.sin(i)))
r = Vector3.new(rx,ry,rz)
rdot = Vector3.new(rdotx,rdoty,rdotz)
return r , rdot
end
end