Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Sep 1, 2018 #1
    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.
  2. jcsd
  3. Sep 3, 2018 #2
    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.
  4. Sep 10, 2018 #3
    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


    ORBIT4 is designed to treat elliptical (two-body) orbits around the sun. It is written in Prime Programming Language (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

    // Data for time 2

    // Data for time 3

    // Data for time 4

    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.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted