Orbit Determination for Dummies Another hillbilly tutorial.

[Jenab]

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

 

[Jenab]

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

 

[Jenab]

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

 

[enigma]

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?

 

[BobG]

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?

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.

 

[Jenab]

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?

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

 

[Jenab]

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

 

[Jenab]

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

 

[Jenab]

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

 

[Bumbleflea]

Not for Dummies!

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

Me->
All of you:->

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:

For i=1 to 3

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?

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)

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:

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)
...

Tau? ... TAU?!? Wtf is Tau?? Dummies want to know!
NC? NC TAU?? Sounds like one of those rap-guys...

 

[dagmare]

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?

 

[Jenab2]

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?

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.

Tau? ... TAU?!? Wtf is Tau?? Dummies want to know!
NC? NC TAU?? Sounds like one of those rap-guys...

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

 

[Jenab2]

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?

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

Jerry Abbott

 

[Bumbleflea]

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

Ah, thanks for your help, Jerry!

 

[xonics]

are these all relating to the aspect of vectors?

 

[Jenab6]

are these all relating to the aspect of vectors?

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

 

[fructivore]

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?

 

[varun atri]

hey guys i am a student of mechanical engineering. I am very much interested in astronomy how can i enter in the astronomy deptt. after my graduation in mechanical engineering

 

[Drakkith]

hey guys i am a student of mechanical engineering. I am very much interested in astronomy how can i enter in the astronomy deptt. after my graduation in mechanical engineering

Please post in the Academic Guidance section of the forum.

 

[varun atri]

hello drakkith where i could find the academic guidance section..

 

[Drakkith]

hello drakkith where i could find the academic guidance section..

Main page, 3rd forum down from the top.