Determining an Orbit by the Method of Gauss

Racist Celestial Mechanics, Lecture 5

Revised: 13 August 2004

Note: The "symbol" font is REQUIRED for the correct presentation of the equations.

When an asteroid is discovered, the first determination of its orbit is often made without any a priori knowledge of its distance from Earth or its position with respect to the sun. What is known initially is only the position of Earth with respect to the sun, and the right ascension and declination of the asteroid at three times of observation.

Karl Friedrich Gauss, a smart German fellow, figured out how to approximate the geocentric and heliocentric distances of the asteroid, based on the apparent curvature of the orbit along its observed path.

My presentation of this "method of Gauss" for determining the orbital elements of asteroids from three angular positions at three known times was adapted from the fifth chapter of The Determination of Orbits, by A. D. Dubyago, as translated from Russian into English by the RAND Corporation. Like Gauss, Dubyago was a White man. I doubt that many Blacks could even understand Dubyago's explanation of Gauss's method. But if you happen to be a nigger, consider yourself free to make the attempt anyway.

Most of the work will be done in the celestial coordinate system, rather than the ecliptic system that I usually prefer. We'll make the rotation to ecliptic coordinates near the end.

Initial Data.

For i=1 to 3

The time of observation: t_i
The position of the sun: Xs_i, Ys_i, Zs_i
The anglular position of the asteroid: a_i, d_i

The position of the sun is in local celestial coordinates, not geocentric celestial coordinates. To get Rs, you would need to subtract the your geocentric position on Earth's surface (in celestial coordinates) from the geocentric position of the sun (in celestial coordinates). If you are White, you should have no trouble handling that.

Some Handy Constants.

k = 0.01720209895 (the mean motion of Earth in radians per day)
A = 0.00577551833 (the number of days required for light to travel 1 AU)

Unit Vector to the Asteroid.

For i=1 to 3

a_i = cos a_i cos d_i
b_i = sin a_i cos d_i
c_i = sin d_i

The position vectors are in local celestial coordinates.

The First Approximation.

t₁ = k { t₃ - t₂ }
t₂ = k { t₃ - t₁ }
t₃ = k { t₂ - t₁ }

n^o₁ = t₁ / t₂
n^o₃ = t₃ / t₂

n₁ = t₁ t₃ {1 + n^o₁} / 6
n₃ = t₁ t₃ {1 + n^o₃} / 6

D = a₂ { b₃ c₁ - b₁ c₃ } + b₂ { c₃ a₁ - c₁ a₃ } + c₂ { a₃ b₁ - a₁ b₃ }

For i=1 to 3

d_i = Xs_i { b₃ c₁ - b₁ c₃ } + Ys_i { c₃ a₁ - c₁ a₃ } + Zs_i { a₃ b₁ - a₁ b₃ }

Rs_i = { Xs_i² + Ys_i² + Zs_i² }^1/2

q_i = -2 { a_i Xs_i + b_i Ys_i + c_i Zs_i }

K₀ = { d₂ - d₁ n^o₁ - d₃ n^o₃ } / D

L₀ = { d₁ n₁ + d₃ n₃ } / D

r₂ : the distance from the sun to the asteroid at time 2.
r₂ : the distance from Earth to the asteroid at time 2.

We must solve these two equations simultaneously for r₂ and r₂:

r₂ = K₀ - L₀ / r₂³

r₂ = { Rs₂² + q₂ r₂ + r₂² }^1/2

This leads to an 8th degree polynomial in r₂.

f{r₂} = r₂⁸ - { Rs₂² + q₂ K₀ + K₀² } r₂⁶ + L₀ { q₂ + 2 K₀ } r₂³ - L₀²

We are interested in the roots of the function f{r₂}. We will use Newton's method to get the answer. The first derivative is

f'{r₂} = 8 r₂⁷ - 6 { Rs₂² + q₂ K₀ + K₀² } r₂⁵ + 3 L₀ { q₂ + 2 K₀ } r₂²

As an initial guess at r₂, we'll set it equal to one.

r_2,(j=0) = 1.0

Repeat over j

r_2,j+1 = r_2,j - f{r_2,j} / f'{r_2,j}

Assign to r₂ the converged value from the loop.

r₂ = r_2,j+1

r₂ = K₀ - L₀ / r₂³

n₁ = n^o₁ + n₁ / r₂³
n₃ = n^o₃ + n₃ / r₂³

Q₁ = n₁ Xs₁ - Xs₂ + n₃ Xs₃ + a₂ r₂
Q₂ = n₁ Ys₁ - Ys₂ + n₃ Ys₃ + b₂ r₂
Q₃ = n₁ Zs₁ - Zs₂ + n₃ Zs₃ + c₂ r₂
Q₄ = a₁ - a₃ c₁ / c₃
Q₅ = b₁ - b₃ a₁ / a₃
Q₆ = b₃ - b₁ a₃ / a₁
Q₇ = a₃ - a₁ c₃ / c₁

r₁ = { [ Q₁ - a₃ Q₃ / c₃ ] / [ n₁ Q₄ ] + [ Q₂ - b₃ Q₁ / a₃ ] / [ n₁ Q₅ ] } / 2
r₃ = { [ Q₂ - b₁ Q₁ / a₁ ] / [ n₃ Q₆ ] + [ Q₁ - a₁ Q₃ / c₁ ] / [ n₃ Q₇ ] } / 2

r₁ = { Rs₁² + q₁ r₁ + r₁² }^1/2
r₃ = { Rs₃² + q₃ r₃ + r₃² }^1/2

Correction for Planetary Abberation.

Light doesn't travel instantaneously. We must correct the initial times for the amount of time it took light to travel from the asteroid to our eyes. This is done only once, after the first approximation, but before the 2nd approximation. It will not be done between any other successive approximations.

For i=1 to 3

t^o_i = t_i - A r_i

The values t^o_i are the adjusted times of observation.

t₁ = k { t^o₃ - t^o₂ }
t₂ = k { t^o₃ - t^o₁ }
t₃ = k { t^o₂ - t^o₁ }

n^o₁ = t₁ / t₂
n^o₃ = t₃ / t₂

Recursive Procedure for Successive Approximations.

For i=1 to 3

x_i = a_i r_i - Xs_i
y_i = b_i r_i - Ys_i
z_i = c_i r_i - Zs_i

h₁ = t₁² / { K₁² [ K₁/3 + ( r₂+r₃ )/2 ] }
h₂ = t₂² / { K₂² [ K₂/3 + ( r₁+r₃ )/2 ] }
h₃ = t₃² / { K₃² [ K₃/3 + ( r₁+r₂ )/2 ] }

For i=1 to 3

j₁ = (11/9) h_i

j₂ = j₁

Repeat

j₃ = j₂

j₂ = j₁ / ( 1 + j₂)

Until | ( j₂ - j₃ ) / j₃ | < 10^-12

y_i = 1 + (10/11) j₂

n₁ = n^o₁ r₂³ { y₂/y₁ - 1 }
n₃ = n^o₃ r₂³ { y₂/y₃ - 1 }

K₀ = { d₂ - d₁ n^o₁ - d₃ n^o₃ } / D

L₀ = { d₁ n₁ + d₃ n₃ } / D

f{r₂} = r₂⁸ - { Rs₂² + q₂ K₀ + K₀² } r₂⁶ + L₀ { q₂ + 2 K₀ } r₂³ - L₀²

f'{r₂} = 8 r₂⁷ - 6 { Rs₂² + q₂ K₀ + K₀² } r₂⁵ + 3 L₀ { q₂ + 2 K₀ } r₂²

r_2,(j=0) = 1.0

Repeat over j

r_2,j+1 = r_2,j - f{r_2,j} / f'{r_2,j}

r₂ = r_2,j+1

r₂ = K₀ - L₀ / r₂³

Q₁ = n₁ Xs₁ - Xs₂ + n₃ Xs₃ + a₂ r₂
Q₂ = n₁ Ys₁ - Ys₂ + n₃ Ys₃ + b₂ r₂
Q₃ = n₁ Zs₁ - Zs₂ + n₃ Zs₃ + c₂ r₂
Q₄ = a₁ - a₃ c₁ / c₃
Q₅ = b₁ - b₃ a₁ / a₃
Q₆ = b₃ - b₁ a₃ / a₁
Q₇ = a₃ - a₁ c₃ / c₁

r₁ = { [ Q₁ - a₃ Q₃ / c₃ ] / [ n₁ Q₄ ] + [ Q₂ - b₃ Q₁ / a₃ ] / [ n₁ Q₅ ] } / 2
r₃ = { [ Q₂ - b₁ Q₁ / a₁ ] / [ n₃ Q₆ ] + [ Q₁ - a₁ Q₃ / c₁ ] / [ n₃ Q₇ ] } / 2

r₁ = { Rs₁² + q₁ r₁ + r₁² }^1/2
r₃ = { Rs₃² + q₃ r₃ + r₃² }^1/2

Take the values for r_i at the end of each approximation and use them as input for the next approximation. Continue repeating this procedure until the values for r_i converge.

Forming the state vector at adjusted time 2.

When we are finished converging the heliocentric distance of the asteroid at the three times of observation (adjusted for planetary abberation), we are ready to form the state vector (position and velocity) of the asteroid at adjusted time 2.

For i=1 to 3

x_i = a_i r_i - Xs_i
y_i = b_i r_i - Ys_i
z_i = c_i r_i - Zs_i

The obliquity of the Earth was (at the time of this writing)

q = 23.439281 degrees = 0.40909263 radians

X_i = x_i
Y_i = z_i sin q + y_i cos q
Z_i = z_i cos q - y_i sin q

The velocity at adjusted time 2.

b₁ = [ t^o₂ - t^o₃ ] / { [ t^o₁ - t^o₂ ] [ t^o₁ - t^o₃ ] }

b₂ = [ 2 t^o₂ - t^o₁ - t^o₃ ] / { [ t^o₂ - t^o₁ ] [ t^o₂ - t^o₃ ] }

b₃ = [ t^o₂ - t^o₁ ] / { [ t^o₃ - t^o₁ ] [ t^o₃ - t^o₂ ] }

Here's a conversion factor,

Vcf = 1731456.8367

When multiplied, it converts a velocity from AU/day to meters/second.

Vx₂ = Vcf { b₁ X₁ + b₂ X₂ + b₃ X₃ }

Vy₂ = Vcf { b₁ Y₁ + b₂ Y₂ + b₃ Y₃ }

Vz₂ = Vcf { b₁ Z₁ + b₂ Z₂ + b₃ Z₃ }

The velocity has been calculated from the derivative of a Lagrange interpolating polynomial (degree two) curvefit to the three points R_i, which are the heliocentric position vectors of the asteroid in ecliptic coordinates at the adjusted times of observation.

The state vector for the asteroid whose orbit you are trying to determine is

[ X₂, Y₂, Z₂, Vx₂, Vy₂, Vz₂ ]

The orbital elements may be calculated from this state vector with the procedure found in Lecture 2.

Worked Example.

This example follows Dubyago's demonstration problem found at the end of Chapter 5 of The Determination of Orbits, in which, however, he makes a number of typographical errors in the presentation of his figures. My own work presents the numbers that I obtained by working the same problem myself. The asteroid for which we seek orbital elements is 1590 Tsiolkovskaja, also known as 1933 NA.

Initial data.

t₁ = JD 2427255.460417 = 1 July 1933, 23h 3m 0s UT
Xs₁ = -0.169709
Ys₁ = +0.919710
Zs₁ = +0.398865
a₁ = 19.46730000 hours
d₁ = -13.86869444 degrees

t₂ = JD 2427283.391181 = 29 July 1933, 21h 23m 18s UT
Xs₂ = -0.600429
Ys₂ = +0.751016
Zs₂ = +0.325697
a₂ = 19.06218056 hours
d₂ = -14.11902778 degrees

t₃ = JD 2427312.342083 = 27 August 1933, 20h 12m 36s UT
Xs₃ = -0.908371
Ys₃ = +0.405220
Zs₃ = +0.175716
a₃ = 18.98696667 hours
d₃ = -15.24394444 degrees

Unit vectors from observer to planet.

a₁ = +0.363835156
b₁ = -0.900093900
c₁ = -0.239697622

a₂ = +0.266215600
b₂ = -0.932536300
c₂ = -0.243937090

a₃ = +0.246531192
b₃ = -0.932786462
c₃ = -0.262929245

First approximation.

t₁ = 0.498016281
t₂ = 0.978484047
t₃ = 0.480467766

n^o₁ = 0.508967195
n^o₃ = 0.491032805

n₁ = 0.060177805
n₃ = 0.059462580

d₁ = 0.015443444
d₂ = 0.018648221
d₃ = 0.017700437

D = 0.001964676

Rs₁ = 1.016740339
Rs₂ = 1.015193850
Rs₃ = 1.010058035

q₁ = 1.970356907
q₂ = 1.879285653
q₃ = 1.296252781

K₀ = 1.067106795
L₀ = 1.008749516

f{r₂} = r₂⁸ - 4.174733956 r₂⁶ + 4.048615418 r₂³ - 1.017575585

f'{r₂} = 8 r₂⁷ - 25.048403737 r₂⁵ + 12.145846255 r₂²

r₂ = 1.898670742

r₂ = 0.919728216

n₁ = 0.517759189
n₃ = 0.499720304

Q₁ = +0.303475173
Q₂ = -0.930010982
Q₃ = -0.255727953
Q₄ = +0.139086706
Q₅ = +0.476529097
Q₆ = -0.322891519
Q₇ = -0.152567065

r₁ = 0.884503909
r₃ = 1.110851828

r₁ = 1.886503769
r₃ = 1.922018156

Correction for planetary abberation.

t^o₁ = JD 2427255.455309
t^o₂ = JD 2427283.385869
t^o₃ = JD 2427312.335667

t₁ = 0.497997293
t₂ = 0.978461559
t₃ = 0.480464267

n^o₁ = 0.508959486
n^o₃ = 0.491040514

Second approximation.

x₁ = +0.491522618
y₁ = -1.715846574
z₁ = -0.610878483

x₂ = +0.845274999
y₂ = -1.608695948
z₂ = -0.550052825

x₃ = +1.182230625
y₃ = -1.441407547
z₃ = -0.467791433

K₁ = 3.801232986
K₂ = 3.732555561
K₃ = 3.766593197

h₁ = 0.005401695
h₂ = 0.021826229
h₃ = 0.005168609

y₁ = 1.005962773
y₂ = 1.023636797
y₃ = 1.005707071

n₁ = 0.061204833
n₃ = 0.059919537

n₁ = 0.517901529
n₃ = 0.499794774

K₀ = 1.067097940
L₀ = 1.020939404

f{r₂} = r₂⁸ - 4.174698414 r₂⁶ + 4.097521445 r₂³ - 1.042317267

f'{r₂} = 8 r₂⁷ - 25.048190484 r₂⁵ + 12.292564335 r₂²

r₂ = 1.896390493

p₂ = 0.917399712

[Skipping the Q's.]

r₁ = 0.882742391
r₃ = 1.107232428

r₁ = 1.884757972
r₃ = 1.918706335

The successive approximations for the heliocentric distance.

The three heliocentric positions in celestial coordinates.

x₁ = +0.490688082
y₁ = -1.713782011
z₁ = -0.610328684

x₂ = +0.844612308
y₂ = -1.606374583
z₂ = -0.549445592

x₃ = +1.181313679
y₃ = -1.437938149
z₃ = -0.466813496

The three heliocentric positions in ecliptic coordinates.

X₁ = +0.490688082
Y₁ = -1.815139083
Z₁ = +0.121737398

X₂ = +0.844612308
Y₂ = -1.692376793
Z₂ = +0.134872344

X₃ = +1.181313679
Y₃ = -1.504970227
Z₃ = +0.143685676

The adjusted times of observation.

t^o₁ = 2427255.455309
t^o₂ = 2427283.385869
t^o₃ = 2427312.335667

The velocity at adjusted time 2.

b₁ = -0.01822231549
b₂ = +0.00126052146
b₃ = +0.01696179403

Vx₂ = +21055.170 m/sec
Vy₂ = +9377.163 m/sec
Vz₂ = +673.258 m/sec

The state vector is now completed. I'll finish on this page with the derivation of the orbital elements, even though the procedure is in Lecture 2.

The orbital elements.

t^o₂ = JD 2427283.386

X₂ = +0.844612308
Y₂ = -1.692376793
Z₂ = +0.134872344

Vx₂ = +21055.170 m/sec
Vy₂ = +9377.163 m/sec
Vz₂ = +673.258 m/sec

r₂ = 1.896233031 AU = 2.836724238E+11 meters
V₂ = 23058.722 m/sec

a = 3.285211608E+11 meters
The semimajor axis: a = 2.1960283 AU

{Note: The book value of the semimajor axis, a, is 2.2300732 AU.}

Angular momentum per unit mass.

hx = -3.596521613E+14 m² sec^-1
hy = +3.397544310E+14 m² sec^-1
hz = +6.515488154E+15 m² sec^-1

h = 6.534245835E+15 m² sec^-1

The eccentricity: e = 0.14387321

{Note: The book value of the eccentricity, e, is 0.15728660.}

The inclination: i = 4.34244 degrees

{Note: The book value of the inclination, i, is 4.34657 degrees.}

Longitude of the ascending node: W = 226.62959 degrees

{Note: The book value of the longitude of the ascending node, W, is 226.70986 degrees.}

True anomaly: f = 21.20895 degrees

Argument of the perihelion: w = 48.73692 degrees

{Note: The book value of the argument of perihelion, w, is 52.63858 degrees.}

ex = 0.1365170037
ey = 0.0454159545

Eccentric anomaly: u = 18.40105 degrees

Mean anomaly: M = 15.79891 degrees

Period: P = 1188.6536 days

Time of perihelion passage: T = JD 2427231.2208

{According to Dubyago, asteroid 1590 Tsiolkovskaja (1933 NA) was at perihelion at JD 2427236.0497, so our figure was about five days early.}

Note: Don't use this set of orbital elements in hopes of locating this asteroid. The data used is over 70 years old, and the orbit calculated isn't good enough to track it for that long. This procedure calculates a preliminary set of orbital elements which needs improving by a more advanced method using a greater number of observations.

Racist Celestial Mechanics
Jerry's Aryan Battle Page