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: a(i), b(i), c(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)
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) { c(1) b(3) - c(3) b(1) } + Ys(i) { a(1) c(3) - a(3) c(1) } + Zs(i) { b(1) a(3) - b(3) a(1) }
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
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.

Example (after Dubyago with improved accuracy).

Initial data. In practice, what you determine observationally about a planet's position is right
ascension and declination. So I'll present those as initial data and calculate a(i), b(i), c(i) from
them. Remember that Xs,Ys,Zs is the position of the sun, while RA and DEC pertain to the
planet whose orbit we're trying to determine. (By the way, the planet is 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
The Lagrange equation is
f(r2) = r(2)^8 - 4.174733956 r(2)^6 + 4.048615418 r(2)^3 - 1.017575585
The first derivative of the Lagrange equation is
f'(r2) = 8 r(2)^7 - 25.048403737 r(2)^5 + 12.145846255 r(2)^2
This actually takes a while to converge with Newton's method, so I'd consider going to Danby's
method for faster convergence. Anyway, it converges to
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
The Lagrange equation is
f(r2) = r(2)^8 - 4.174698414 r(2)^6 + 4.097521445 r(2)^3 - 1.042317267
The first derivative of the Lagrange equation is
f'(r2) = 8 r(2)^7 - 25.048190484 r(2)^5 + 12.292564335 r(2)^2
Newton's method converges to...
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
Higher Approximations. You've seen 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.

Successive approximations for the distances.

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.

Composing the state vector at time 2.

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

The orbital elements.

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.

the unit vectors are:

for i = 1 .. 3
a(i) = cos(DEC(i))*cos(RA(i))
b(i) = cos(DEC(i))*sin(RA(i))
c(i) = sin(DEC(i))
next i
but before calculating that you need to convert hours to radians and deg to radians using:
for a given value in hours:
Value_in_rad = Value_in_h * PI/12.0
for a given value in deg Value_in_rad = Value_in_red * PI/180