Orbit Determination for Dummies! Another hillbilly tutorial.

  1. Jenab

    Jenab 115
    Science Advisor

    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.

    MORE ORBIT DETERMINATION IN POSTS TO FOLLOW.


    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.

    Jerry Abbott
     
    Last edited: Jul 31, 2004
  2. jcsd
  3. Jenab

    Jenab 115
    Science Advisor

    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.

    Jerry Abbott
     
    Last edited: Jul 26, 2004
  4. Jenab

    Jenab 115
    Science Advisor

    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

    Jerry Abbott
     
    Last edited: Jul 31, 2004
  5. enigma

    enigma 1,815
    Staff Emeritus
    Science Advisor
    Gold Member

    Your speed of light is wrong. Maybe it's in AU/second?

    Light goes 1 AU in ~8 minutes, IIRC.

    For those of us who don't have the book, what are you trying to find?
     
  6. BobG

    BobG 2,357
    Science Advisor
    Homework Helper

    If he's using the normal definition of AU, then speed of light in AU/sec is .002004. At least in the same area code.

    Mean motion is pretty close.

    I'd suggest using the International Earth Rotation Service's constants:
    http://hpiers.obspm.fr/eop-pc/models/constants.html

    Their job is measuring earth rotation parameters, so their constants are a little more accurate and current than some of the older books.
     
  7. Jenab

    Jenab 115
    Science Advisor

    Sorry, I got in a hurry and wrote down the wrong thing.

    A is the number of days required by light to travel 1 AU. Instead of the speed of light, it's sort of the reciprocal of the speed of light, in days/AU.

    I'm getting to the orbital elements of this asteroid, using the Method of Gauss on angle-only observation data. Dubyago's a clear read on this subject, but he made a lot of typos in his worked example problem. So I don't rely on his answers.

    Jerry Abbott
     
    Last edited: Jul 28, 2004
  8. Jenab

    Jenab 115
    Science Advisor

    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.

    Jerry Abbott
     
    Last edited: Aug 4, 2004
  9. Jenab

    Jenab 115
    Science Advisor

    We have reached the end of this dissertation. I must now go feed my goats and start fixing a kettle of stew for my supper.

    Jerry Abbott
     
  10. Jenab

    Jenab 115
    Science Advisor

    I have a FreePascal program with the above "method of Gauss" procedure online at

    http://www.jabpage.org/features/gauss.zip

    It's about 26 kB for the zip file. It contains the source code, the executable and the input data file. Output is to a DOSPROMPT window.

    This program is for WIN32 type OS, such as WIN98 and such.

    Jerry Abbott
     
  11. Not for Dummies!

    Hello, I'm Bumbleflea & this is my 1st post. I've come here as a layman seeking knowledge:

    Me->:confused:
    All of you:->:rolleyes::rolleyes::rolleyes::rolleyes:

    Please bear with me. First, what I'm interested in most is understanding the basic physics of our solar system, which led me to this sticky post.

    My first impression is that this isn't written for dummies; it's written for students of astrophysics. Point in case:
    What the heck is 'i'?? A point in time? Are we simply talking about iterations of a loop, as in BASIC programming? If so, then why is the scope limited to 1<i<3?
    Okay, so we're talking about one/two/three points in time. Got it.
    I know what an 'abc' unit vector is & constants are fine with me, but then we get to:

    Tau? ... TAU?!? Wtf is Tau?? Dummies wanna know!
    NC? NC TAU?? Sounds like one of those rap-guys....:bugeye:
     
  12. 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

    am i right?
     
    Last edited: Oct 23, 2007
  13. i is an index variable. Part of the input for this algorithm are the geocentric right ascensions and declinations for an asteroid (or planet) at THREE DIFFERENT TIMES.

    t1, RA1, DEC1
    t2, RA2, DEC2
    t3, RA3, DEC3

    Or...

    For i=1 to 3: t(i), RA(i), DEC(i)

    If these indexed variables occur in a long expression, and if we must cycle over each time, it saves typing to use the indices.

    The taus are time differences with changed units. The index on the tau is the index missing from the two t's being subtracted (with the earlier t being subtracted from the later one, to keep all the differences positive). The NC and etcetera are geometric quantities, some of which have to do with Kepler's equal areas in equal times law.

    Jerry Abbott
     
  14. If your calculator or computer assumes that the trig function arguments are radians, you do need to convert them.

    Jerry Abbott
     
  15. Ah, thanks for your help, Jerry!
     
  16. are these all relating to the aspect of vectors?
     
  17. Certainly, there is a lot of vector arithmetic used in celestial mechanics work. We do vector subtractions and additions, we rotate from one coordinate system to another, we convert from spherical to rectangular and vice versa. We sometimes have to find a vector's direction and magnitude from two separate calculations and then combine them by multiplying the scalar by the unit vector. Being handy with vector arithmetic (and with matrix arithmetic) is necessary for seeing in your mind what the math is all about.

    Jerry Abbott
     
  18. Hello,

    I am developing a simple orbital mechanics simulator as a software project and am approaching the problem from a slightly different perspective from the one described here. Specifically I have information on satellite position and velocity in spherical polar coordinates (radius, latitude, longitude; radial speed, latitudinal speed, longitudinal speed) rather than Cartesian coordinates.

    My intuition is that this should make orbital calculations simpler and so I shouldn't have to back-convert to Cartesian before working the problem out. But my weak math brain can't wrap itself around how to adapt some of the hillbilly equations to a polar system.

    Any suggestions?
     
  19. hey guys i am a student of mechanical engineering. I am very much intrested in astronomy how can i enter in the astronomy deptt. after my graduation in mechanical engineering
     
  20. Drakkith

    Staff: Mentor

    Please post in the Academic Guidance section of the forum.
     
  21. hello drakkith where i could find the academic guidance section..
     
Know someone interested in this topic? Share a link to this question via email, Google+, Twitter, or Facebook

Have something to add?