To help me in understanding this procedure, I read the fifth chapter of The Determination of Orbits by A. D. Dubyago, as translated from Russian into English by the RAND Corporation. Much of the nomenclature used in this exposition is patterned after Dubyago.

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 unit vector in the direction of the planet: RA(i), DEC(i)

The position vectors are in local celestial coordinates.

Handy Constants.

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

Unit vector in the direction of the planet with rectangular components.

a(i) = cos RA(i) cos DEC(i)
b(i) = sin RA(i) cos DEC(i)
c(i) = sin DEC(i)

First Approximation.

tau(1) = k {t(3)-t(2)}
tau(2) = k {t(3)-t(1)}
tau(3) = k {t(2)-t(1)}

nc(1) = tau(1)/tau(2)
nc(3) = tau(3)/tau(2)

nu(1) = tau(1) tau(3) {1 + nc(1)} / 6
nu(3) = tau(1) tau(3) {1 + nc(3)} / 6

D = a(2) { b(3) c(1) - b(1) c(3) } + b(2) { c(3) a(1) - c(1) a(3) } + c(2) { a(3) b(1) - a(1) b(3) }

For i=1 to 3
d(i) = Xs(i) { b(3) c(1) - b(1) c(3) } + Ys(i) { c(3) a(1) - c(1) a(3) } + Zs(i) { a(3) b(1) - a(1) b(3) }

For i=1 to 3
Re(i) = { Xs(i)^2 + Ys(i)^2 + Zs(i)^2 }^0.5

For i=1 to 3
q(i) = -2 { a(i) Xs(i) + b(i) Ys(i) + c(i) Zs(i) }

K0 = { d(2) - d(1) nc(1) - d(3) nc(3) } / D

L0 = { d(1) nu(1) + d(3) nu(3) } / D

Let r(2) be the distance from the sun to the planet at time 2
Let p(2) be the distance from Earth to the planet at time 2.

We must solve these two equations simultaneously for r(2) and p(2):

p(2) = K0 - L0 / r(2)^3

r(2) = { Re(2)^2 + q(2) p(2) + p(2)^2 }^0.5

This leads to an 8th degree polynomial in r(2):

f{r(2)} = r(2)^8 - { Re(2)^2 + q(2) K0 + K0^2 } r(2)^6 + L0 {q(2) + 2 K0} r(2)^3 - L0^2
f{r(2)} = 0

We will use Newton's method to get the answer.

df/dr(2) = 8 r(2)^7 - 6 { Re(2)^2 + q(2) K0 + K0^2 } r(2)^5 + 3 L0 {q(2) + 2 K0} r(2)^2

As an initial guess at r(2)...

r2(j=0) = 1.0

Then...

Repeat over j

r2(j+1) = r2(j) - f{r2(j)} / df/dr2(j)

Until abs{ r2(j+1) - r2(j) } < 1.0E-12

Assign to r(2) the converged value from the loop:

r(2) = r2(j+1)

p(2) = K0 - L0 / r(2)^3

n(1) = nc(1) + nu(1) / r(2)^3
n(3) = nc(3) + nu(3) / r(2)^3

Dubyago left a system of equations unsolved in his chapter, no doubt figuring that his readers wouldn't mind the exercise. I figured that you probably would mind, so I solved it for p(1) and p(3).

Q1 = n(1) Xs(1) - Xs(2) + n(3) Xs(3) + a(2) p(2)
Q2 = n(1) Ys(1) - Ys(2) + n(3) Ys(3) + b(2) p(2)
Q3 = n(1) Zs(1) - Zs(2) + n(3) Zs(3) + c(2) p(2)
Q4 = a(1) - a(3) c(1) / c(3)
Q5 = b(1) - b(3) a(1) / a(3)
Q6 = b(3) - b(1) a(3) / a(1)
Q7 = a(3) - a(1) c(3) / c(1)

p(1) = { [ Q1 - a(3) Q3 / c(3) ] / [ n(1) Q4 ] + [ Q2 - b(3) Q1 / a(3) ] / [n(1) Q5] } / 2
p(3) = { [ Q2 - b(1) Q1 / a(1) ] / [ n(3) Q6 ] + [ Q1 - a(1) Q3 / c(1) ] / [n(3) Q7] } / 2

r(1) = { Re(1)^2 + q(1) p(1) + p(1)^2 }^0.5
r(3) = { Re(3)^2 + q(3) p(3) + p(3)^2 }^0.5

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 planet 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
tc(i) = t(i) - A p(i)

Where A is the time (in days) required for light to travel 1 astronomical unit.

tau(1) = k {tc(3)-tc(2)}
tau(2) = k {tc(3)-tc(1)}
tau(3) = k {tc(2)-tc(1)}

nc(1) = tau(1)/tau(2)
nc(3) = tau(3)/tau(2)

Recursive Procedure for Successive Approximations.

For i=1 to 3
x(i) = a(i) p(i) - Xs(i)

y(i) = b(i) p(i) - Ys(i)
z(i) = c(i) p(i) - Zs(i)

K(1) = { 2 [ r(2) r(3) + x(2) x(3) + y(2) y(3) + z(2) z(3) ] }^0.5
K(2) = { 2 [ r(1) r(3) + x(1) x(3) + y(1) y(3) + z(1) z(3) ] }^0.5
K(3) = { 2 [ r(1) r(2) + x(1) x(2) + y(1) y(2) + z(1) z(2) ] }^0.5

h(1) = tau(1)^2 / { K(1)^2 [ K(1)/3 + ( r(2)+r(3) )/2 ] }
h(2) = tau(2)^2 / { K(2)^2 [ K(2)/3 + ( r(1)+r(3) )/2 ] }
h(3) = tau(3)^2 / { K(3)^2 [ K(3)/3 + ( r(1)+r(2) )/2 ] }

For i=1 to 3
Begin
Queasy1 = (11/9) h(i)
Queasy2 = Queasy1
Repeat
Queasy3 = Queasy2
Queasy2 = Queasy1 / ( 1 + Queasy2)
Until abs( Queasy2 - Queasy3 ) / Queasy3 < 1E-12;
Yippy(i) = 1 + (10/11) Queasy2
End

n(1) = nc(1) yippy(2) / yippy(1)
n(3) = nc(3) yippy(2) / yippy(3)

nu(1) = nc(1) r(2)^3 { yippy(2) / yippy(1) - 1 }
nu(3) = nc(3) r(2)^3 { yippy(2) / yippy(3) - 1 }

K0 = { d(2) - d(1) nc(1) - d(3) nc(3) } / D

L0 = { d(1) nu(1) + d(3) nu(3) } / D

p(2) = K0 - L0 / r(2)^3

r(2) = { Re(2)^2 + q(2) p(2) + p(2)^2 }^0.5

f{r(2)} = r(2)^8 - { Re(2)^2 + q(2) K0 + K0^2 } r(2)^6 + L0 {q(2) + 2 K0} r(2)^3 - L0^2
f{r(2)} = 0

df/dr(2) = 8 r(2)^7 - 6 { Re(2)^2 + q(2) K0 + K0^2 } r(2)^5 + 3 L0 {q(2) + 2 K0} r(2)^2

The customary initial guess at r(2)...

r2(j=0) = 1.0

Repeat over j

r2(j+1) = r2(j) - f{r2(j)} / df/dr2(j)

Until abs{ r2(j+1) - r2(j) } < 1.0E-12

Assign to r(2) the converged value from the loop:

r(2) = r2(j+1)

p(2) = K0 - L0 / r(2)^3

Q1 = n(1) Xs(1) - Xs(2) + n(3) Xs(3) + a(2) p(2)
Q2 = n(1) Ys(1) - Ys(2) + n(3) Ys(3) + b(2) p(2)
Q3 = n(1) Zs(1) - Zs(2) + n(3) Zs(3) + c(2) p(2)
Q4 = a(1) - a(3) c(1) / c(3)
Q5 = b(1) - b(3) a(1) / a(3)
Q6 = b(3) - b(1) a(3) / a(1)
Q7 = a(3) - a(1) c(3) / c(1)

p(1) = { [ Q1 - a(3) Q3 / c(3) ] / [ n(1) Q4 ] + [ Q2 - b(3) Q1 / a(3) ] / [n(1) Q5] } / 2
p(3) = { [ Q2 - b(1) Q1 / a(1) ] / [ n(3) Q6 ] + [ Q1 - a(1) Q3 / c(1) ] / [n(3) Q7] } / 2

r(1) = { Re(1)^2 + q(1) p(1) + p(1)^2 }^0.5
r(3) = { Re(3)^2 + q(3) p(3) + p(3)^2 }^0.5

Repeat this whole section (Recursive Procedure for Successive Approximations) until your values for r(i) converge.

Higher Approximations.

You see how it works. You just keep taking p(i) from the end of the latest approximation and plugging it back into the top of the next approximation and repeating the math with the improved numbers. You keep repeating the recursive section until all the values for r(i) have converged.

Jerry Abbott

Example (after Dubyago with improved accuracy).

Initial data.

The "planet" we will use is the asteroid 1590 Tsiolkovskaja, also known as 1933 NA.

t(1) = JD 2427255.460417 = 1 July 1933, 23h 3m 0s UT
Xs(1) = -0.169709
Ys(1) = +0.919710
Zs(1) = +0.398865
RA(1) = 19.46730000 hours
DEC(1) = -13.86869444 degrees

t(2) = JD 2427283.391181 = 29 July 1933, 21h 23m 18s UT
Xs(2) = -0.600429
Ys(2) = +0.751016
Zs(2) = +0.325697
RA(2) = 19.06218056 hours
DEC(2) = -14.11902778 degrees

t(3) = JD 2427312.342083 = 27 August 1933, 20h 12m 36s UT
Xs(3) = -0.908371
Ys(3) = +0.405220
Zs(3) = +0.175716
RA(3) = 18.98696667 hours
DEC(3) = -15.24394444 degrees

Unit vectors from observer to planet.

a(1) = +0.363835156
b(1) = -0.900093900
c(1) = -0.239697622

a(2) = +0.266215600
b(2) = -0.932536300
c(2) = -0.243937090

a(3) = +0.246531192
b(3) = -0.932786462
c(3) = -0.262929245

First Approximation.

tau(1) = 0.498016281
tau(2) = 0.978484047
tau(3) = 0.480467766

nc(1) = 0.508967195
nc(3) = 0.491032805

nu(1) = 0.060177805
nu(3) = 0.059462580

d(1) = 0.015443444
d(2) = 0.018648221
d(3) = 0.017700437

D = 0.001964676

Re(1) = 1.016740339
Re(2) = 1.015193850
Re(3) = 1.010058035

q(1) = 1.970356907
q(2) = 1.879285653
q(3) = 1.296252781

K0 = 1.067106795
L0 = 1.008749516

f(r2) = r(2)^8 - 4.174733956 r(2)^6 + 4.048615418 r(2)^3 - 1.017575585

f'(r2) = 8 r(2)^7 - 25.048403737 r(2)^5 + 12.145846255 r(2)^2

r(2) = 1.898670742

p(2) = 0.919728216

n(1) = 0.517759189
n(3) = 0.499720304

Q1 = +0.303475173
Q2 = -0.930010982
Q3 = -0.255727953
Q4 = +0.139086706
Q5 = +0.476529097
Q6 = -0.322891519
Q7 = -0.152567065

p(1) = 0.884503909
p(3) = 1.110851828

r(1) = 1.886503769
r(3) = 1.922018156

Correction for planetary abberation.

tc(1) = JD 2427255.455309
tc(2) = JD 2427283.385869
tc(3) = JD 2427312.335667

tau(1) = 0.497997293
tau(2) = 0.978461559
tau(3) = 0.480464267

nc(1) = 0.508959486
nc(3) = 0.491040514

Second Approximation.

x(1) = +0.491522618
y(1) = -1.715846574
z(1) = -0.610878483

x(2) = +0.845274999
y(2) = -1.608695948
z(2) = -0.550052825

x(3) = +1.182230625
y(3) = -1.441407547
z(3) = -0.467791433

K(1) = 3.801232986
K(2) = 3.732555561
K(3) = 3.766593197

h(1) = 0.005401695
h(2) = 0.021826229
h(3) = 0.005168609

Yippy(1) = 1.005962773
Yippy(2) = 1.023636797
Yippy(3) = 1.005707071

nu(1) = 0.061204833
nu(3) = 0.059919537

n(1) = 0.517901529
n(3) = 0.499794774

K0 = 1.067097940
L0 = 1.020939404

f(r2) = r(2)^8 - 4.174698414 r(2)^6 + 4.097521445 r(2)^3 - 1.042317267

f'(r2) = 8 r(2)^7 - 25.048190484 r(2)^5 + 12.292564335 r(2)^2

r(2) = 1.896390493

p(2) = 0.917399712

[Skipping the Q's.]

p(1) = 0.882742391
p(3) = 1.107232428

r(1) = 1.884757972
r(3) = 1.918706335

Jerry Abbott

The 1st approximation (derived explicitly in previous post) is
p(1,2,3) = 0.884503909, 0.919728216, 1.110851828
r(1,2,3) = 1.886503769, 1.898670742, 1.922018156

The 2nd approximation (derived explicitly in previous post) is
p(1,2,3) = 0.882742391, 0.917399712, 1.107232428
r(1,2,3) = 1.884757972, 1.896390493, 1.918706335

The 3rd approximation is
p(1,2,3) = 0.882277763, 0.917293056, 1.107206810
r(1,2,3) = 1.884297496, 1.896286050, 1.918682897

The 4th approximation is
p(1,2,3) = 0.882223313, 0.917244422, 1.107137725
r(1,2,3) = 1.884243532, 1.896238425, 1.918619694

The 5th approximation is
p(1,2,3) = 0.882212065, 0.917240187, 1.107134081
r(1,2,3) = 1.884232385, 1.896234278, 1.918616361

The 6th approximation is
p(1,2,3) = 0.882210524, 0.917239076, 1.107132612
r(1,2,3) = 1.884230857, 1.896233190, 1.918615016

The 7th approximation is
p(1,2,3) = 0.882210241, 0.917238946, 1.107132476
r(1,2,3) = 1.884230577, 1.896233062, 1.918614892

The 8th approximation is
p(1,2,3) = 0.882210199, 0.917238919, 1.107132442
r(1,2,3) = 1.884230536, 1.896233036, 1.918614861

The 9th approximation is
p(1,2,3) = 0.882210192, 0.917238915, 1.107132438
r(1,2,3) = 1.884230529, 1.896233032, 1.918614857

The 10th approximation is
p(1,2,3) = 0.882210191, 0.917238914, 1.107132437
r(1,2,3) = 1.884230527, 1.896233032, 1.918614856

The 11th approximation is
p(1,2,3) = 0.882210191, 0.917238914, 1.107132437
r(1,2,3) = 1.884230527, 1.896233032, 1.918614856

And we're done approximating distances.

Jerry Abbott

The positions at times 1,2,3.

For i=1 to 3
x(i) = a(i) p(i) - Xs(i)
y(i) = b(i) p(i) - Ys(i)
z(i) = c(i) p(i) - Zs(i)

Rotate the position vectors into ecliptic coordinates.

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

c1 = [ t(2) - t(3) ] / { [ t(1) - t(2) ] [ t(1) - t(3) ] }
c2 = [ 2 t(2) - t(1) - t(3) ] / { [ t(2) - t(1) ] [ t(2) - t(3) ] }
c3 = [ t(2) - t(1) ] / { [ t(3) - t(1) ] [ t(3) - t(2) ] }

Here's a conversion factor,

Velcf = 1731456.8367

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

Vx = Velcf { c1 X(1) + c2 X(2) + c3 X(3) }
Vy = Velcf { c1 Y(1) + c2 Y(2) + c3 Y(3) }
Vz = Velcf { c1 Z(1) + c2 Z(2) + c3 Z(3) }

The velocity has been calculated from the derivative of a Lagrange interpolating polynomial (degree two) curvefit to the three points r(i) we calculated in previous posts, except the r(i) vectors were rotated first from heliocentric celestial to heliocentric ecliptic coordinates.

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

[ X(2), Y(2), Z(2), Vx, Vy, Vz ]

Plugging in the numbers...

The values for p(i) in the last approximation are

p(1) = 0.882210191
p(2) = 0.917238914
p(3) = 1.107132437

Here are the three heliocentric position vectors of the planet at the three selected times of observation, in celestial coordinates. (See the initial data section of the top post for a(i), b(i), c(i), Xs(i), Ys(i), Zs(i).)

x(1) = +0.490688082
y(1) = -1.713782011
z(1) = -0.610328684

x(2) = +0.844612308
y(2) = -1.606374583
z(2) = -0.549445592

x(3) = +1.181313679
y(3) = -1.437938149
z(3) = -0.466813496

Here are the three heliocentric position vectors of the planet at the three selected times of observation, in ecliptic coordinates.

X(1) = +0.490688082
Y(1) = -1.815139083
Z(1) = +0.121737398

X(2) = +0.844612308
Y(2) = -1.692376793
Z(2) = +0.134872344

X(3) = +1.181313679
Y(3) = -1.504970227
Z(3) = +0.143685676

Here are the values for r(i) for the last approximation, so you can check them with the magnitude of the vectors above.

r(1) = 1.884230527
r(2) = 1.896233032
r(3) = 1.918614856

The times of observation, corrected for planetary abberation, are

tc(1) = 2427255.455309
tc(2) = 2427283.385869
tc(3) = 2427312.335667

c1 = -0.01822231549
c2 = +0.00126052146
c3 = +0.01696179403

The sun-relative velocity at t(2), referred to ecliptic coordinates, is

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

Jerry Abbott

We have a heliocentric vector of position, a sun-relative velocity, and a time, from which we want to determine the elements of the orbit.

The time is: tc(2)
The heliocentric position vector is: X(2), Y(2), Z(2)
The sun-relative velocity is: Vx, Vy, Vz

Important note. In some of the equations, it will be necessary to apply distance in astronomical units, while in other equations, they must be applied in meters. I leave it to you to figure out which is required to keep the units consistent. And likewise for angles in either radians or degrees.

r = [ x(2)^2 + y(2)^2 + z(2)^2 ]^0.5

V = [ Vx^2 + Vy^2 + Vz^2 ]^0.5

GMsun = 1.32712440018E+20 m^3 sec^-2

1 AU = 1.49597870691E+11 meters

We solve for the semimajor axis of the orbit by putting MKS values for (r) and (V) into the Vis Viva equation:

a = [ 2/r - V^2 / GMsun ]^-1

We solve for the components of the orbit's angular momentum vector (per unit mass) with a cross product.

hx = Y(2) Vz - Z(2) Vy
hy = Z(2) Vx - X(2) Vz
hz = X(2) Vy - Y(2) Vx

h = {hx^2 + hy^2 + hz^2}^0.5

We solve for the orbital eccentricity:

e = { 1 - h^2 / (a GMsun) }^0.5

We solve for the inclination of the orbit to the ecliptic:

i = pi/2 - arcsin(hz/h)

We solve for the longitude of the ascending node:

L = arctan2(hX , -hY)

The arctan2 function is the arctangent adjusted for quadrant by the arguments. It is defined more fully in my Hillbilly Tutorial on Transfer Orbits. The 1st argument is proportional to the sine of the angle, while the 2nd argument is proportional to the cosine of the angle.

We solve for the true anomaly:

fx = h^2 / { r(2) GMsun } - 1
fy = h { X(2) Vx + Y(2) Vy + Z(2) Vz } / { r(2) GMsun }

f = arctan2( fy , fx )

We solve for the argument of the perihelion:

cos(f+w) = { X(2) cos L + Y(2) sin L } / r(2)
sin(f+w) = Z(2) / { r(2) sin i }

f + w = arctan2{ sin(f+w) , cos(f+w) }

w = (f + w) - f

If w<0 then correct it to the interval [0, 2 pi radians).

We solve for the eccentric anomaly.

ex = 1 - r(2) / a

ey = { X(2) Vx + Y(2) Vy + Z(2) Vz } / { a GMsun }^0.5

u = arctan2(ey,ex)

We solve for the mean anomaly.

M = u - e sin u

The eccentric anomaly, u, must be entered to the equation immediately above in radians! M will result in radians. You can convert it to degrees for display purposes.

We solve for the period.

P = 365.256898326 a^1.5

We solve for the time of perihelion passage.

T = tc(2) - P M / (2 pi radians)

Plugging in the numbers...

t(2) = JD 2427283.386

X(2) = +0.844612308
Y(2) = -1.692376793
Z(2) = +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
a = 2.1960283 AU

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

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

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

e = 0.14387321

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

i = 4.34244 degrees

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

L = 226.62959 degrees

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

f = 21.20895 degrees

w = 48.73692 degrees

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

ex = 0.1365170037

ey = 0.0454159545

u = 18.40105 degrees

M = 15.79891 degrees

P = 1188.6536 days

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.

Jerry Abbott

Calculated Orbital Elements of 1590 Tsiolkovskaja (1933 NA)
(Epoch 1933.)

a = 2.1960283 AU
e = 0.14387321
i = 4.34244 degrees
L = 226.62959 degrees
w = 48.73692 degrees
T = JD 2427231.2208

There. Now you, too, can figure preliminary orbits with nothing more than a telescope, a timepiece and a calculator. However, remember that the accuracy of the elements of an orbit determined this way is very sensitive to the precision with which the observations are made. You would have to use several triplet sets, involving different observations RA(i) and DEC(i), before you could have confidence in your preliminary orbit.

The improvement of a preliminary orbit is a more complicated affair, which requires a large number of observations as input for a reduction with a least-squares Taylor series kind of method in which partial derivatives of RA and DEC with respect to each orbital element go flying every which way. But if finding orbits is what you wanted to do, the procedure given above will get you started.

Jerry Abbott