Document First Posted July 1, 2003
Last Revised June23, 2004

Moments of the Distance from the Force Center in a Two-Body Kepler Orbit
                         (Or "Just How Far Away is the Sun Anyway?")
 

Al Lehnen Jeremy Kessenich
Mathematics Department Mathematics Department
Madison Area Technical College Mount Horeb High School
3550 Anderson Street 53704 305 S. 8th Street
Madison, WI 53704 Mt. Horeb, WI 53572
(608) 246-6567 (608) 437-2400 ext. 2115
alehnen@madisoncollege.edu kessenichjeremy@mhasd.k12.wi.us
http://my.execpc.com/~aplehnen/al.htm

Table of Contents:

Abstract
1.Introduction
2.The Average over Arc Length
3.The Average over True Anomaly
4. The Average over Eccentric Anomaly
5. The Average over Time
6. A Comparison of the Different Methods of Averaging
7. A Comparison between the Moments Averaged over Time and Astronomical Data
8. Conclusion
References
 

Abstract:
This work examines the mean distance from the force center ("sun") at a focus to a planet in a Kepler elliptical orbit. Four different methods of averaging are presented: the average over arc length, the average over the central angle about the force center (the "true anomaly"), the average over the central angle about the center of the ellipse (the "eccentric anomaly"), and the average over time. For each of the averaging methods the analysis actually develops explicit formulas for the mean of the distance raised to an arbitrary real power. These formulas are expressed as functions of the eccentricity of the orbit. Connections and comparisons between the different averages are developed. In particular, for integer exponents greater than -2, the time averaged moment of the distance is always given exactly by a polynomial in the eccentricity squared. Closed form expressions are given for the averages over true and eccentric anomalies for any integer exponents. Results for arbitrary real exponents and the results for the average over arc length are presented as explicit power series in the eccentricity squared. Finally, a comparison is made between the exact time averaged moments of a Kepler orbit and astronomical data obtained from the United States Naval Observatory web site. The agreement between observation and a Kepler orbit model is analyzed and it appears that the difference is a function of the exponent on the distance.
 

1.Introduction

In 1609 Johannes Kepler published his New Astronomy in which he formulated the first two of his famous three laws of Planetary motion [10, p. 122]. This was followed by the third law which appeared in the Harmonice Mundi in 1618. According to the first law the planets follow elliptic orbits with the sun at a focus. Unlike a circle, on an elliptical orbit the distance from the sun to the planet is not a constant. It varies from its minimum value at perihelion to its maximum value at aphelion. A natural question to ask then is "What is the average distance from the sun to the planet?". For the case of the earth this value is often [4, p. 245, 9] stated as 149,597,890 km or 92,955,819 miles. Recently there seems to have been a renewed interest in the problem of the average distance in a Kepler orbit and the definition of the astronomical unit [6, 19, 28]. But what is meant by the "average"? In this work we present four different but related "averages" and generalize the analysis to compute these averages for the distance raised to any power.

One way to analyze the problem is to take the center of the ellipse as the origin and align the x-axis parallel to the major axis so that perihelion has the coordinates (a, 0). The semi-major axis is a and the sun will then have the coordinates  , where is the eccentricity of the ellipse. This is illustrated in Figure 1. From familiar results of analytic geometry the semi-minor axis b is given by  and the ellipse can be represented parametrically by the following equations.

with  . (1)

In the Kepler orbit the parametric central angle E is sometimes called [26, p. 89], the "eccentric anomaly". Designate r as the distance from the focus at  to the point generated by the angle E. An application of the distance formula gives the surprisingly simple result that

. (2)

The ellipse can also be represented in polar coordinates as

with  .

The angle is measured from the x axis about the focus at  and is called the "true anomaly". The distance from the "sun" to the position of perihelion is r0 . The major axis is the sum of r at with r at , so that

.

It follows that  and the distance from the sun to the planet can be written as

. (3)

Figure 1

2. The Average over Arc Length

The mean value of r averaged over "the points" on the ellipse can be calculated by a rather simple procedure. From the symmetry of the ellipse about the x-axis, the average of r over the first and second quadrants equals the average of r over the entire ellipse. Now consider points on the ellipse in the first quadrant generated by central angles  with i ranging from 1 to m. The reflection about the y axis of these points are generated by the central angles  . This is illustrated in Figure 2.

The average value of r over the 2m points above the x-axis is then given by  . The symmetry of the angles and the definition of an ellipse require that  , a result which also follows immediately from equation (2). Thus, the "geometric" average distance from a focus to points on an ellipse is just the semi-major axis a. For this reason some authors [26, p. 87] refer to a as the mean distance of the orbiting body.
 
 

Figure 2

To generalize and extend this result to arbitrary powers of r, the element of arc length along the ellipse is required. From equation (1) this is calculated as

.

The "geometric" average of rn, the n’th moment of r, will thus be defined as the average of rn with respect to arc length.

(4)

It is sufficient to use an upper bound of  in the integrals since the ellipse is symmetric about the x-axis. The denominator in equation (4) can be expressed in terms of a complete elliptic integral of the second kind [27, pp. 517-518, 22]. An alternative formulation is to represent the integral as an infinite series in the eccentricity squared. For the binomial expansion of gives the following absolutely convergent infinite series.

(5)

The "double factorial", defined as  , is used to make formulas more concise. Integration by parts gives the standard reduction formula [20, p. T-3]

, which in turn implies the recursion that  . By an induction argument one concludes that

. (6)

Combining equations (5) and (6) gives the explicit formula,

Here, F is Gauss’s hypergeometric function [1, p. 556, 25, 27, p. 281]. On the domain  the integrals of odd powers of cosine are zero by symmetry. Therefore, the numerator of equation (4) can be expanded as a series where only even powers of cosine appear.

(7) When n is a positive integer this sum will terminate and the coefficient of the integral in the sum over j is simply the binomial coefficient  . Using equations (5) and (6), the integral in the sum over j that appears in equation (7) can also be expressed in terms of a hypergeometric function. Combining these results the n’th moment of r has the following explicit representations.  


 

If n is a positive integer the upper limit on the sum over j is .
Here designates the greatest integer function of x .

An explicit power series in the eccentricity for the nth moment of r would facilitate both the computation and the interpretation of results. To obtain such an expression expand the hypergeometric function that appears in the following series.

Using the Cauchy procedure [2, pp. 384-385] the product of the two power series in can be expressed as a single power series. This gives the following "simplified" result.

, where

.

Thus, a rather long but straightforward calculation has established that
. (8)
Now expand as a series in and multiply the result with the power series in the numerator of equation (8). This gives the explicit formula that , (9)
with
  This equation reproduces the previously derived result that  . As another example in the use of this formula consider the variation in the values of rn as the planet moves around its elliptical orbit. The variance is a measure of this spread and it is computed as
.
From equation (9) this has the approximation that . Thus, the standard deviation of rn can be computed to lowest order in the eccentricity by the formula
. (10)
Actually, this result is not at all surprising. From equation (2),  . So equation (10) can be interpreted as the familiar relation between the root mean square and amplitude of a sine wave.

For the special case of n = 1, a more careful use of equation (8) gives the result that

. This means that the standard deviation of r averaged over arc length is approximately the focal length of the ellipse divided by the square root of 2.
 

3. The Average over True Anomaly ( )

A second way to calculate the mean of rn over an ellipse would be to sample rn uniformly over the angle . From equation (3) this average is computed by the following formula.

Again the integral is defined only on because of the symmetry of the ellipse about the x-axis. For any real number p let the integral be defined as

, (11) so that . (12) If n is a positive integer the binomial theorem and equation (6) can be used to obtain a polynomial in the eccentricity squared for  .
(13)

The problem of evaluating  when n is a positive integer is solved by noting some fundamental symmetries. First the substitution  demonstrates that  . A second, less obvious symmetry is stated as a theorem.

Theorem 1: For any real number p, .

Proof: Consider for the transformation (motivated by the Weierstrass substitution for rational functions of sine and cosine [20, p. 570])

.

The differential is then given by  . The standard trigonometry identities that  and  give after some manipulation the following rather simple differential relation.

Now to evaluate  the expression  needs to be expressed in terms of  z.

Hence,  . The limits of integration transform as follows:

.
All of the "pieces are assembled" below.

Thus, the transformation confirms the stated symmetry.

From Theorem 1 and equation (13) for any positive integer n .(14)

From equations (12) through (14) the angle averaged moments of the distance from the focus for any positive integer n reduce to the following closed form expressions.

(15)

(16)

When n is not a negative integer the binomial expansion and equation (6) can be used to express as an infinite series in the eccentricity squared.
This gives a general formula for the moment averaged over true anomaly valid for any real value of n.

(17)

In addition for any real n, Theorem 1 and equation (12) imply that

. (18)

The hypergeometric function, although not required, can be used to derive these same results. Indeed,  can itself be expressed as a hypergeometric function.

 Equation (17) now becomes
.

In this formulation, Theorem 1 and the symmetry expressed in equation (18) are just applications of the Euler transformation
[1, p. 559, 8, p. 1069, 24].

The results for small integer values of n are shown in Table 1. These formulas are not new. Indeed, some of them can be found in the existing literature [16, p. 54, 17, p. 63, 23, 28]. However, it appears that the complete and rather simple characterization of  for integer n is not well known.

Table 1: Mean Moments of the Distance on an Ellipse from the Focus Averaged over the Angle


 


Moment
n = 1
n = 2
n = 3
n = 4

For the special case of  n = 1, the mean value of r is just the semi-minor axis,  .The mean of r averaged over being less than a is certainly reasonable. Measured about the focus, half of the angles are on one side of the latus rectum where r is less than  . Thus "small" values of r are "over-sampled" when compared to the average taken over arc length. The standard deviation of r averaged over true anomaly is given by the formula

.

To lowest order in the eccentricity this is the same answer as the average over arc length.

4. The Average over Eccentric Anomaly (E)

Probably the easiest way to calculate the mean of rn is to sample over the angle E. From equation (2) and the symmetry of the ellipse about the x-axis this average is computed by the formula,

.

From equation (11) and the symmetry that , it follows that

.

Theorem 1 implies a fundamental symmetry between the mean of r -n and the mean of r n-1.

(19)

When n is any positive integer this symmetry and the binomial theorem establish the following formulas.

For any real value of n the infinite series in for results in the following expansion
(20)


The relations between the average over eccentric anomaly and the average over true anomaly

are obtained from equations (12) and (19).

(21)
The results for small integer values of n are shown in Table 2. Smart gives the results for n = -1 and n = -2 [15, p. 36], while the n = 1 case is a worked example in Szebehely’s text [16, p. 53].

Table 2: Mean Moments of the Distance on an Ellipse from the Focus Averaged over the Angle E


 


Moment
n = 1
n = 2
n = 3
n = 4
1

The mean value of r averaged over eccentric anomaly is identical to the mean value of r averaged over arc length. The standard deviation of r averaged over eccentric anomaly is given exactly by

.

5. The Average over Time

For most people, once they have thought about it, the "average" distance of the earth from the sun would probably be calculated by measuring the distance every day for one full year and then dividing the sum of these measurements by the number of days in a year. That is, the average would be "over time", not arc length or angle about the sun. The result of a time average depends on the dynamics of the problem. If the earth orbited the sun at constant linear speed, then the average over time and the average over arc length would be identical. Similarly, if the earth orbited the sun at constant angular speed, then the average over time and the average over would be the same. However, for planets in elliptical orbits neither of these speeds is a constant.

The dynamics this work assumes is that of two spherically symmetric mass distributions interacting with each other through the inverse square law force of gravitational attraction. As is shown in many standard texts this results in the distance between the centers of the two bodies following the orbit of a conic section [7, pp. 58-80, 20, pp. 857-865, 26, pp. 76-89]. For planetary speeds at perihelion which are less than the escape velocity, , the orbits are ellipses (this is Kepler’s First Law) with an eccentricity given by

. (22)

Here is the speed of the orbiting planet at perihelion, G is the "universal gravitational constant", M is the mass of the sun (or the "attractor") and m is the mass of the planet. The presence of m may come as a surprise to some readers. It is due to the separation of the motion of the center of mass from the two-body problem [7, pp. 58-59, 26, pp. 76-77]. This replaces the mass of the planet, m, in the Newton’s Second Law (F = ma) with the "reduced mass"  . Since the sun is more than 300,000 times more massive than the earth, for the earth-sun system M + m is often replaced by just M .

Suppose along the orbit there is a quantity, W, that can be written as an explicit function of  , i.e., , then the time-average of W over one full revolution is given by the integral.

Here T is the period or time required for one full revolution. It is given by the formula

. (23)

This is really just the statement of Kepler’s Third Law. Since the equation of the elliptical orbit is given by equation (3) it is convenient to transform the variable of integration from t to  .

The gravitational force is an example of a "central force". This requires that the quantity  be a constant. This is in essence Kepler’s Second Law. The initial conditions at perihelion lead to the result that  . From this relation and equation (23) calculating the time average of rn reduces to evaluating the following angle integral.

Using the eccentricity as given by equation (22) results in a further simplification.

     (24)

Thus, for a Kepler elliptical orbit there is a fundamental symmetry between the moments of r averaged over time and the moments of r averaged over true anomaly. This symmetry is summarized by the following brief formula, which is valid for any real exponent n.

(25)

In turn, this relation and equation (21) connect the moments of r averaged over time and the moments of r averaged over eccentric anomaly.

The symmetry of the average over true anomaly stated by equation (18) and the connection between the average over time and the average over true anomaly stated in equation (25) combine to reveal a corresponding symmetry present in the average over time.

Pasternack demonstrated in 1937 that there is an analagous symmetry between the expectation values of rn and r-(n+3)  for the solutions of the non-relativistic Schrodinger equation for the hydrogen atom [12].

From equations (15) and (25) if n is an integer and  , then

. (26)
As this formula makes explicit, when n is an integer greater than negative one, the time averaged moment of rn is always a polynomial in the eccentricity squared [5]. The simple form of this general result does not however seem to be widely known.

From equations (16) and (25) if n is an integer and  , then

, (27) while, for n any real number, using equations (17) and (25) gives

. (27)

The resulting formulas for some small integer values of n are shown in Table 3.

Table 3: Mean Moments of the Distance on an Ellipse from the Focus Averaged over Time ( t )


 


Moment
n = 1
n = 2
n = 3
n = 4
1

Some of these results can be found in the literature [11, p. 51, 13, p. 38, 15, p. 41, 16, p. 54, 17, p. 63, 23, 28]. It should be noted that in most of these references the results would be inferred not as a time average, but rather as an average over a variable called the "mean anomaly" (for a definition see:  [11, p. 14, 13, pp. 24-25, 15, p. 18, 16, p. 47, 18, p. 152, 26, p. 89]).

However, as Serafin [14] has noted, the expectation value of a quantity as a function of time is identical to the expectation value of that same quantity as a function of mean anomaly. Surprisingly, the most concise and complete treatment of the time average of rn we found is given not in an astronomy text, but in a treatise [3, pp. 143-145] by Max Born on the “Old Quantum Theory”. If the complete and rather simple characterization of  as stated in equations (26) through (28) are not new, then they are at least not very well known.

For the special case of n = 1, the mean value of r is  . This is larger than a, the average over arc length. Since the planets move faster near the sun and slower when far away, they spend more time at long distances than at short distances. Thus, the time average must be somewhat larger than the spatial average. The desire to quantify this difference between the time average of r and a was the original motivation of this entire analysis.

The standard deviation of the time average of r is given by the formula

.

Thus, to first order in the eccentricity the four different ways of averaging all give the same result (the focal length over the square root of two) for the standard deviation of r. The corresponding equivalence for rn will be demonstrated in section 6.

The result in Table 3 for  is a specific example of a more general result for dynamical systems known as the "Virial Theorem"[7, pp. 69-71]. The total energy is the sum of kinetic and potential energy terms and is a constant along a given orbit. Evaluating the total energy, ET, at perihelion and using equation (22) gives the following result.

For a force that varies as the inverse square, the Virial Theorem states that the time-averaged potential energy is twice the total energy. Hence,  or  .
 

6. A Comparison of the Different Methods of Averaging
 

The original motivation for this work was the average distance (n = 1) from a planet to the sun. Table 4 summarizes the results already presented for this important case.

Table 4: Means and Standard Deviations of the Distance on an Ellipse from the Focus


 


 
Mean of 
Standard Deviation of 
Average over Arc Length
1
Average over True Anomaly
Average over Eccentric Anomaly
1
Average over Time

Quite a few relations between the mean moments averaged over true anomaly, eccentric anomaly, and time have already been derived. However, it is very easy to get them confused. Table 5 provides a concise and explicit summary of the important exact symmetries of these three averaging methods.

Table 5: Exact Symmetries of the Moments of the Distance on an Ellipse from the Focus


 


An interesting application of the four averaging methods is the location of the mean position of a planet in a Kepler orbit. As elementary calculations based on equation (1) readily confirm, the mean position when averaged over either arc length or eccentric anomaly is the origin, i.e., the center of the ellipse. The average over true anomaly requires a little more effort.


The average value of y vanishes due to the symmetry of the ellipse about the x-axis, while the following trick establishes the average value of x.
So,
.

Hence, for small eccentricities, the mean position of a planet averaged over true anomaly is approximately half way between the center of the ellipse and the focus (sun). Equation (24) provides the starting point to calculate the mean position averaged over time.


Again by symmetry, the average value of y is zero. Repeating the same trick as before computes the average of x.

For a Kepler orbit the time-averaged position is exactly half way between the center of the ellipse and the “empty” focus opposite the force center. This result is given both by Plummer [13, p. 39] and Born [3, p. 145]. The distance between the mean position averaged over true anomaly and the mean position averaged over time is larger than the focal length of the ellipse but approaches the focal length for small eccentricities.

Further comparisons among the four different averages will be made at small values of the eccentricity so that terms of order and higher can be neglected. This is done for the following three reasons.

                It makes the analysis particularly easy.
                It is relevant. In our solar system all of the planets have orbits with eccentricities less than 0.25 [9, 20, p. 864].
                Differences between the three of the four methods of averaging are already apparent in the terms of order  .

For any real value of n the moment equation

, (29)

follows from equation (9). An identical result is obtained from equation (20) for the average over eccentric anomaly. The binomial expansion applied to equations (17) and (28) gives

(30)

and

. (31)

From these equations a quick calculation of  confirms that all four methods of averaging give a variance of rnequal to  and a standard deviation of  given by . This implies that for small eccentricities  and  have the same standard deviation.

Comparing equations (29) and (31) gives an approximate symmetry between the time average and the average over arc length.

.

Since the average over arc length and the average over eccentric anomaly agree to order , this result is just a restatement of the exact symmetry that exists between the time average and the average over eccentric anomaly.

The average of rn over "underestimates" the "true" or arc length average of rn , while the average of rn over time "overestimates" the answer. Therefore, one might expect the averages of these two estimates to be close to the arc length average. Using all three of the equations (29) through (31) demonstrates that this is indeed the case.

.

We therefore have the truly "tongue twisting" conclusion that for small eccentricities the average of the average of rn averaged over true anomaly and the average of rn averaged over time is the average of rn averaged over arc length!
 

7. A Comparison between the Moments Averaged over Time and Astronomical Data

A final point of interest is to compare  calculated for a Kepler orbit with observations of the earth’s actual motion about the sun. Using the MICA (Multi-Year Interactive Computer Almanac) located on the United States Naval Observatory web site [21], the daily distance from the earth to the sun in astronomical units (AU) was computed for the calendar year June 28, 2002 through June 28, 2003. On each date calculation was performed for 00:00 UT (“Universal Time”). Figure 3 displays the results of these calculations.

Figure 3

Yearly means and population standard deviations for  were computed using the 365 values of r from June 28, 2002 through June 27, 2003. To perform these calculations it is necessary to determine a value for a, the semi-major axis of the earth’s orbit. It is often claimed that a = 1 AU. However, there exists some confusion (certainly on our part) on this point. In his text Taff gives a value of a = 1.00000023 AU [18, p. 191]. The comprehensive Allen's Astrophysical Quantities states both that a = 1 AU [4, top of p. 12] and that a = 1.00000105726665 AU [4, bottom of p. 12]! These differences are of course quite slight, but for small eccentricities the discrepancy between  and 1 will also be small. Finally, we decided that the best approach was to be self-consistent within the MICA calculations. A numerical search of MICA determined that aphelion occurred on July 06, 2002 at 03:27 UT and perihelion occurred on January 04, 2003 at 04:53 UT. The computed values of r at these times give a = 1.000004292 AU and this is the value we used.

In order to compute for a Kepler orbit the eccentricity must be specified. The value for the eccentricity of the earth’s orbit found in many elementary texts is 0.0167 [20, p. 864]. Reference 4 quotes the extremely precise value of  [4, p. 245], while the JPL’s web site states a value of 0.01671022. We eventually settled on  for reasons that we discuss below. Of course, all of these values give essentially the same results for the average of rn over time of a Kepler orbit. The computed values are displayed in Tables 6 and 7 for non-zero integer values of n between – 4 and 4.

Table 6: Comparison for the Earth of MICA Calculations with results of a Kepler Orbit


 


 
n = 1
n = 2
n = 3
n = 4
Yearly Mean of 
1.0001199916
1.0003793140
1.0007779881 
1.0013160643
Yearly Standard Deviation of 
0.0118032332
0.0236070423 
0.0354145115
0.0472287252
1.0001393415
1.0004180244
1.0008360780
1.0013935604
0.0118034767
0.0236079815
0.0354165982
0.0472324107

 

Table 7: Comparison for the Earth of MICA Calculations with results of a Kepler Orbit


 


 
n =– 1
n =– 2
n =– 3
n =– 4
Yearly Mean of 
 1.0000193471
1.0001780702
1.0004762355
1.0009139385
Yearly Standard Deviation of 
0.0118057423
0.0236170805
0.0354371035
0.0472689034
1
1.0001393706
1.0004181701
1.0008364859
0.0118055330
0.0236162080
0.0354351128
0.0472653379

Let  represent the yearly average of  based on the MICA calculations and represent the fractional deviation between this mean and  .
The values shown for both  and  are close to 1, so that the values of  are quite small. It is really the agreement between  and  that measures the how well the Kepler two-body orbit fits the MICA results. From Tables 6 and 7, this agreement, while fair, is not entirely convincing. In contrast, the agreement between the standard deviations of the Kepler orbit and those of the MICA calculations seems much better, but then this is to be expected. The standard deviation of the MICA results is given by

.
In terms of the alpha function this can be expressed as
.
Since the alphas are small, the observed standard deviation will be very well approximated by the Kepler orbit result,
,
which in turn is well approximated by  .
The differences between the values of  and  show significant non-random discrepancies between the Kepler orbit time-averages and the MICA yearly averages. A graph of  for integer values of n ranging from –5 to 6 was generated and is shown in Figure 4. A least squares linear fit of the results is also shown.

Figure 4

The agreement should be exact at n = 0. Hence, the value of the eccentricity was varied to make the regression intercept vanish. This “fitted” value of the eccentricity worked out to be 0.0166938. To the numerical accuracy of the calculation this value of the eccentricity gave an “exact” linear fit with a correlation coefficient equal to –1. The calculated slope of the regression line is –0.0000193469866 . Curiously enough, for this value of the eccentricity and n between –5 and 6, a value of
a = 0.999984944700 AU brings the Kepler orbit time-averages and the MICA yearly averages into near exact agreement, i.e., to within a relative error of  or less. However, this value of a is not self-consistent within the MICA calculated r’s which clearly have a semi-major axis greater than 1 AU.

The linear or “near-linear” deviation of the Kepler orbit results from the MICA calculations appears rather startling. We therefore attempted to understand it better. The average given by a definite integral assumes continuous and uniform sampling over the domain of integration. The average of the MICA results was based on 365 discrete scores. Thus, there is a “continuity” correction. This could be reduced by say downloading MICA calculations for every hour in a year rather than every day, but we were eager to avoid this complication! We therefore decided to turn the problem around. Rather than approximate the discrete average by an integral, we decided to use the discrete MICA values of r to perform a numerical quadrature of  and then compare this result directly against  . Since a full year needs to pass for the calculation of  , two additional MICA calculations are needed. One is obviously for June 28, 2003, the second is a bit more problematic. In the spirit of self-consistency it would be nice to define T directly from the MICA results. Unfortunately, there are already enough deviations from a Kepler orbit to cause r to be non-periodic over a mean tropical year. The MICA calculations have r at 22:40 on June 25, 2003 matching the value of r at 00:00 on June 28, 2002. We therefore used the published value that the number of days in a year is given by [7]. Therefore, the right end point of the domain of integration requires r to be calculated at 06:09 on Jun 28, 2003. Thus, to perform the numerical integration there are 367 discrete values of r computed by MICA. For j less than 366 define Rj  to be the value of  at  and let R366 be the value of  at t = T .
The trapezoid rule approximation to  is then given by the formula,

.

A Simpson’s rule approximation requires the messy job of fitting a parabola through the last three R values on the last two unequal intervals,  and  . The approximation to the integral is the rather complicated looking formula shown below.

The results of the two approximations are displayed in Table 8. From the convergence of the two methods neither increased sampling of r values nor more sophisticated quadrature methods would significantly improve the answers.

Table 8: Numerical Approximations to  for the MICA Generated Values of r


 

Approximation n = 1 n = 2  n = 3 n = 4

Trapezoid Rule
1.0001316895 1.0004028096 1.0008133848 1.0013634686

Simpson’s Rule
1.0001316896 1.0004028098 1.0008133852 1.0013634691

Trapezoid Rule
1.0000077460 1.0001549613 1.0004417092 1.0008680819

Simpson’s Rule
1.0000077459 1.0001549611 1.0004417089 1.0008680815

The fractional deviations of  from  for integer values of n ranging from –5 to 6 are displayed in Figure 5. The discrepancies have been reduced by slightly more than a factor of 2 from the discrepancies associated with the discrete averages. Nevertheless, there remains a nearly linear (the correlation coefficient is -0.999823639 ) deviation of the MICA calculations from the Kepler two-body results. We are uncertain whether this is a “real” linear deviation present in the actual orbit of the earth, an artifact of the MICA calculations, a flawed assumption in our treatment of the results of the MICA calculations, or a combination of these factors. In any case, it can not be explained away as just a “continuity” correction.

Figure 5

Despite our fixation in trying to explain small differences, we were actually impressed by how well the Kepler orbit averaged over time models the MICA results (reality?). After all, the Kepler two-body solution ignores the gravitational forces between the planets and between planets and their satellites. Furthermore, the relevant mass distributions are not perfectly spherical. Eventually even general relativity makes corrections [26, p. 389].
 

8. Conclusion

This work has revisited one of the most famous and influential problems in the history of science. The primary question, "What is the average value of the distance between a planet and the sun raised to a power?", is both easy to state and to understand. Surprisingly, within the constraints of the Kepler two-body orbit, the problem can be satisfactorily solved using only the methods typically taught in the introductory calculus sequence. In particular, the systematic characterization of the time averaged moments of rn for integer n seems especially simple and appealing. Similarly, the symmetries inherent in the four different methods of averaging do not, initially at least, appear self-evident. The precise statement and demonstration of these symmetries would by itself seem to be a topic of interest. Given the importance and reputation of the Kepler problem it would be extremely surprising if any of these results are really new. Nevertheless, the detailed solutions presented here may help motivate introductory calculus students and their instructors. There is always a need for interesting but tractable examples. Just the existence of multiple ways of calculating the average distance is an “eye-opener” for some students. Finally, the systematic examination of the function, while a much more challenging exercise in celestial mechanics than this study, might be a problem worth further investigation.
 

Acknowledgements

We would like to thank Professor Charles Goebel of the University of Wisconsin-Madison Physics Department for many helpful comments. In particular, it was Professor Goebel who pointed out that our Theorem 1 is a special case of one of the Euler transformations of the hypergeometric function. He also suggested to us the problem of the mean position. The reference to the work of Max Born was kindly provided by Professor Turgay Uzer of the Physics Department at Georgia Institute of Technology. Professor Lorenzo Curtis of the Department of Physics and Astronomy at the University of Toledo provided references to his own work as well of that of Pasternack and Williams.
 

References
 

1.   M. Abramowitz and I. Stegun, editors, Handbook of Mathematical Functions with Formulas, Graphs, and
      Mathematical Tables, Tenth ed., National Bureau of Standards Applied Mathematics Series - 55, 1972.
      Several online versions of this reference are freely available. One is at
      http://www.convertit.com/Go/Photonicsonline/Reference/AMS55.ASP .

2.   R. Bartle, The Elements of Real Analysis, John Wiley, 1964.

3.    M. Born, J. Fisher, and D. Hartree, The Mechanics of the Atom, London, G. Bell and Sons, 1927.

4.   A. Cox, editor, Allen's Astrophysical Quantities, Fourth ed., Springer-New York, 2000.

5.    L. Curtis, Classical mnemonic approach for obtaining hydogenic expectation values of rP, Physical Review A, 43, 568-569 (1991).

6.    D. Deever, The Average Distance of the Earth from the Sun, College Math Journal, 30, 218-220 (1999).

7.   H. Goldstein, Classical Mechanics, Addison-Wesley, 1965.

8.   I. Gradshteyn, J. Ryzhik, and A. Jeffrey, editor, Table of Integrals, Series, and Products, Fifth ed., Academic Press,
      1994.

9.   The Jet Propulsion Laboratory of the California Institute of Technology web sites:
      http://sse.jpl.nasa.gov/planets/charchart.cfm and http://ssd.jpl.nasa.gov/astro_constants.html

10.   A. Koestler, The Watershed, Anchor Books, 1960.

11.   J. Kovalevsky, Introduction to Celestial Mechanics, D. Reidel Publishing, 1967.

12.    S. Pasternack, On the Mean Value of rs for Keplerian Systems, Proceedings of the National Academy of  Science USA 23, 91-94 (1937).

13.   H. Plummer, An Introductory Treatise of Dynamical Astronomy, Cambridge University Press, 1918.

14.   R. Serafin, The Mean Anomaly in Elliptic Motion as Random Variable, Celestial Mechanics, 21, 351-356 (1980).

15.   W. Smart, Celestial Mechanics, Longmans Green and Company, London, 1953.

16.   V. Szebehely, Adventures in Celestial Mechanics, A First Course in the Theory of Orbits,
        University of Texas Press, Austin, 1989.

17.   L. Taff, Celestial Mechanics: A Computational Guide for the Practitioner, John Wiley, New York, 1985.

18.   L. Taff, Computational Spherical Astronomy, Wiley Interscience, New York, 1981.

19.    D. Teets, Transits of Venus and the Astronomical Unit, Mathematics Magazine, 76, 335-348 (2003).

20.   G. Thomas with revisions by R. Finney, M. Weir and F. Giordano, Calculus Early Transcendentals, updated Tenth
        Ed., Addison-Wesley, 2003.

21.   The United States Naval Observatory web sites: http://aa.usno.navy.mil/ and
      http://aa.usno.navy.mil/cgi-bin/aa_micaform2?calc=3&ZZZ=END .

22.   E. Weissteins, World of Mathematics, The Complete Elliptic Integral of the Second Kind entry
        at the Wolfram Research web site:
      http://mathworld.wolfram.com/CompleteEllipticIntegraloftheSecondKind.html .

23.   E. Weissteins, World of Physics, The Elliptical Orbit entry at the Wolfram Research web site:
      http://scienceworld.wolfram.com/physics/EllipticalOrbit.html . As of this date, December 2003, there are discrepancies
      between the formulas we derived for the average of  rn over true anomaly and some of those stated on this web site.

24.   E. Weissteins, World of Mathematics, The Euler’s Hypergeometric Transformations entry at the Wolfram Research web
        site: http://mathworld.wolfram.com/EulersHypergeometricTransformations.html .

25.   E. Weissteins, World of Mathematics, The Hypergeometric Function entry at the Wolfram Research web site:
      http://mathworld.wolfram.com/HypergeometricFunction.html .

26.   E. Whittaker, A treatise on the Analytical Dynamics of Particles and Rigid Bodies, Fourth Edition, Cambridge
        University Press, 1965.

27.   E. Whittaker and G. Watson, A Course of Modern Analysis, Cambridge University Press, 1969.

28.    D. Williams, Average distance between a star and a planet in an eccentric orbit, American Journal of Physics, 71, 1198-1200 (2003).