- #1

- 82

- 22

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

.

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.

P = 365.256898326 a^(1.5)

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 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.

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.

θ = 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.

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

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

z'' = z'''

x' = x''

y' = y'' cos i

z' = y'' sin i

x = x' cos Ω − y' sin Ω

y = x' sin Ω + y' cos Ω

z = z'

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.

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

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

Vz'' = Vz'''

Vx' = Vx''

Vy' = Vy'' cos i

Vz' = Vy'' sin i

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.

.

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.

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

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ᵢ₊₁

r = a (1 − e cosh u)

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

x''' = r cos θ

y''' = r sin θ

z''' = 0

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

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

z'' = z'''

x' = x''

y' = y'' cos i

z' = y'' sin i

x = x' cos Ω − y' sin Ω

y = x' sin Ω + y' cos Ω

z = z'

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.

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

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

Vz'' = Vz'''

Vx' = Vx''

Vy' = Vy'' cos i

Vz' = Vy'' sin i

Vx = Vx' cos Ω − Vy' sin Ω

Vy = Vx' sin Ω + Vy' cos Ω

Vz = Vz'

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.

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.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.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'''].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.

Last edited: