Selected Questions- Answers From All Experts Astronomy Forum (Celestial Mechanics)

Question: I am interested in how spacecraft get to other planets and how astronomers are able to predict future planet positions to do this. So basically can you provide some background details on what makes up celestial mechanics?  (P.S. I am a math major so don't hold back on math details! :) ) - Barbara, Santa Monica CA

Answer:

The name "celestial mechanics"  conveys exactly what the subject embodies: Mechanics applied to the dynamics of celestial objects. A typical diagram from my own Celestial Mechanics notes when I took it at Univ. of South Florida, is shown below:














The diagram shows assorted "orbital elements" for a planet of mass m2 (the Sun is m1) and this can be used as a basis for orbital energy analysis and also to predict future positions. Energy constants in celestial mechanics are very useful for quickly coming to terms with specific properties of an orbit such as shown in the accompanying sketch- designating a generic orbit in x-y-z space. In the diagram, w   is the argument of the perihelion, W is the longitude of the ascending node, f is the true anomaly and i is the inclination of the orbit. The critical or key parameter here is h, the angular momentum vector for the orbiting system. It may be useful here to refer to the diagram (b) in the figure used for angular momentum vector (z) in an atomic system:

Getting specific, assuming r and r' are r (radius vector) and d r/dt, respectively, the magnitude h, of the angular momentum vector is:

h = r x r’ =

(y z’ - z y’)

(z x’ - x z’) = (c1 c2 c3)

(x y’ - y z’)

so (r x r’) = (c1/ h, c2/ h, c3/h)

and inserting variables one finds:

c1/ h = sin W sin (i)

c2/ h = - cos W sin (i)

c3/h = cos(i)

Now since the inclination of Earth's orbit to the ecliptic  (i) is known (23.5 deg) and therefore cos(i) can be determined, then sin(i) can be as well. Also h can be determined, since: h = c3 / cos(i) = (GMm a (1 – e ²) ^{1 /2} where all the constants are known (a = semi-major axis of orbit, e = eccentricity of orbit)

(The energy equation is: ½V ² - m  /r = c, and the c's - energy integration constants- are  found as shown below)  Also, h can be determined, since: h = c3 / cos(i) .  (Also h =  [c1 ² +  c2 ²   +  c3 ²]  ^½)   We also know  W =  11.26 deg.

To fix ideas, W  is one of the key orbital elements (difference between mean anomaly  M and argument of the perihelion,  w )  where M can be obtained from a table based on observations, and w can be obtained using a Fourier expansion of the mean anomaly, M, e.g.:

 w  = M + (2e – e³ / 4) sin M + 5 e²/4 sin 2M + ... etc.

Once W is known, and c1 and c2 and c3 are known, the astronomer is ready to compute the position of a planet, say Jupiter, forty or so years in the future.  The actual calculations are formidable but can be found in a good celestial mechanics text, or using the orbital element formulae can be worked out using a computer.

The energy integration constants  (c1, c2, c3) factor into a host of different problems but let's consider just one relatively simple application, say to compute a planet's velocity at the perihelion and aphelion points of orbit.  The diagram below can be used:


Here,   the points A and P denote the aphelion (farthest point) and perihelion (closest point) to the Sun (S), respectively. We let   V_A, V_P  be the respective velocities at those orbital extrema. As may be deduced here, points A and P are the only ones in the whole orbit for which the velocities are truly tangential or at right angles to the radius vectors for those positions. Consequently, we can write: 

V = (2π/T) r

where r is the radius vector at the point, and T is the period. 

 If Kepler's 2nd (equal areas)  law holds at every point (i.e. equal areas swept out in equal intervals of time) we also have:

r² (2π/T) = h 

where 'h' is a constant ('specific relative angular momentum') which is twice the rate of area description (i.e. by the radius vector). Thus, if the radius vector is r1, then h = 2A1, when A1 = π(r1)². Hence, at aphelion and perihelion only we have: 

V = h/r    for which r = a(1   +  e)  OR:   r =  a (1  -  e)

For the perihelion velocity we have:

V_P = h/ a(1 - e)

where a is the semi-major axis, and e is the eccentricity.

For the velocity at aphelion: 

V_A = h/ a (1 + e)

Then the ratio of velocities is:

(V_P/V_A) = (1 + e)/ (1 - e)

The correct energy equation can be written:

½V² - m  / r =   a 

where a  is an energy integration constant.   More conventionally it is written  without  the energy constant:


V ² = m (2/r - 1/a)


which is also known as the "vis viva" equation, one of the most important in celestial mechanics. 

Since for any bound system of masses m1 and m2, m = G (m1 + m2), where G is the Newtonian gravitational constant (G = 6.7 x 10^-11 Nm²/kg²) then if we know V_P and V_A, along with a and e, we can compute  a, viz. 

a      = ½V_P² - m /a(1 - e) 

at perihelion, and 

a = ½V_A² - m /a(1 +  e) 

at aphelion. 

Consider for the Earth-Sun system :

m  = 1.33 x 10²⁰ Nm²/kg

We already know G and m1= 1.99 x 10³⁰ kg (Sun's mass) and m2 = 6.4 x 10²⁴ kg, (Earth's mass)

Also: a (semi-major axis)  = 1.496 x 10¹¹ m

Then h = + [m a(1 - e²)]^½ =    4.46 x 10¹⁵ N-m/kg = 

4.46 x 10¹⁵ J/kg 

The energy constant a  =   - m/ 2a   for an elliptical orbit

So:   

a  =    - (1.33 x 10²⁰ Nm²/kg) / 2 (1.496 x 10¹¹ m)  =   

-4.45 x 10⁸ m²/s² 

The eccentricity of the orbit e, can now be obtained from:

e =   [1   +   (2 h ²  a )/ m ² ] ^½ =

[1   +    2(4.46 x 10¹⁵ J/kg ) ² (-4.45 x 10⁸ m²/s²)/ (1.33 x 10²⁰ Nm²/kg)² ] ^½ =

0.016

(Which is confirmed on consulting a planetary data table)

What about the velocities at perihelion and aphelion?

Since we have obtained h and e,    the velocity at perihelion is easy to calculate, e.g. 

V_P = h/ a(1- e) =

4.46 x 10¹⁵ J/kg / [1.496 x 10¹¹ m(1 - 0.016)]

V_P = 3.03 x 10⁴ m/s  = 30, 300 m/s

and the velocity at aphelion can be obtained using:

V_A = h/ a(1 + e) =

4.46 x 10¹⁵ J/kg / [1.496 x 10¹¹ m(1 + 0.016)]

V_A = 2.93 x 10⁴ m/s    = 29, 300 m/s

Let's take another illustration from the Pluto-Charon system.   The orbit of Charon with respect to Pluto, has an eccentricity e = 0.0020. The semi-major axis of the orbit is 19, 450 km. The mass of Pluto = 1.27 x 10²² kg and the mass ratio (Charon to Pluto) is found to be m(c)/m(P) = 0.12.  From this  data, find: 

a) The mass of Charon

b) The ratio of the velocity of Charon at perihelion to aphelion

c) The period of Charon, and its velocity 

a) The mass is found using the mass ratio. Since m(c)/m(P) = 0.12 and m(P) =1.27 x 10²² kg , then:

m(C) = 0.12 (1.27 x 10²² kg ) = 1.52 x 10²¹ kg

b) This ratio is obtained from:

(V_P/V_A)  = (1 + e)/ (1 - e)

and we know, e = 0.0020, so:

(V_P/V_A)  = (1 + 0.0020)/ (1 - 0.0020) =

 (1.0020)/ (0.998) = 1.004

Therefore:  (V_P/V_A)  = 1.004

c) The period is obtained from: T = 2π (a³/m)^½

where a = 19, 450 km = 1.94 x 10⁷ m

For preserving unit consistency this is one time we decline to use units of AU and years, and instead use the distance in meters to conform with the required units for G. Then: T =

2π [(1.94 x 10⁷ m)³/[6.7 x 10^-11 Nm²/kg²) (1.422 x 10²² kg)]^½

T = 5.5 x 10⁵ s = 6.37 days


The velocity is:

V = 2πa/T = 2π( 1.94 x 10⁷ m )/(5.5 x 10⁵ s)=

 221.6 ms^-1


This is all meant to show the power of celestial mechanics at even a basic level.   Of course, much much more can be done such as determining the geocentric coordinates of Jupiter on July 20, 2048.  But that exercise is far beyond the scope of this response and requires knowing not only all the relevant orbital parameters of each planet (Earth, Jupiter) but also their mean motions and then using the equations needed such as the expansion of mean anomaly, e.g.

w  = M + (2e – e³ / 4) sin M + 5 e²/4 sin 2M + ...

Precision targeting of a planet for future spacecraft arrival, also necessitates computation of likely perturbations such as shown below:

Perturbations are important because they show how the positions of some given planet are influenced by the gravitational attraction of another. Unless one then reckons the differential location into computations, he will find a spacecraft destined for some destination planet will likely miss it entirely.

In the illustration shown  above we seek to quantitatively estimate the effect of Jupiter on the orbital position of Earth. The respective masses are: m1(Sun), m2(Earth) and m3 (Jupiter) and we assign relative radius vectors that approximate to the actual distance in AU or astronomical units where 1 AU = 149 million kilometers. Thus, r = 1.0 and r3 (for Jupiter) = 5.  Also, D defines the separation between Earth and Jupiter and is the key parameter for estimating the magnitude of the perturbation.

The angle S separating the r and r3 vectors can be anything but for working purposes we use S = 120 degrees, which yields a value: cos(S) = cos (120) =      - ½. This will be useful for computing the first three Legendre polynomials which factor into the kind of generic perturbations I'm looking at in this example.

The first three Legendre polynomials for this context are:

P 0  = 1

P 1 = cos (S)

P 2 = ½ (3 cos ² (S) – 1)

Note that by comparing  P 1  above to  P 1 (x) earlier we can easily discern:

x <-> cos (S)

In other words, x in the Rodrigues’ formula,  e.g.

d/dx [(1 – x ² ) d/dx P n  (x)] + n(n + 1)P n (x)

is now replaced by cos (S) for the purposes of this perturbation estimate.

As we are perturbing an inner planet by an outer one, we are expanding (1/ D) using the Legendre polynomials. Computing the Legendre polynomials for an Angle S = 120 deg, one finds: P 0 = 1,  P 1 = - ½ and  P 2 = -1/8

The perturbing function is then easily computed, viz.:
  
1/D  = 1/ {1 + (1/5) 2    – 2(1/5) (-½)} ^1/2  

And   1/D = 0.18

 is now replaced by cos (S) for the purposes of this perturbation estimate.  

Again, this is all intended as a general illustration of how celestial mechanics works in the situation of sending a spacecraft to a plant.