Astrodynamics
Calculating the state vector from the orbital elements and an instant of time.
Posts  1 - 1  of  1
Jenab6
We'll use meters for distances (a, r), days for intervals of time (P), Julian date format for moments of time (t,T) and radians for angles.

The Keplerian orbital elements.
[ a, e, i, Ω, ω, T ] @ time = t

The solar gravitational parameter, GM.
GM = 1.32712440018e20 m³ sec⁻²

The period of the orbit, P, days.
P = (π / 43200) sqrt[ a³ / (GM) ]

The mean anomaly, m, of the object in the orbit at time = t.
m₀ = (t − T) / P
m = 2π [m₀ − integer(m₀)]

The eccentric anomaly, u, of the object in the orbit at time = t. (Danby's method is shown.)

u = m
Repeat...
. U₀ = u
. F₀ = U₀ − e sin U₀ − m
. F₁ = 1 − e cos U₀
. F₂ = e sin U₀
. F₃ = e cos U₀
. D₁ = −F₀ / F₁
. D₂ = −F₀ / [ F₁ + D₁ F₂ / 2 ]
. D₃ = −F₀ / [ F₁ + D₁ F₂ / 2 + (D₂)² F₃ / 6 ]
. u = U₀ + D₃
UNTIL |u−U₀| 0 then angle = +π/2
. if x = 0 and y = 0 then angle = 0
. if x = 0 and y 0 and y > 0 then angle = atn(y/x)
. if x > 0 and y < 0 then angle = atn(y/x) + 2π
Arctan = angle

Rotate the triple-prime position vector about the z''' axis by the argument of the perihelion, ω.
x'' = x''' cos ω − y''' sin ω
y'' = x''' sin ω + y''' cos ω
z'' = z''' = 0

Rotate the double-prime position vector about the x'' axis by the inclination, i.
x' = x''
y' = y'' cos i
z' = y'' sin i

To find the heliocentric position of the object in the orbit at time t, rotate the single-prime position vector about the z' axis by the longitude of the ascending node, Ω.
x = x' cos Ω − y' sin Ω
y = x' sin Ω + y' cos Ω
z = z'

k = sqrt { GM / [ a (1 − e²) ] }
k is a speed in meters per second.

The canonical velocity vector of the object in its orbit at time t.
Vx''' = −k sin θ
Vy''' = k (e + cos θ)
Vz''' = 0

Rotate the triple-prime velocity vector about the z''' axis by the argument of the perihelion, ω.
Vx'' = Vx''' cos ω − Vy''' sin ω
Vy'' = Vx''' sin ω + Vy''' cos ω
Vz'' = Vz''' = 0

Rotate the double-prime velocity vector about the x'' axis by the inclination, i.
Vx' = Vx''
Vy' = Vy'' cos i
Vz' = Vy'' sin i

To find the sun-relative velocity in ecliptic coordinates, at time t, rotate the single-prime velocity vector about the z' axis by the longitude of the ascending node, Ω.
Vx = Vx' cos Ω − Vy' sin Ω
Vy = Vx' sin Ω + Vy' cos Ω
Vz = Vz'

The state vector of the object in the orbit at time t.
[ x, y, z, Vx, Vy, Vz ]

Jerry Abbott
Save
Cancel
Reply
 
x
OK