Can you find the Keplerian elements of this asteroid's orbit?
Find the semimajor axis, the eccentricity, the inclination to the ecliptic, the longitude of the ascending node, the argument of the perihelion, and the time of perihelion passage for this asteroid.
t₁ = JD 2457204.625
X⊕₁ = +0.155228396 AU
Y⊕₁ = −1.004732775 AU
Z⊕₁ = +0.00003295786 AU
α₁ = 20h 46m 57.02s
δ₁ = −27°41′33.9″
t₂ = JD 2457214.625
X⊕₂ = +0.319493277 AU
Y⊕₂ = −0.965116604 AU
Z⊕₂ = +0.0000311269 AU
α₂ = 20h 39m 57.10s
δ₂ = −28°47′21.5″
t₃ = JD 2457224.625
X⊕₃ = +0.4747795623 AU
Y⊕₃ = −0.8983801739 AU
Z⊕₃ = +0.00002841127 AU
α₃ = 20h 31m 22.81s
δ₃ = −29°49′22.7″
t₄ = JD 2457234.625
X⊕₄ = +0.616702829 AU
Y⊕₄ = −0.8063620175 AU
Z⊕₄ = +0.00002486325 AU
α₄ = 20h 22m 06.57s
δ₄ = −30°41′57.3″
tᵢ is the time of observation number i
X⊕ᵢ , Y⊕ ᵢ , Z⊕ᵢ is the position vector of Earth in heliocentric ecliptic coordinates
αᵢ is the geocentric right ascension of the asteroid at time tᵢ
δᵢ is the geocentric declination of the asteroid at time tᵢ
The source's author's name in Russian is Александр Дмитриевич Дубяго.
I had to trim away the results of the intermediate calculations because the length limit for answers is 10000 characters. However, you can find the completely worked out example problem at
- MorningfoxLv 74 months agoFavourite answer
Sounds like a standard application of Gauss’ Method. Yes, I could do it, but it’s a bit long for this space. Takes about 14 steps. Or, you could cheat, by recognizing that this is Ceres and looking up the elements for that.
- 4 months ago
Morningfox is right. The example that I used is Ceres. The method for orbit determination that I'm using here, though, is a variation on the three-observation method of Gauss that doesn't have a peculiar weakness that that method has.
Specifically, the three-point method of Gauss won't work if the apparent path of the asteroid across Earth's sky lacks sufficient curvature. That lack of curvature causes a certain trapezoidal volume to approach zero, and that volume is a factor in the denominator of certain coefficients in that 8-degree polynomial that you must solve for the asteroid's heliocentric and geocentric distances, which, initially, you don't know.
So here it is: a four-observation method.
Determining the Keplerian Elements of an Elliptical Orbit from Four Observations
Presented hereafter is a method for determining a preliminary heliocentric orbit from four geocentric directions of a sun-orbiting object at four distinct times of observation. I take this method after that presented in chapter six of "The Determination of Orbits," by A.D. Dubyago, who credits it to
The German mathematician Karl Friedrich Gauss
The German astronomer Julius Bauschinger
The Soviet celestial mechanic Mikhail Fedorovich Subbotin
The initial data
t₁, X⊕₁, Y⊕₁, Z⊕₁, α₁, δ₁
t₂, X⊕₂, Y⊕₂, Z⊕₂, α₂, δ₂
t₃, X⊕₃, Y⊕₃, Z⊕₃, α₃, δ₃
t₄, X⊕₄, Y⊕₄, Z⊕₄, α₄, δ₄
The times of observation, tᵢ, are given in Julian date format. The vectors [X⊕ᵢ,Y⊕ᵢ,Z⊕ᵢ] are the positions of the Earth in heliocentric ecliptic coordinates, with the components being in astronomical units. The αᵢ are the geocentric right ascensions for Ceres. The δᵢ are the geocentric declinations for Ceres.
The time intervals should be about 0.5% to 1% of the object's period (estimated as 8 to 16 days for a main belt asteroid), should be near opposition with the sun, but should NOT span an apside of the object's orbit. The precision in right ascension should be 0.01 seconds or better, and the precision in declination should be 0.1 arcseconds or better.
The Earth's orbital mean motion,
κ = 0.01720209895 radians/day
We will use a single value for the obliquity of the ecliptic to transform all four of the observation angles from celestial coordinates to ecliptic coordinates. First, we find the middle of the observation time window in units of 10000 years since 1 January 2000, and then we use that number to find the obliquity using the 10-degree polynomial fit of J. Laskar.
t = (t₁ + t₄)/2
T = (t − 2451545)/3652500
The obliquity in seconds of arc is
ε" = 84381.448 − 4680.93 T − 1.55 T² + 1999.25 T³ − 51.38 T⁴
− 249.67 T⁵ − 39.05 T⁶ + 7.12 T⁷ + 27.87 T⁸ + 5.79 T⁹ + 2.45 T¹⁰
The obliquity in radians is
ε = (π / 648000) ε"
Although Earth's axial tilt (i.e. the obliquity of the ecliptic) does change over time, it changes so slowly that the difference will almost always be negligible across the t₁ to t₄ time window. However, in the rare case when this isn't true, separate evaluations of the obliquity will have to be made for each time of observation.
The geocentric positions of the sun in celestial coordinates are (for i = 1 to 4)
Xᵢ = −X⊕ᵢ
Yᵢ = −Y⊕ᵢ cos ε + Z⊕ᵢ sin ε
Zᵢ = −Y⊕ᵢ sin ε − Z⊕ᵢ cos ε
The geocentric unit vectors in the direction of the target object, in celestial coordinates, are (for i = 1 to 4)
aᵢ = cos αᵢ cos δᵢ
bᵢ = sin αᵢ cos δᵢ
cᵢ = sin δᵢ
The squares of the Sun-Earth distances at times t₁ and t₄ are
R₁² = X₁² + Y₁² + Z₁²
R₄² = X₄² + Y₄² + Z₄²
The values of 2Rᵢ cos Θᵢ at times t₁ and t₄ are
2R₁ cos Θ₁ = −2 ( a₁X₁ + b₁Y₁ + c₁Z₁ )
2R₄ cos Θ₄ = −2 ( a₄X₄ + b₄Y₄ + c₄Z₄ )
where Θ is the supplementary angle to the sun-Earth-object angle.
We find some time differences:
τ₁ = κ (t₄−t₂)
τ₂ = κ (t₂−t₁)
τ₃ = κ (t₄−t₁)
τ₄ = κ (t₄−t₃)
τ₅ = κ (t₃−t₁)
We find this pair of determinants:
Φ = a₂b₄−b₂a₄
φ = a₃b₄−b₃a₄
We calculate some intermediate quantities:
A = ( a₁b₂ − b₁a₂ ) / Φ
B = ( a₂Y₁ − b₂X₁ ) / Φ
C = ( b₂X₂ − a₂Y₂ ) / Φ
D = ( a₂Y₄ − b₂X₄ ) / Φ
A' = ( a₁b₃ − b₁a₃ ) / φ
B' = ( a₃Y₁ − b₃X₁ ) / φ
C' = ( b₃X₃ − a₃Y₃ ) / φ
D' = ( a₃Y₄ − b₃X₄ ) / φ
And then more intermediate quantities:
E = τ₁/τ₂
F = (4/3) τ₁τ₃
G = AE
H = F(A−G)
I = 4Aτ₁²
K = E (B + C) + C + D
L = F (B − C + D − K)
M = 4 (Bτ₁² + τ₁τ₂C)
E' = τ₄/τ₅
F' = (4/3) τ₄τ₃
G' = A'E'
H' = F'(A'−G')
I' = 4A'τ₄²
K' = E' (B' + C') + C' + D'
L' = F' (B' − C' + D' − K')
M' = 4 (B'τ₄² + τ₄τ₅C')
Make initial guesses for the sun-object distance r₁ at time t₁ and for the sun-object distance r₄ at time t₄. For main belt asteroids, a reasonable initial guess for both times is 2.75 AU. Then use the loop below to converge, by successive approximations, to the true sun-object distances, r, and for the Earth-object distances, ρ, distances at times t₁° and t₄° (i.e. the times for the first and fourth observations, corrected for the speed of light travel time).
O = 9.999e+99
N = r₁ + r₄
while |N−O|/N > 1ᴇ-11 do
ξ = (r₁ + r₄)⁻³
η = (r₄ − r₁) / (r₁ + r₄)
P = G + ξH + ηξI
Q = K + ξL + ηξM
P' = G' + ξH' + ηξI'
Q' = K' + ξL' + ηξM'
ρ₁ = (Q'−Q)/(P−P')
ρ₄ = Pρ₁ + Q
r₁ = √[R₁² + (2R₁cos Θ₁)ρ₁ + ρ₁²]
r₄ = √[R₄² + (2R₄cos Θ₄)ρ₄ + ρ₄²]
O = N
N = r₁ + r₄
The positions of the object at times t₁ and t₄ in geocentric celestial coordinates are
x₁ = a₁ρ₁ − X₁
y₁ = b₁ρ₁ − Y₁
z₁ = c₁ρ₁ − Z₁
x₄ = a₄ρ₄ − X₄
y₄ = b₄ρ₄ − Y₄
z₄ = c₄ρ₄ − Z₄
The reciprocal of the speed of light
ç = 0.00577551833 days/AU
The sun's gravitational parameter
μ = 1.32712440018ᴇ20 m³ sec⁻²
The conversion factor from AU to meters is
U = 1.495978707ᴇ11
The conversion factor from AU/day to m/sec is
β = 1731456.8368
Correcting times of observation for planetary abberation.
t₁° = t₁ − çρ₁
t₄° = t₄ − çρ₄
The nominal time associated with the forthcoming state vector is
t₀ = ½ (t₁° + t₄°)
The nominal heliocentric distance of the object at time t₀ is
r₀ = ½ (r₁ + r₄)
Find the heliocentric position vector [x',y',z'] for the object at time t₀ in celestial coordinates.
x" = ½ (x₁ + x₄)
y" = ½ (y₁ + y₄)
z" = ½ (z₁ + z₄)
r" = √[(x")²+(y")²+(z")²]
x' = (r₀/r") U x"
y' = (r₀/r") U y"
z' = (r₀/r") U z"
Find the sun-relative velocity vector for the object at time t₀ in celestial coordinates.
S = √[ (x₄ − x'/U)² + (y₄ − y'/U)² + (z₄ − z'/U)² ]
s = √[ (x'/U − x₁)² + (y'/U − y₁)² + (z'/U − z₁)² ]
Ψ = S + s
ψ = √[(x₄−x₁)² + (y₄−y₁)² + (z₄−z₁)²]
Vx' = β (Ψ/ψ) (x₄−x₁) / (t₄°−t₁°)
Vy' = β (Ψ/ψ) (y₄−y₁) / (t₄°−t₁°)
Vz' = β (Ψ/ψ) (z₄−z₁) / (t₄°−t₁°)
The object's position in heliocentric ecliptic coordinates
x₀ = x'
y₀ = y' cos ε + z' sin ε
z₀ = −y' sin ε + z' cos ε
The object's sun-relative velocity in ecliptic coordinates
Vx₀ = Vx'
Vy₀ = Vy' cos ε + Vz' sin ε
Vz₀ = −Vy' sin ε + Vz' cos ε
The object's speed relative to the sun
V₀ = √[(Vx₀)² + (Vy₀)² + (Vz₀)²]
The semimajor axis of the object's orbit, in AU
a = (2/r₀ − V₀²/μ)⁻¹ / U
The angular momentum per unit mass in the object's orbit
hx = y₀ Vz₀ − z₀ Vy₀
hy = z₀ Vx₀ − x₀ Vz₀
hz = x₀ Vy₀ − y₀ Vx₀
h = √[(hx)² + (hy)² + (hz)²]
The eccentricity of the object's orbit
e = √[1 − h²/(aμU)]
The inclination of the object's orbit
i = arccos(hz/h)
The longitude of the ascending node of the object's orbit
Ω' = arctan(−hx/hy)
if hy>0 then Ω = Ω' + π
If hy<0 and hx<0 then Ω = Ω' + 2π
The true anomaly at time t₀
sin θ₀ = h ( x₀ Vx₀ + y₀ Vy₀ + z₀ Vz₀ ) / (r₀μ)
cos θ₀ = h²/(r₀μ) − 1
θ₀' = arctan( sin θ₀ / cos θ₀ )
If cos θ₀ < 0 then θ₀ = θ₀' + π
If cos θ₀ > 0 and sin θ₀ < 0 then θ₀ = θ₀' + 2π
The sum of the true anomaly at time t₀ and the argument of the perihelion of the object's orbit
sin(θ₀+ω) = z₀ / (r₀ sin i)
cos(θ₀+ω) = ( x₀ cos Ω + y₀ sin Ω ) / r₀
(θ₀+ω)' = arctan[ sin(θ₀+ω) / cos(θ₀+ω) ]
If cos(θ₀+ω) < 0 then θ₀ = θ₀' + π
If cos(θ₀+ω) > 0 and sin(θ₀+ω) < 0 then (θ₀+ω) = (θ₀+ω)' + 2π
The argument of the perihelion of the object's orbit
ω' = (θ₀+ω) − θ₀
if ω'<0 then ω=ω'+2π else ω=ω'
The eccentric anomaly of the object at time t₀
cos u₀ = 1 − r₀/(aU)
sin u₀ = (x₀ Vx₀ + y₀ Vy₀ + z₀ Vz₀) / √(aμU)
u₀' = arctan( sin u₀ / cos u₀ )
If cos u₀ < 0 then u₀ = u₀' + π
If cos u₀ > 0 and sin u₀ < 0 then u₀ = u₀' + 2π
The mean anomaly of the object at time t₀
M₀ = u₀ − e sin u₀
The period of the object's orbit in days
P = 365.256898326 a¹·⁵
The object's time of perihelion passage
T = t₀ − PM₀/(2π)
These are the results for the example problem.
Orbital elements for Ceres (as calculated)
a = 2.76694735 AU
e = 0.076026341
i = 10.5918141°
Ω = 80.3183813°
ω = 72.6265868°
T₀ = JD 2456552.87
Elements for Ceres from JPL Horizons for JD 2457219.5 (16 July 2015)
a = 2.76800849 AU
e = 0.075773402
i = 10.5922177°
Ω = 80.3268351°
ω = 72.6626266°
T = JD 2456552.64Source(s): The Determination of Orbits, by A.D. Dubyago