# A Determination of a Preliminary Orbit from Four Observations

## Recommended Posts

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.

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¹⁰

ε = (π / 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₄
endwhile

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π)

Example problem. Find the orbit of Ceres.

Observation #1
t₁ = JD 2457204.625
X⊕₁ = +0.155228396 AU
Y⊕₁ = −1.004732775 AU
Z⊕₁ = +0.00003295786 AU
α₁ = 20h 46m 57.02s
δ₁ = −27°41′33.9″

Observation #2
t₂ = JD 2457214.625
X⊕₂ = +0.319493277 AU
Y⊕₂ = −0.965116604 AU
Z⊕₂ = +0.0000311269 AU
α₂ = 20h 39m 57.10s
δ₂ = −28°47′21.5″

Observation #3
t₃ = JD 2457224.625
X⊕₃ = +0.4747795623 AU
Y⊕₃ = −0.8983801739 AU
Z⊕₃ = +0.00002841127 AU
α₃ = 20h 31m 22.81s
δ₃ = −29°49′22.7″

Observation #4
t₄ = JD 2457234.625
X⊕₄ = +0.616702829 AU
Y⊕₄ = −0.8063620175 AU
Z⊕₄ = +0.00002486325 AU
α₄ = 20h 22m 06.57s
δ₄ = −30°41′57.3″

Converting the observation angles to radians:

α₁ = 5.44084723
δ₁ = −0.483329666

α₂ = 5.41030979
δ₂ = −0.502468171

α₃ = 5.37290956
δ₃ = −0.520509058

α₄ = 5.33245865
δ₄ = −0.53580299

Continuing through the procedure,

t = 2457219.63
T = 0.001553628

X₁ = −0.155228396
Y₁ = +0.921851498
Z₁ = +0.399597005

X₂ = −0.319493277
Y₂ = +0.885503087
Z₂ = +0.383841559

X₃ = −0.474779562
Y₃ = +0.824271594
Z₃ = +0.357299982

X₄ = −0.616702829
Y₄ = +0.739843884
Z₄ = +0.320703493

a₁ = +0.589463371
b₁ = −0.660726081
c₁ = −0.464730008

a₂ = +0.563195268
b₂ = −0.671477524
c₂ = −0.4815901

a₃ = +0.532276132
b₃ = −0.685093499
c₃ = −0.497321844

a₄ = +0.49965704
b₄ = −0.69978587
c₄ = −0.510531662

R₁² = 1.03358381
R₄² = 1.03054208

2R₁ cos θ₁ = 1.772595
2R₄ cos θ₄ = 1.97920299

τ₁ = 0.344041979
τ₂ = 0.17202099
τ₃ = 0.516062969
τ₄ = 0.17202099
τ₅ = 0.344041979

Φ = −0.058607619
φ = −0.030167527

A = +0.404275125
B = −7.08013785
C = +4.84883362
D = −0.043927505

A' = +1.72864027
B' = −12.7399767
C' = +3.76138574
D' = +0.951283093

E = +2
F = +0.236729767
G = +0.80855025
H = −0.095703956
I = +0.191407912
K = +0.342297675
L = −2.91537363
M = −2.20429551

E' = +0.5
F' = +0.118364883
G' = +0.864320133
H' = +0.102305152
I' = +0.204610303
K' = +0.223373365
L' = −1.86702289
M' = −0.617533885

r₁ = 2.75 (initial guess)
r₄ = 2.75 (initial guess)

1st approximation
ρ₁ = 1.97723208
ρ₄ = 1.9223289
r₁ = 2.90652064
r₄ = 2.92071388

2nd approximation
ρ₁ = 2.001149
ρ₄ = 1.94460337
r₁ = 2.93008666
r₄ = 2.94292188

3rd approximation
ρ₁ = 2.00417695
ρ₄ = 1.94741724
r₁ = 2.93307059
r₄ = 2.94572742

4th approximation
ρ₁ = 2.00455348
ρ₄ = 1.94776713
r₁ = 2.93344165
r₄ = 2.94607627

5th approximation
ρ₁ = 2.0046002
ρ₄ = 1.94781054
r₁ = 2.93348769
r₄ = 2.94611956

6th approximation
ρ₁ = 2.00460599
ρ₄ = 1.94781593
r₁ = 2.9334934
r₄ = 2.94612492

7th approximation
ρ₁ = 2.00460671
ρ₄ = 1.9478166
r₁ = 2.93349411
r₄ = 2.94612559

8th approximation
ρ₁ = 2.0046068
ρ₄ = 1.94781668
r₁ = 2.9334942
r₄ = 2.94612567

9th approximation
ρ₁ = 2.00460681
ρ₄ = 1.94781669
r₁ = 2.93349421
r₄ = 2.94612568

10th approximation (final, converged)
ρ₁ = 2.00460681
ρ₄ = 1.94781669
r₁ = 2.93349421
r₄ = 2.94612568

HEC positions in AU at t₁ & t₄

x₁ = +1.33687069
y₁ = −2.2463475
z₁ = −1.33119794

x₄ = +1.58994315
y₄ = −2.10289848
z₄ = −1.31512559

The reciprocal of the speed of light
ç = 0.00577551833 days/AU

Aberration corrections to time

çρ₁ = 0.011577643 days
t₁° = 2457204.61

çρ₄ = 0.011249651 days
t₄° = 2457234.61

Epoch of state vector
t₀ = 2457219.61 JD

The object's state vector in heliocentric ecliptic coordinates

x₀ = +1.46520344 AU
y₀ = −2.52458426 AU
z₀ = −0.349479243 AU
Vx₀ = +14610.4367 m/s
Vy₀ = +7967.42879 m/s
Vz₀ = −2442.63758 m/s

The object's distance from the sun at t₀
r₀ = 2.93980995 AU

Sun-relative speed
V₀ = 16819.9661 m/s

The sun's gravitational parameter
μ = 1.32712440018ᴇ20 m³ sec⁻²

The semimajor axis of the object's orbit
a = 2.76694735 AU

The angular momentum per unit mass in the object's orbit

hx = +1.3390648ᴇ15 m²/sec
hy = −2.2844842ᴇ14 m²/sec
hz = +7.26435028ᴇ15 m²/sec

h = 7.39026848ᴇ15 m²/sec

The eccentricity of the object's orbit
e = 0.076026341

The inclination of the object's orbit
i = 10.5918141°

The longitude of the ascending node of the object's orbit
Ω = 80.3183813°

The true anomaly at time t₀
θ₀ = 147.669798°

The sum of the true anomaly at time t₀ and the argument of the perihelion of the object's orbit
(θ₀+ω) = 220.296384°

The argument of the perihelion of the object's orbit
ω = 72.6265867°

The eccentric anomaly of the object at time t₀
u₀ = 145.259666°

The mean anomaly of the object at time t₀
M₀ = 142.777370°

The period of the object's orbit
P = 1681.12408 days

Times of perihelion passage

T₀ = 2456552.87
T₁ = 2458234.01 = T₀+P

A summary of the calculation results and a comparison with JPL's numbers

Orbital elements (as calculated)

a = 2.76694735 AU
e = 0.076026341
i = 10.5918141°
Ω = 80.3183813°
ω = 72.6265868°
T₀ = JD 2456552.87
T₁ = JD 2458234.01

Orbital elements (From the JPL Small-Body Database)

a = 2.76916515 AU
e = 0.076009027
i = 10.5940672°
Ω = 80.3055309°
ω = 73.5976947°
T = JD 2458238.75

Edited by David Sims

## Create an account

Register a new account

• ### Similar Content

• #### Observing partners!

By maw lod qan,

• 0 replies
• 66 views
• #### MAY 15TH - WL - AR3007 PLUS FOUR, SN1O5

By paulastro,

• 0 replies
• 67 views
• #### Adapters four thirds and micro four thirds lenses?

By PatG,

• 7 replies
• 396 views
• #### Mars & Saturn: Images Four Years Apart

By orion25,

• 8 replies
• 327 views
• #### APRIL 17TH - AR2993, AR2994 PLUS FOUR

By paulastro,

• 0 replies
• 132 views
×
×
• Create New...