How To Reduce Keplerian Orbital Elements and a time to Heliocentric Position and Velocity
Definitions
a : semimajor axis of orbit, positive for elliptical orbits, negative for hyperbolic orbits
e : eccentricity of orbit, between 0 and 1 for elliptical orbits, greater than 1 for hyperbolic orbits
i : inclination of the orbital plane to the ecliptic plane
Ω : longitude of the ascending node; the angle, subtended at the sun, measured in the ecliptic, from the vernal equinox to the place where the orbital plane intersects the ecliptic with the motion of the orbiting body having a Z component in the direction of the north ecliptic pole
ω : argument of the perihelion; the angle, subtended at the sun, measured in the plane of the orbit and in the direction of motion, from the ascending node to the perihelion of the orbit
T : the time of perihelion passage, typically presented in Julian Date format, but can be translated into a decimal calendar date
As examples, here are the elements of Earth's orbit, for epoch JD 2458792.5 (or midnight beginning 5 November 2019)
a = 0.9999951820728348 AU
e = 0.01674899215492258
i = 0.02633205404161869°
Ω = 176.9917546445248°
ω = 286.0839149800637°
T = 2458852.774528838694 = 6h 35m 19s on 4 Jan 2020 UTC
and here are the elements of 2I/Borisov orbit, for epoch JD 2458792.5
a = −0.8513198164554499
e = 3.357068272255771
i = 44.05161909545966°
Ω = 308.1483096529710°
ω = 209.1213073058442°
T = 2458826.048866978846 = 13h 10m 22s on 8 Dec 2019 UTC
FINDING THE POSITION AND VELOCITY FROM THE ORBITAL ELEMENTS AND A SPECIFIED TIME...
The solar gravitational parameter,
GM = 1.32712440018e+20 m³ sec⁻²
The conversion factor from astronomical units to meters,
U = 1.495978707e+11 meters/AU
Let the specified time, t, be 8h 52m on 11 December 2019, or
t = JD 2458828.86944
For ELLIPTICAL orbits...
Calculate the period of the orbit, P, in days
P = (π/43200) √[(aU)³/GM]
Calculate the mean anomaly in the orbit at time t
M = 2π(t−T)/P
if M<0 then M=M+2π
Obtain an initial approximation for the eccentric anomaly, u
u = (e − e³/8 + e⁵/192) sin(M) + (e²/2 − e⁴/6) sin(2M)
+ (3e³/8 − 27e⁵/128) sin(3M) + (e⁴/3) sin(4M)
Refine the eccentric anomaly to the exact value
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 + F₄D₂²/6)
u = u + D₃
Until |u−U|<0.0000000001
Find the canonical position of the object in its own orbit
x''' = a[cos(u)−e]
y''' = a sin(u) √(1−e²)
Find the true anomaly, θ, in radians
θ' = arctan(y'''/x''')
if x'''>0 and y'''>0 then θ=θ'
if x'''<0 then θ=θ'+π
if x'''>0 and y'''<0 then θ=θ'+2π
Rotate the canonical position vector into heliocentric ecliptic coordinates
x'' = x''' cos(ω) − y''' sin(ω)
y'' = x''' sin(ω) + y''' cos(ω)
x' = x''
y' = y'' cos(i)
z' = z'' sin(i)
x = x' cos(Ω) − y' sin(Ω)
y = x' sin(Ω) + y' cos(Ω)
z = z'
Find the heliocentric distance, r
r = √(x² + y² + z²)
Find the heliocentric ecliptic longitude, λ
λ' = arctan(y/x)
if x>0 and y>0 then λ=λ'
if x<0 then λ=λ+π
if x>0 and y<0 then λ=λ+2π
Find the heliocentric ecliptic latitude, β
β = arcsin(z/r)
Find the canonical velocity of the object in its own orbit
K = √{GM/[aU(1−e²)]}
Vx''' = −K sin(θ)
Vy''' = +K [e + cos(θ)]
Rotate the canonical velocity vector into heliocentric ecliptic coordinates
Vx'' = Vx''' cos(ω) − Vy''' sin(ω)
Vy'' = Vx''' sin(ω) + Vy''' cos(ω)
Vx' = Vx''
Vy' = Vy'' cos(i)
Vz' = Vz'' sin(i)
Vx = Vx' cos(Ω) − Vy' sin(Ω)
Vy = Vx' sin(Ω) + Vy' cos(Ω)
Vz = Vz'
Find the sun-relative speed, V
V = √(Vx² + Vy² + Vz²)
(elliptical orbit treatment ends here)
For HYPERBOLIC orbits...
Calculate the hyperbolic mean anomaly of the orbit,
M = 0.01720209895 (t−T) √[1/(−a)³]
Initialize the variable for the eccentric anomaly, u, as zero.
u = 0
U = 9
Refine the eccentric anomaly to the exact value
While |u−U|>0.0000000001 do
U = u
F₁ = e sinh(u) − u − M
F₂ = e cosh(u) − 1
F₃ = e sinh(u)
F₄ = e cosh(u)
D₁ = −F₁/F₂
D₂ = −F₁/(F₂ + D₁F₃/2)
D₃ = −F₁/(F₂ + D₁F₃/2 + F₄D₂²/6)
u = u + D₃
end while
Find the true anomaly, θ
if M≥0 then
θ = arccos{[e−cosh(u)]/[e cosh(u)−1]}
else
θ = −arccos{[e−cosh(u)]/[e cosh(u)−1]}
Find the heliocentric distance, r
r = a[1−e cosh(u)]
Calculate the canonical position of the object in its own orbit
x''' = r cos θ
y''' = r sin θ
Rotate the canonical position vector into heliocentric ecliptic coordinates
x'' = x''' cos(ω) − y''' sin(ω)
y'' = x''' sin(ω) + y''' cos(ω)
x' = x''
y' = y'' cos(i)
z' = z'' sin(i)
x = x' cos(Ω) − y' sin(Ω)
y = x' sin(Ω) + y' cos(Ω)
z = z'
Find the heliocentric distance, r
r = √(x² + y² + z²)
Find the heliocentric ecliptic longitude, λ
λ' = arctan(y/x)
if x>0 and y>0 then λ=λ'
if x<0 then λ=λ+π
if x>0 and y<0 then λ=λ+2π
Find the heliocentric ecliptic latitude, β
β = arcsin(z/r)
Find the canonical velocity of the object in its own orbit
K = √[GM/(−aU)]
Vx''' = (a/r) K sinh(u)
Vy''' = −(a/r) K √(e²−1) cosh(u)
Rotate the canonical velocity vector into heliocentric ecliptic coordinates
Vx'' = Vx''' cos(ω) − Vy''' sin(ω)
Vy'' = Vx''' sin(ω) + Vy''' cos(ω)
Vx' = Vx''
Vy' = Vy'' cos(i)
Vz' = Vz'' sin(i)
Vx = Vx' cos(Ω) − Vy' sin(Ω)
Vy = Vx' sin(Ω) + Vy' cos(Ω)
Vz = Vz'
Find the sun-relative speed, V
V = √(Vx² + Vy² + Vz²)
(hyperbolic orbit treatment ends here)