How to get state vector in 4-observation method of Gauss

[Jenab2]

In the four-observation method of Gauss for orbit determination, the right ascension and declination of an asteroid is observed at specified times, and the heliocentric position of Earth is obtained from tables (or JPL Horizons) for those same times.

I can follow the procedure to the point where the asteroid's heliocentric and geocentric distances converge as the result of successive approximations, leaving me with position vectors for the asteroid in ecliptic coordinates, sun-to-asteroid (x,y,z) and Earth-to-asteroid (X,Y,Z), for time 1 and for time 4.

From there, however, my textbook goes into some weird process for obtaining the state vector that is to be used in finding the elements.

I managed to obtain a state vector by correcting the times of observation for planetary aberration, then using the average position (r₁+r₄)/2 and a numerical derivative for the velocity (r₄−r₁)/(t₄−t₁). That works fairly well, though it seems crude, as if I'm not using all my information to best advantage. I'm thinking that the textbook probably gives a more accurate way of obtaining a state vector for the asteroid. I just don't understand what's going on. Like I said, it looks weird.

 

[Jenab2]

I understand the problem better now. The reason for the crudeness of my earlier method for composing the position portion of the state vector is that I was assigning to it the midpoint of a secant line that cut across the orbit. In reality, this midpoint only provides the direction. The distance should be taken as the average between the magnitudes of r₁ and r₄. When I corrected the magnitude of the heliocentric position vector to the curved path of the orbit, I had a small but significant improvement in the orbital elements as a result.

It might be worthwhile to determine whether the velocity part of the state vector can be similarly adjusted, hopefully resulting in a further improvement in the preliminary orbit.

 

[Jenab2]

Okay, after some research I've answered my own question. The reference material that I used is The Determination of Orbits by A.D. Dubyago, especially chapter six.

I used the method presented by Dubyago, which is a form of the method of Gauss, in my program ORBIT4. Actually, I followed Dubyago only until I had the range information. After that, I went my own way with a simpler route of march toward the orbital elements, which I understood better.

There is no need to have a priori range information. The method of Gauss ferrets it out by successive approximations.

ORBIT4 is available for free download from

https://my.cloudme.com/jenab6/ORBIT4

ORBIT4 is designed to treat elliptical (two-body) orbits around the sun. It is written in Prime Programming Langauge (PPL) for the HP Prime Calculator or its emulators. The input, in this case for the asteroid 1 Ceres, is in the program (inline code), and looks like this:

// Data for time 1
L1:=
{2457204.625,
0.155228396,
−1.004732775,
0.00003295786,
HMS→(20°46′57.02″),
HMS→(−27°41′33.9″)};

// Data for time 2
L2:=
{2457214.625,
0.319493277,
−0.965116604,
0.0000311269,
HMS→(20°39′57.10″),
HMS→(−28°47′21.5″)};

// Data for time 3
L3:=
{2457224.625,
0.4747795623,
−0.8983801739,
0.00002841127,
HMS→(20°31′22.81″),
HMS→(−29°49′22.7″)};

// Data for time 4
L4:=
{2457234.625,
0.616702829,
−0.8063620175,
0.00002486325,
HMS→(20°22′06.57″),
HMS→(−30°41′57.3″)};

The first number (after the open curly bracket) is the time of observation in Julian Date. The observations should be separated in time by somewhere between 0.5% to 1.0% of the object's orbital period. The observation interval should be reasonably near the opposition of the asteroid with the Sun, but it should not span an apside of the asteroid's orbit.

The 2nd number is the X component of the Earth's position in heliocentric ecliptic coordinates. The 3rd number is the Y component. The 4th number is the Z component. These can be obtained from JPL Horizons.

The 5th number, inside the HMS operand parentheses, is the asteroid's geocentric right ascension in HH°MM'SS.SS" format. Right ascension should be accurate to 0.01 seconds of time.

The 6th number, also inside HMS operand parantheses, is the asteroid's geocentric declination in degrees, arcminutes, arcseconds format. Declination should be accurate to 0.1 arcsec.

If observational data aren't available, then test data for the RA & DEC of known asteroids or planets can be obtained from JPL Horizons.

The output for the data (shown in the PPL code above) follows:

ORBIT4 by David Sims
Method of Gauss with four observed positions to find the Keplerian orbital elements. User provides input by adjusting inline data.

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

Successive approximations
r₁ 2.90652064 r₄ 2.92071388
r₁ 2.93008666 r₄ 2.94292188
r₁ 2.93307059 r₄ 2.94572742
r₁ 2.93344165 r₄ 2.94607627
r₁ 2.93348769 r₄ 2.94611956
r₁ 2.9334934 r₄ 2.94612492
r₁ 2.93349411 r₄ 2.94612559
r₁ 2.9334942 r₄ 2.94612567
r₁ 2.93349421 r₄ 2.94612568
r₁ 2.93349421 r₄ 2.94612568

Heliocentric distances in AU at t₁ & t₄
r₁ 2.93349421
r₄ 2.94612568

Geocentric distances in AU at t₁ & t₄
ρ₁ 2.00460681
ρ₄ 1.94781669

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

Aberration corrections to time
Aρ₁ 0.011577643 days
Aρ₄ 0.011249651 days

Epoch of state vector & obliquity
t₀ 2457219.61 JD
ε₀ 0.409057547 radians

HEC state vector
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

Heliocentric distance
r 2.93980995 AU

Sun−relative speed
v 16819.9661 m/s

True anomaly 147.669798°
Ecc. anomaly 145.259666°
Mean anomaly 142.77737°

Period of orbit 1681.12408 days

Orbital elements
a 2.76694735 AU
e 0.076026341
i 10.5918141°
Ω 80.3183813°
ω 72.6265868°
T 2456552.87 JD
T+P 2458234 JD

For comparison, here is the official NASA/JPL orbit for Ceres:

a 2.768008676 AU
e 0.075773357
i 10.59221734°
Ω 80.32683297°
ω 72.66267214°
T 2456552.644 JD
P 1682.091 days

The source code can be found on my LiveJournal website.