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
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
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
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.
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
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
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
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
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: 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....
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?
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
If your calculator or computer assumes that the trig function arguments are radians, you do need to convert them. Jerry Abbott
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
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?
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