Position and Velocity in Heliocentric Ecliptic Coordinates

[Jenab2]

There have been many questions on this forum about celestial mechanics in general, and concerning position and velocity in an orbit in particular. So I offer this post as a summary and reference.

Here's a method for finding heliocentric position and sun-relative velocity in ecliptic coordinates from the Keplerian orbital elements (a, e, i, Ω, ω, T) at a time of interest, t, for both elliptical and hyperbolic orbits relative to the sun.

a : semimajor axis
* here usually in astronomical units

e : eccentricity

i : inclination to the ecliptic
* entered in degrees, usually converted to radians

Ω : longitude of the ascending node
* entered in degrees, usually converted to radians
* Ω is an angle, subtended at the sun, measured in the ecliptic plane in the counterclockwise direction from the perspective of someone positioned far along the North Ecliptic Pole, from a ray extending toward the Vernal Equinox to another ray extending toward the place where the orbit crosses the ecliptic going from south to north.

ω : argument of the perihelion
* entered in degrees, usually converted to radians
* ω is an angle, subtended at the sun, measured in the plane of the orbit, in the angular direction of motion, from a ray extending toward the ascending node to another ray extending toward the perihelion of the orbit.

T : time of perihelion passage
* here usually in Julian date, with differences in time being usually in days

.

For ELLIPTICAL ORBITS.

An elliptical orbit has a positive semimajor axis and an eccentricity strictly between zero and one. It is a bound orbit with a period. The planets, the asteroids, and some comets are in elliptical orbits.

The period of the orbit, P.

P = 365.256898326 a^(1.5)

The mean anomaly, m.

m''' = (t−T)/P
m'' = m''' − integer(m''')
If m'' is less than 0, then m'=m''+1 else m'=m''
m = 2πm'

The eccentric anomaly, u.

The mean anomaly must be in radians when finding the eccentric anomaly.

u₀ = m + (e − e³/8 + e⁵/192) sin(m) + (e²/2 − e⁴/6) sin(2m) + (3e³/8 − 27e⁵/128) sin(3m) + (e⁴/3) sin(4m)

i = 0

Repeat
i = i+1
E = uᵢ − e sin uᵢ − m
F = 1 − e cos uᵢ
G = e sin uᵢ
H = e cos uᵢ
A = −E/F
B = −E/(F + ½ AG)
C = −E/(F + ½ AG + ⅙ B²H)
uᵢ₊₁ = uᵢ + C
Until |uᵢ₊₁−uᵢ| < 1e-12

u = uᵢ₊₁

The eccentric anomaly, u, will be in radians.
* u is an angle, subtended at the geometric center of the elliptical orbit (not at the sun), measured in the plane of the orbit in the direction of motion, from a ray extending through the orbit's perihelion, to another ray extending through a point on a circle that circumscribes the orbital ellipse which is in the +y direction of the position of the body in its orbit at time t.

The canonical heliocentric position vector [x''',y''',z'''].

x''' = a (cos u − e)
y''' = a sin u √(1− e²)
z''' = 0

The canonical coordinate system has the sun at the origin, the XY plane in the plane of the orbit, with the orbit's perihelion being on the +x axis, and the +z axis such that an observer somewhere along it would observe a counterclockwise direction of motion.

The true anomaly, θ.

θ = arctan( y''' , x''' )

* θ is an angle, subtended at the sun, measured in the plane of the orbit, in the angular direction of motion, from a ray extending toward the perihelion of the orbit to another ray extending toward the orbiting body's position at time t.

Rotate the coordinates from canonical to ecliptic.

x'' = x''' cos ω − y''' sin ω
y'' = x''' sin ω + y''' cos ω
z'' = z'''

x' = x''
y' = y'' cos i
z' = y'' sin i

Position in heliocentric ecliptic coordinates [x,y,z].

x = x' cos Ω − y' sin Ω
y = x' sin Ω + y' cos Ω
z = z'

The canonical heliocentric velocity vector [Vx''',Vy''',Vz'''].

Vx''' = (−29784.6918 m/s) √{1/[a(1−e²)]} sin θ
Vy''' = (29784.6918 m/s) √{1/[a(1−e²)]} (e + cos θ)
Vz''' = 0

Enter the semimajor axis, a, in astronomical units. The dimension of length has already been extracted from the square root.

Rotate the coordinates from canonical to ecliptic.

Vx'' = Vx''' cos ω − Vy''' sin ω
Vy'' = Vx''' sin ω + Vy''' cos ω
Vz'' = Vz'''

Vx' = Vx''
Vy' = Vy'' cos i
Vz' = Vy'' sin i

Sun-relative velocity in ecliptic coordinates [Vx,Vy,Vz].

Vx = Vx' cos Ω − Vy' sin Ω
Vy = Vx' sin Ω + Vy' cos Ω
Vz = Vz'

Sun-relative speed, v.
v = √[(Vx)²+(Vy)²+(Vz)²]
v = (29784.6918 m/s) √(2/r−1/a)

Enter a,r in astronomical units. The dimension of length has already been extracted from the square root.

.

For HYPERBOLIC ORBITS.

A hyperbolic orbit has a negative semimajor axis and an eccentricity greater than one. It is unbound and not periodic, but one-pass only. Some comets are in hyperbolic orbits.

The hyperbolic mean anomaly, m.

m = 0.01720209895 (t−T) √[1/(−a)³]

The hyperbolic eccentric anomaly, u.

The hyperbolic mean anomaly must be in radians when finding the hyperbolic eccentric anomaly.

u₀ = 0

i = 0

Repeat
i = i+1
E = e sinh uᵢ − uᵢ − m
F = e cosh uᵢ − 1
G = e sinh uᵢ
H = e cosh uᵢ
A = −E/F
B = −E/(F + ½ AG)
C = −E/(F + ½ AG + ⅙ B²H)
uᵢ₊₁ = uᵢ + C
Until |uᵢ₊₁−uᵢ| < 1e-12

u = uᵢ₊₁

Distance from the sun, r.

r = a (1 − e cosh u)

True anomaly, θ.

θ = arccos{ (e − cosh u) / (e cosh u − 1) }

The canonical heliocentric position vector [x''',y''',z'''].

x''' = r cos θ
y''' = r sin θ
z''' = 0

Rotate the coordinates from canonical to ecliptic.

x'' = x''' cos ω − y''' sin ω
y'' = x''' sin ω + y''' cos ω
z'' = z'''

x' = x''
y' = y'' cos i
z' = y'' sin i

Position in heliocentric ecliptic coordinates [x,y,z].

x = x' cos Ω − y' sin Ω
y = x' sin Ω + y' cos Ω
z = z'

The canonical heliocentric velocity vector [Vx''',Vy''',Vz''').

Vx''' = (29784.6918 m/s) (a/r) √(−1/a) sinh u
Vy''' = (−29784.6918 m/s) (a/r) √[(1−e²)/a] cosh u
Vz''' = 0

Enter a,r in astronomical units. The dimension of length has already been extracted from the square root.

Rotate the coordinates from canonical to ecliptic.

Vx'' = Vx''' cos ω − Vy''' sin ω
Vy'' = Vx''' sin ω + Vy''' cos ω
Vz'' = Vz'''

Vx' = Vx''
Vy' = Vy'' cos i
Vz' = Vy'' sin i

Sun-relative velocity in ecliptic coordinates [Vx,Vy,Vz].

Vx = Vx' cos Ω − Vy' sin Ω
Vy = Vx' sin Ω + Vy' cos Ω
Vz = Vz'

Sun-relative speed, v.

v = √[(Vx)²+(Vy)²+(Vz)²]
v = (29784.6918 m/s) √(2/r−1/a)

Enter a,r in astronomical units. The dimension of length has already been extracted from the square root.

 

[Jenab2]

You can also go the other way.

Beginning with a position vector and a velocity vector in heliocentric ecliptic coordinates, at a time t, for an object in orbit around the sun, you can derive its Keplerian orbital elements (a, e, i, Ω, ω, T). The combination of the vectors of position and of velocity in an orbit is sometimes called the "state vector."

If you are given the position in spherical ecliptic coordinates (r = distance, λ = ecliptic longitude, β = ecliptic latitude), then do this to get the position in rectangular ecliptic coordinates:

x = r cos λ cos β
y = r sin λ cos β
z = r sin β

Let's assume that x, y, and z are in meters, and that Vx, Vy, Vz are the components of the sun-relative velocity in ecliptic coordinates in meters per second.

Let t be the time-of-interest in Julian date.

The solar gravitational parameter.

GM = 1.32712440018e20 m³ sec⁻²

The heliocentric distance, r (in meters), of the object at time t.

r = √( x² + y² + z² )

The sun-relative speed, v (in m/sec), of the object at time t.

v = √[ (Vx)² + (Vy)² + (Vz)² ]

The semimajor axis, a, of the orbit, meters.

a = 1 / { 2/r − v²/(GM) }

You can convert the semimajor axis to astronomical units.
1 AU = 1.49597870700e11 meters.

The angular momentum per unit mass (m²/sec), h, in the orbit.

hx = y Vz − z Vy
hy = z Vx − x Vz
hz = x Vy − y Vx

h = √[ (hx)² + (hy)² + (hz)² ]

The eccentricity, e, of the orbit.

e = √[ 1 − h² / (GMa) ]

The inclination, i, of the orbit.

i = ArcCos( hz / h )

The longitude of the ascending node, Ω, of the orbit.

Ω = Arctan ( hx , −hy )

Definition of the two-argument Arctan function (result in radians).

atn(z) = single argument Arctan function of the argument z.

Function Arctan( y , x )
if x = 0 and y > 0 then angle = +π/2
if x = 0 and y = 0 then angle = 0
if x = 0 and y < 0 then angle = −π/2
if x < 0 then angle = atn(y/x) + π
if x > 0 and y > 0 then angle = atn(y/x)
if x > 0 and y < 0 then angle = atn(y/x) + 2π
Arctan = angle

Note: the math to this point works for either elliptical or hyperbolic orbits. However, some of the stuff to follow works only for elliptical orbits. You'll know that the orbit is a hyperbola if (1) the semimajor axis is negative and (2) the eccentricity is greater than one.

The true anomaly, θ, of the object in the orbit at time t.

cos θ = h²/(rGM) − 1
sin θ = h (x Vx + y Vy + z Vz) / (rGM)
θ = Arctan ( sin θ , cos θ )

Argument of the perihelion, ω, of the orbit.

cos ω'' = (x cos Ω + y sin Ω) / r
If sin i = 0 then sin ω'' = (y cos Ω − x sin Ω) / r
If sin i ≠ 0 then sin ω'' = z / (r sin i)
ω'' = Arctan( sin ω'' , cos ω'' )
ω' = ω'' − θ
If ω' > 0 then ω = ω'
If ω' < 0 then ω = ω' + 2π

The eccentric anomaly, u, of the object in the orbit at time t.

sin u = (r/a) sin θ / √(1−e²)
cos u = (r/a) cos θ + e
u = Arctan( sin u , cos u )

The mean anomaly, m, of the object in the orbit at time t.

m = u − e sin u

Note: u must be entered in radians, and m will return in radians.

The period of the orbit, P, days.

P = (π / 43200) √[ a³/(GM) ]

The mean motion, μ, in the orbit, radians per day.

μ = 2π / P

The time of perihelion passage, T, of the object in the orbit, Julian Date.

T = t − m/μ

The Keplerian orbital elements.

[ a, e, i, Ω, ω, T ]

 

[Jenab2]

Converting from heliocentric ecliptic coordinates to geocentric celestial coordinates.

Given the heliocentric positions of Earth [x₀,y₀,z₀] and of another planet in ecliptic coordinates [x,y,z] for the same moment of time, t, expressed as a Julian Date, we will find the geocentric position of the planet in celestial coordinates; i.e., we want the planet's right ascension and declination.

dx = x − x₀
dy = y − y₀
dz = z − z₀

We find the obliquity of Earth, ε, at time t.

T = t − 2451545.0

ε = 23.4392911° − 3.562266e-7 T − 1.22848e-16 T² + 1.03353e-20 T³

This equation gives the obliquity in degrees.

dx' = dx
dy' = dy cos ε − dz sin ε
dz' = dy sin ε + dz cos ε

The distance between Earth and the planet, dr, at time t is found from

dr = √{ (dx)² + (dy)² + (dz)² }

The planet's geocentric right ascension in decimal hours, α, at time t is found from

α = (12/π) Arctan( dy' , dx' )

The planet's geocentric declination in decimal degrees, δ, at time t is found from

δ = (180/π) Arcsin( dz' / dr )

 

[Jenab2]

Converting Calendar Date into Julian Date.

Year = the year as a four digit integer
Month = the month, an integer from 1-12
Day = the day of the month, an integer

Q₁ = integer( ( Month − 14 ) / 12 )
Q₂ = integer( ( 1461 ( Year + 4800 + Q₁ ) ) / 4 )
Q₃ = integer( ( 367 ( Month − 2 − 12 Q₁ ) ) / 12 )
Q₄ = integer( ( Year + 4900 + Q₁ ) / 100 )
Q₅ = integer( ( 3 Q₄ ) / 4 )

J = Q₂ + Q₃ − Q₅ + Day − 32075.5

J : Julian Date of midnight, Greenwich Mean Time (GMT), beginning the calendar date of Day/Month/Year.

Hour = hours in the GMT of the time of interest
Min = minutes in the GMT of the time of interest
Sec = seconds in the GMT of the time of interest

t = J + Hour/24 + Min/1440 + Sec/86400

t : Julian date of the time of interest