Celestial Mechanics

Author : J. B. Calvert
Pure Mathematical Physics
 Contents

i. Introduction
ii. Newton's Laws
iii. Keplerian Orbits
iv. Coordinates and Orbits
v. Cometary Orbits
vi. Determination of Orbits
vii. Earth Satellites
viii. Other Orbit Lore
ix. References

Pure Mathematical Physics

 1

Introduction
The grand title for this article conceals the fact that
all I want to do here is to show how to use published
orbital elements to find the location of a planet in the
solar system, and to provide a good explanation for
what is involved in the process, including the
elements of Newtonian mechanics. I have derived all
the results used, with one or two minor exceptions,
but it is not necessary to understand the derivations
to use the results. However, it certainly helps to
know what you are doing.
Johannes Kepler (1571-1630) discovered that the
planets moved in elliptical orbits, and his three laws
permit the calculation of planetary position once the
orbit is known. Although correct, the loss of ideal
circular motion as a fundamental (though impossible
to apply) concept was a disappointment. This was
changed into beautiful triumph in 1686 by Isaac
Newton (1642-1727), who showed that the observed

Pure Mathematical Physics

 2

orbital motion was a consequence of fundamental

principles of universal application.
Newton's formulation of mechanics, which involved
the new concepts of mass and force, was subjected
to intense, but sterile, criticism by Ernst Mach
(1838-1916) and others, which did not change its
application to the slightest degree, and shed no light
on its fundamentals. The ideas of Albert Einstein
(1879-1955), on the other hand, modified the
fundamentals essentially, especially in the theory of
gravity and the kinematics of motion, making the
theory able to explain the smallest discrepancies
with the earlier predictions. This theory of relativity
has also been abundantly proved by observation.
Except for Newton and Einstein, all the vast amount
of work in mechanics has been the elaboration and
application of the theory, which has proved in every
respect correct. Particularly notable is the work of
Joseph Louis Lagrange (1736-1813) in The
Beautiful Theory, where forces are replaced by work
and energy, and the laws of mechanics can be
expressed as problems in maxima and minima. This
branch of mechanics is usually called analytical
mechanics, and mechanics with forces is called

Pure Mathematical Physics

 3

vector mechanics, from the convenient use of that

formalism.
Newton's theory of motion, in its original form, is
still the basis of engineering mechanics, and a
fundamental part of engineering training. Indeed, it
is still used in celestial mechanics. The only aspect
needing relativistic corrections so far has been the
time effects in the Global Positioning System. In the
solar system, velocities are very comfortably less
than the speed of light, 3 x 108 m/s, although
distances are great. Light requires about 8 minutes to
pass from the sun to the earth, and gravity travels at
the same speed. However, the sun does not care a
great deal about the force exerted on it by the earth,
while the earth moves in a static gravitational field,
so the effects of the finite speed of gravity are not
evident. These effects have, however, recently been
measured and verified, so it is known that gravity is
not instantaneous. We will not need to consider any
relativistic effects in the sequel, leaving these tiny
effects to another article.

Pure Mathematical Physics

 4

Newton's Laws
Newton's laws of motion are: (1) a body under no
forces moves uniformly in a straight line; (2) the
acceleration of a body is equal to the force acting on
it divided by its mass; and (3) the mutual forces
acting between two bodies are equal and opposite.
Newton stated the laws in this form so that these
axioms would use only common notions and not
involve any results following from them, in the spirit
of Euclidean geometry. I have used different words
in an attempt at conciseness, but the import is the
same.
The first law gives an operational method for
determining an inertial system, a way of describing
the position of a mass point in space as a function of
time. In the absence of gravitating masses (such as
the earth) it is found that a Euclidean space is such a
system, in which a uniform time can be
approximated by periodic events (such as the

Pure Mathematical Physics

 5

oscillations in lasers). Strictly (and uselessly)

speaking, the definition is circular, but in practice it
is clearly possible to establish inertial systems at
least approximately, and the results obtained agree
perfectly with observation.
The second law defines force as the cause of a time
rate of change of velocity, by dv/dt = f/m. Again, the
definition is logically faulty, but can be realized with
great accuracy in practice. Whatever causes a change
in velocity in an inertial system is a force. Masses
can be compared dynamically, or more easily by the
gravitational forces on them. A consistent array of
masses and forces can be assembled, which always
gives the observed motion. When the velocity is
known, the position can be found from the
kinematical relation dr/dt = v. Relations not
involving force we call kinematical, those involving
force, dynamical.
The third law guarantees that the total linear
momentum and total angular momentum will be
conserved; that is, that they will remain constant
during the motion. These quantities are defined later,
and their conservation is the motive for the third law,
which cannot be stated in a general form in other

Pure Mathematical Physics

 6

than an extremely tedious and unilluminating way.

Of course, not all forces obey it (magnetic forces
between currents are an example), but nevertheless
momenta are conserved. The conservation of
momentum is an important fundamental of
mechanics; the third law is a way of introducing it in
terms of forces. In celestial mechanics, all the forces
are central (act along the line joining two interacting
masses) and so the simple statement is sufficient.
The path traced by the tip of the velocity vector v as
time elapses is called the hodograph, and the path
traced by the tip of the position vector r is called the
trajectory. An orbit is a special kind of trajectory, a
closed (or almost closed) curve. For more
information on the hodograph, see The Hodograph.
These curves are useful in describing, visualizing
and analyzing the motion.
The momentum is the product p = mv. The second
law can be written dp/dt = f, so that force is the rate
of change of momentum. If two particles m and m'
interact, then f' = -f from the third law, so that
(d/dt)(p + p') = f + f' = 0, or dP = 0, where P = p + p'
is the total momentum. The two particles, taken as a
whole, constitute an isolated system, on which no

Pure Mathematical Physics

 7

net force acts, so the rate of change of momentum is

zero, and thus the momentum is
constant.
The angular momentum with respect to
a point O is the vector J = r x p, where r
is the radius vector from O to any point
on the line of action of p = mv. In the
figure, J points out of the page. If we
use r x p instead of p in the equations of
the preceding paragraph, we can prove
that the total angular momentum is conserved. The
conservation of linear and angular momentum gives
us six scalar constants of the motion, useful in
analyzing and understanding the motion. The
existence of these constants is no trivial matter, but a
beautiful and fundamental result of the theory.
The second law can be formally integrated to give p
= ∫ fdt = I, where I is a vector called the impulse, the
integral of the force with respect to time. By taking
moments about point O, we obtain a similar
expression for the angular impulse. Of course, the
limits on the integrals must be handled properly.

Pure Mathematical Physics

 8

If the second law is scalar-multiplied by v, we find

mv·(dv/dt) = f·v, or (d/dt)(mv2/2) = f·v = -dW/dt.
This introduces two new scalar quantities, the kinetic
energy T = mv2/2, and the work -W = ∫ f·dr. I am
sorry for the minus sign, but it simply means that
work is being done on the particle, while W is
conventionally the work done by the particle. If W is
a function of position, then in vector notation f = -
grad W. In one dimension, f = - dW/dx, which may
be more familiar. What we have shown is that
(d/dt)(T + W) = 0, or T + W = constant, a new
conservation law. If W is a conservative quantity
(that is, its value is independent of path and depends
only on position), it is called the potential energy V.
Then, the total energy E = T + V is conserved in the
motion. Energy is a widely used and misused
quantity that could be discussed in great detail, but
only its bare definition is required here, and the
reader is assumed to be familiar with it.
The introduction of energy permits the use of
generalized coordinates, and the derivation of
equations of motion by the Lagrange procedure,
which uses the Lagrangian function L = T - V. This
facilitates the solutions of very many problems,
since we are liberated from carrying around a basket

Pure Mathematical Physics

 9

of vector components. For more information, follow

the link in the Introduction to the article
on Lagrangian mechanics.
For vector mechanics, it is convenient
to have expressions for the velocity and
acceleration components in polar
coordinates. These are derived in the
References, but will be summarized
here. v = (dr/dt)r' + (rdθ/dt)θ'. r' and θ' are unit
vectors pointing radially outward, and tangentially in
the anticlockwise sense, with respect to the position
vector r, as shown in the diagram. The rectangular
unit vectors i and j are also shown. Note that dr'/dθ
= θ' and dθ'/dθ = -r', while the derivatives with
respect to r are zero. dv/dt = [d2r/dt2 - r(dθ/dt)2]r' +
(1/r)(d/dt)(r2dθ/dt)θ'.
The reader should carefully note the distinction
between, for example, d2r/dt2 and d2r/dt2. The first is
the second derivative of a scalar quantity, a simple
thing. The second is a vector, and since the direction
of r can change, the complete expression is
somewhat elaborate, as we have just seen. It is
essential to distinguish carefully between scalars and
vectors, which are arrays of scalars.

Pure Mathematical Physics

 10

These relations are very useful. For example,

suppose that the force is central, so that the
tangential component of the acceleration is zero.
This means that (d/dt)(r2dθ/dt) = 0, or r2dθ/dt =
constant = h. Now, half of r times r dθ is the area dA
swept out by the radius vector in dt, so h is 2dA/dt.
Also, mrdθ/dt is the component of the momentum
perpendicular to the radius vector, so mr2dθ/dt = mh
is the angular momentum. We then have that the
angular momentum L = mh is a constant. Also, this
means that the motion takes place in a plane, and L
is perpendicular to this plane. We have found out a
lot here with very little effort. The constancy of areal
velocity is one of the most useful properties of
orbital motion, not just a curiosity. It relates the rate
of change of the angle θ to the radius r.

Pure Mathematical Physics

 11

Keplerian Orbits

Kepler's three laws of planetary motion, deduced by

prolonged and tedious consideration of the observed
position of Mars, are: (1) the planets move in
ellipses with the sun at one focus; (2) the areas swept
out by the radius vector in equal time intervals are
equal; and (3) The cubes of the mean distances (half
the major axis of the orbit) are proportional to the
squares of the periodic times. These laws are
sufficient to determine the position of a planet at any
later time if its position is known at one time, and
the dimensions and orientation of the orbit are
known.
For a deeper understanding, and the power to attack
an arbitrary problem in orbital motion, such as the
movement of earth satellites, we should consider the
dynamics, on the basis of Newton's theory. The
force acting between two spherical, radially
symmetric bodies of masses M and m a distance r
apart is, by Newtonian gravitation, GMm/r2 in

Pure Mathematical Physics

 12

magnitude, and is directed along the line joining the

centres of the two bodies. G is the Newtonian
gravitational constant, 6.67259 x 10-11 m3/kg-s.
Since the force is central, the angular momentum is
conserved, and the bodies revolve about one another
in a plane with the centre of mass (the barycentre)
fixed. The motion can be analyzed as the rotation of
a reduced mass Mm/(M + m) about the centre of
gravity. The potential energy V is -GMm/r, and is
negative because the particles attract. For the total
energy to be zero, T = mv2/2 = GMm/r, or v =
√(2GM/r). This is the escape velocity from the
particle of mass M at the distance r. At the surface of
the earth, the escape velocity is 11,200 m/s or 25,200
mph. This is also the velocity of a meteor that comes
in from infinity. The escape velocity from the sun at
the distance of the earth is 42,000 m/s or about
95,000 mph. It is not easy to get away
from the earth; it is even harder to
escape from the sun.
If we consider a circular orbit of radius
a, as shown in the figure at the left,
then GMm/a2 = [Mm/(M +
m)]a(dθ/dt)2, since v = a(dθ/dt) and the
centripetal acceleration is v2/a. This is the equation

Pure Mathematical Physics

 13

of motion for the relative motion, with r the distance

between M and m. The center of gravity of M and m
moves with constant velocity, not affected by the
relative motion. Thus, G(M + m) = a3(dθ/dt)2 =
4π2(a3/P2), where P is the period of the motion. This
is a derivation of Kepler's third law for a special
case, but it is a general result when a is the mean
distance of any orbit. The quantity G(M + m) is
usually called μ, so that r3/P2 = μ/4π2. If r is in
astronomical units (radius of the earth's orbit, about
150 x 106 km), and the period P is in years, then
G(M + m) = 4π2 = μ. The quantity GM is called the
Gaussian gravitational constant k2 if time is in days
instead of years, so we get approximately k2 =
4π2/(365.25)2, and k = 0.017202. A more precise
value is 0.01720209895. Untangling G from GM
means knowing the mass of the sun accurately,
which is impossible, so celestial mechanics uses k
instead of G. The value should be corrected for the
finite mass of the earth, but since the ratio of the
sun's to the earth's mass is 332,946 the correction is
small.
The polar equation of an ellipse, with the origin in a
focus, is r = p/(1 + e cos θ). The eccentricity of the
ellipse is e, and p is a length called the semi-latus-

Pure Mathematical Physics

 14

rectum. In connection with orbits, the angle θ is

called the true anomaly. In an elliptical orbit, 0 ≤ e <
1. When e = 0, we have a circular orbit r = p. The
same formula gives a parabola for e = 1, and a
hyperbola for e > 1. In an ellipse, the maximum and
minimum radii are p/(1 - e) and p/(1 + e),
respectively, and the corresponding points on the
orbit are called aphelion and perihelion,
respectively. The distance from perihelion to
aphelion is the major axis of the ellipse, twice the
mean distance a. Therefore, 2a = p/(1 - e) + p/(1 + e)
= 2p/(1 - e2), or p = a(1 - e2). The perihelion distance
is a(1 - e), the aphelion distance a(1 + e). The
distance from the center to a focus is a - a(1 - e) = ea
= c. This defines the eccentricity as e = c/a. Since the
ellipse is the locus of points the sum of whose
distances to the two foci is a constant, 2a, the length
of the radius from a focus to the end of the minor
axis is a. Since this is the hypotenuse of a right
triangle, a2 = c2 + b2, where b is the semiminor axis.
Hence b = a √(1 - e2). The reader can show that the
average value of r over one revolution is a. This is
about all we need to know about the geometry of the
ellipse for present purposes. More information can
be found in The Ellipse.

Pure Mathematical Physics

 15

If we know the orbit is an ellipse, and the areal

velocity is constant, we should be able to prove that
the force is central and varies as the inverse square
of the distance. On the other hand, if we know that
the force is central and inverse-square, then we
should be able to prove that the orbit is an ellipse
and the areal velocity is constant. Either of these
things can be done fairly easily. In the first case,
differentiate the expression for r with respect to
time, obtaining dr/dt = (e/p) sin θ 2A, where A is the
areal velocity (this was called dA/dt above; now A
represents the derivative). Then differentiate again to
find d2r/dt2 = (4A2e/pr2) cos θ. The tangential
acceleration is zero from the constancy of the areal
velocity. The radial acceleration is ar = (4A2e/pr2)
cos θ - 4A2/r3 = - 4A2/pr2, which shows that the force
is inverse square, and we have already shown that it
is central. In fact, GM = 4A2/p.
In the second case, the equations of motion in the
plane of the orbit are d2r/dt2 - r(dθ/dt)2 = -GM/r2 and
rd2θ/dt + 2(dr/dt)(dθ/dt) = 0. Now, (d/dt)(r2dθ/dt) =
2r(dr/dt)(dθ/dt) + r2d2θ/dt2 = 0, so the areal velocity
A = (1/2)r2(dθ/dt) = constant. Using this result, the
radial equation becomes d2r/dθ2 - (2/r)(dr/dθ)2 - r = -
r2GM/4A2. Making the substitution u = 1/r, this

Pure Mathematical Physics

 16

equation can be thrown into the form d2u/dθ2 + u =

GM/4A2. This second-order linear equation with
constant coefficients is easy to solve. We find u = 1/r
= GM/4A2 + C cos θ. This is the equation of a conic
section, as we have seen above. If p = 4A2/GM and e
= pC, we have r = p/(1 + e cos θ), which is what we
wanted to prove.
Now we can find the components of the velocity.
The transverse component is rdθ/dt = 2A/r = 2A(1 +
e cos θ)/p, where A is the areal velocity. The radial
component is dr/dt = 2A(e/p)sin θ. The square of the
velocity is the sum of the squares, or v2 = 4A2[(1 +
e2 + 2e cos θ]/p2] = 4A2[(1 - e2)/p2 + 2/pr] = -4A2/pa
+ 4A2/pr, from which the kinetic energy per unit
mass (i.e., setting m = 1 for simplicity) is: v2/2 = -
GM/2a + GM/r, since 4A2/p = GM. Rearranging,
v2/2 - GM/r = -GM/2a, which is to say, T + V = E,
with T = v2/2, V = -GM/r and E = -GM/2a. The total
energy is negative (i.e., the orbiting particle is
"trapped") and a function of the mean distance a
only. The eccentricity is determined by the angular
momentum h through e2 = 1 - h2/aGM. When h =
√aGM, the orbit is circular.

Pure Mathematical Physics

 17

The quantity a, called the "mean distance," is not the

average value of the radius vector r. It is the average
of the perihelion and aphelion distances, however.
The time average value of 1/r turns out to be 1/a,
which is its real definition. A proof will not be
supplied here until I find a simple one. This means
that the time average potential energy is -GM/a, and
the time average kinetic energy is -GM/2a + GM/a =
+GM/2a. For a circular orbit, the kinetic and
potential energies are constant. The time average
value of r is a(1 + e2/2), and the average of r over the
true anomaly is b, the minor axis of the orbit. The
average of 1/r over the true anomaly is 1/p.
We obtained the relation between the mean distance
a and the orbital period P for the special case of a
circular orbit. A more general proof is as follows.
The areal velocity dA/dt = h/2 = constant, so
integrating from t = 0 to t = P we find A = hP/2 =
πab. We also know that p = a(1 - e2) = h2/GM.
Therefore, P2 = 4π2a2b2/h2 = [4π2/GM]a3, since b2/p
= a. This is the desired general result. If we know the
period and mean distance of any orbit, we can
calculate the gravitational constant GM, which in the
general case is G(M + m).

Pure Mathematical Physics

 18

Coordinates and Orbits

The angular coordinates that

are used to describe directions
in space are shown in the figure
at the left. The equator and
ecliptic planes are shown. The
dihedral angle between these
two planes is ε = 23° 26'
21".412, or thereabouts. The
equatorial coordinates are the
right ascension α and the declination δ, while the
analogous ecliptic coordinates are the longitude λ
and the latitude β. The relations between the
coordinates are best found by considering the
rectangular components of the vector OP, considered
as of unit length, and then performing the rotation
about the x-axis that takes one system into the other.
We will not require these relations at the present
time.
The position of a body in the solar system can be
specified in terms of a gigantic rectangular
coordinate system centered on the sun. The x-axis is

Pure Mathematical Physics

 19

always taken as pointing to the vernal equinox, a

direction specified by the line of intersection of the
orbital plane of the earth (the ecliptic plane) with a
plane parallel to the equator of the earth, the
equatorial plane. Both contain the centre of the sun.
This is made interesting by the precession of the
normal to the equatorial plane around the pole of the
ecliptic. In the diagram, this is the movement of z'
(the pole of the equator) around z (the pole of the
ecliptic). Fortunately, this is not too rapid, taking
about 26,500 years for one revolution. The
equatorial plane has to be considered because our
earth-based measurements use it as a reference,
when directions are given in terms of right ascension
and declination. If we know a vector whose
components are referred to this plane, we can find
the right ascension and declination of its direction.
When the earth is located in its orbital plane on the
line of intersection of the ecliptic and equatorial
planes, the sun is seen in front of the vernal equinox
about March 21, climbing from below the equator to
above it. At this time, when we look at the sun
(careful!) we are looking in the direction of the x-
axis of the heliocentric coordinate system.

Pure Mathematical Physics

 20

The y-axis of the ecliptic coordinates is in the plane

of the ecliptic, at right angles to the x-axis, and in
the direction in which the planets move, in which a
screw would advance in the direction of the north
pole of the earth. The ecliptic z-axis is then
perpendicular to the x- and y-axes, and makes a
right-handed coordinate system with them. The
equatorial rectangular coordinates are defined
analogously. If we represent these coordinates by
primes, x and x' are the same, while x',y' is rotated
with respect to x,y by the dihedral angle between the
planes. The rotation from the z' axis to the z axis is a
right-handed rotation about the positive x,x' axis.
Longitudes are measured clockwise (eastwards) on
the ecliptic plane from the vernal equinox, from 0°
to 360°. The longitude of perihelion, ψ, is the angle
measured clockwise from the vernal equinox to the
radius passing through perihelion. (A different
symbol, a kind of ω, is used conventionally, but it is
not available here in HTML.) The longitude of
perihelion of the earth's orbit is currently about 103°.
This locates the major axis of the orbit and the
direction from which the true anomaly is measured.

Pure Mathematical Physics

 21

Other planets move on

orbital planes inclined to the
ecliptic by the inclination, i.
They intersect the ecliptic
on a line called the line of
nodes. The longitude of the
ascending node, Ω, is the
angle measured clockwise
from the vernal equinox to the line of nodes, on the
end where the orbit is climbing above the ecliptic
plane, called the ascending node. The node 180°
away is, not surprisingly, the descending node. Note
in the figure that the traditional symbol for the
ascending node is used, which we represent in text
by Ω, which looks similar. The symbol for the
descending node is inverted, with the two little loops
on the top. The longitude of perihelion in this case is
the longitude of the ascending node, plus the angle ω
in the orbital plane from the node to the perihelion.
This may be a little odd, but it does locate the
perihelion. That is, ψ = Ω + ω. These quantities are
illustrated in the diagram.
Now that the orbital plane and the direction of the
perihelion have been specified by the three orbital
elements i, Ω and ψ, we proceed to specify the shape

Pure Mathematical Physics

 22

of the orbit by its mean distance a and its

eccentricity e. The position of the body in its orbit at
a specified time is given by a sixth element, the
mean longitude at the epoch. The "epoch" is just the
reference time assumed, usually by its Julian date.
The mean distance a determines the orbital period P
through Kepler's third law. Therefore, only one of a
or P need be specified. The mean daily motion, 360°
divided by the orbital period, is often given in place
of P. Orbital elements for planets, asteroids and
comets can be found online at the Jet Propulsion
Laboratory website whose link is given in the
references.
The constancy of the areal
velocity is used to locate the
body in its orbit at a
specified time. This is called
Kepler's Problem, which
Kepler also solved. The
solution is illustrated in the
figure. A reference circle of
radius a is drawn around the
orbital ellipse. Point Q is
projected onto the auxilary circle by a perpendicular
to the major axis to point Q'. The line CQ' then

Pure Mathematical Physics

 23

defines the angle E at the centre of the circle, called

the eccentric anomaly. We can show that it is related
to the radius vector r and the true anomaly θ by r =
a(1 - e cos E) and tan (θ/2) = √[(1 + e)/(1 - e)] tan
E/2. The shaded sector PFQ is the area swept out by
the radius vector, and the area is proportional to the
time. The area FPQ' is proportional to the area FPQ,
since it is merely magnified vertically in the ratio a/b
= 1/√(1 - e2). Indeed, the total area of the orbit is
πab, while the total area of the circle is πa2.
Therefore, this larger area FPQ' increases uniformly
with time as well. Since the total area is πa2, the area
swept out is given by (πa2/P)t, where t is the time
since perihelion passage. This area can also be found
in terms of E, since it is the area of the whole sector
PCQ' less the area of the triangle CQ'F. The area of
the sector is πa2(E/2π) = (a2/2)E. The area of the
triangle is (1/2)(a sin E)(ae) = (a2/2)e sin E.
Therefore, we have E - e sin E = 2πt/P = M, called
Kepler's Equation. M is called the mean anomaly
and increases proportionally to time after perihelion
passage. This is really a beautiful result, allowing us
to find the position in orbit in terms of a uniformly
increasing quantity.

Pure Mathematical Physics

 24

Kepler's Equation can be easily solved by iteration.

If it is written E = M + e sin E, we find the iteration
equation En+1 = M + e sin En (e is the eccentricity
here, not the base of natural logarithms). Starting
with E0 = M, a few interations are enough to get
good precision, especially when e is small. The
iterations are very easy, and can be performed on a
pocket calculator. For e = 0.1, six iterations give 8
digits. The Newton-Raphson method can also be
used. This method finds the roots of the equation
f(x) = 0 by improving an initial guess x0 by the
formula x n+1 = x 0 - f(xn )/f'(x n ). A sketch will reveal
the reasoning behind this formula. In this case, f(E)
= M - E + e sin E = 0 is the function, and f'(E) = e
cos E - 1. Fewer iterations are required, but each
requires more work, than in the case of simple
iteration. In practice, an expansion in powers of e is
often used for small e, but the iteration methods are
no more trouble, especially for a computer.
Although the determination of position in orbit may
sound complicated, it is really quite straightforward,
and after you have done it a few times, it will be
easy. Let's find the position of the earth on 2003
April 2. The figures come from the Astronomical
Almanac for 2003, p. E4. On Jan 1, JD 245 2640.5,

Pure Mathematical Physics

 25

the mean longitude of the earth was 100.2440°, and

the longitude of perihelion was 103.0°. April 2 is 91
days later, so M = 100.2440 - 103.0 + 91 x 0.9856 =
86.93° = 1.5172 radians. A first approximation to E
is M + e sin M, and e = 0.0167, so the correction is
0.01668 radians, and E = 1.5339 radians or 87.89°.
There is no need to iterate again. The mean distance
can be taken as 1.00000, so r = 0.9994 and θ =
88.84°. The earth's longitude is 103.0 + 88.84 =
191.84°. On p. C8, the solar ephemeris gives
longitude 11.8506° and r = 0.9994. The solar
longitude is 180° less than the earth's, or 11.84°
according to our calculation. In both cases, we are
using the mean ecliptic of date, and came quite close
using Keplerian orbital elements instead of the more
elaborate calculations used to calculate the solar
ephemeris. Precise orbits require the consideration
of planetary perturbations, and this is a very difficult
subject indeed.
The orbit of the moon is an example, since it is
subject to perturbations by the earth's spheroidal
shape and the sun that make its orbit vary in a very
complicated way. In the case of the earth's orbit, we
really calculate the position of the earth-moon
barycentre. This point is 4671 km from the centre of

Pure Mathematical Physics

 26

the earth, about two-thirds of the way to the surface.

The mean distance of the moon is 3.844 x 105 km,
and its sidereal period is 27.321661 days, from
which G(M + m) can be calculated and used to relate
mean distance and period for other earth satellites
after correcting for the mass of the moon. The
moon's radius is 1738 km, and its mass 7.3483 x
1022 kg. m/M = 0.0123, which is not quite
negligible. The eccentricity of the moon's orbit is
about 0.0549, and its inclination 5.145396°. The
plane of the orbit and the position of perigee change
rather rapidly. The small eccentricity of the earth's
orbit suggests that the earth-moon system did not
suffer any external disturbance during its formation,
and in particular any collision which created the
moon. This was a one-time event, so speculations on
what happened can never be proved, since no similar
state is known, or is ever likely to become known.
Except for Mercury and Pluto, orbital inclinations
are less than 4°, and eccentricities less than 0.0936,
the value for Mars. Kepler was fortunate to have
picked Mars for study because of its relatively large
eccentricity, about twice those of Jupiter and Saturn.
Venus's orbit is practically circular, with e = 0.0067.
Mercury not only has a large inclination, 7.005°, but

Pure Mathematical Physics

 27

an eccentric orbit, with e = 0.2056. Only Pluto has a

greater inclination and eccentricity, but it is probably
a special case. Observations of Mercury and Pluto
were not available to Kepler, who had to work with
Tycho's naked-eye observations.
The angle made by r with the x-axis in our problem
was 191.84°, so the coordinates will be x = 0.9994
cos 191.84° = -0.97814, and y = 0.9994 x sin
191.84° = -0.20506, in astronomical units. In the
more general case of a planet in an inclined orbital
plane, the radius is first projected on the line of
nodes [a = r cos (ω + θ)] and on a perpendicular
(dip) line [b = r sin (ω + θ)]. Then z = b sin i, y = b
cos i cos Ω + a sin Ω, x = a cos Ω - b cos i sin Ω. It
is just a matter of projecting r on the coordinate
axes. The necessary formulas are given on page E4
of the almanac. There is an expansion for the true
anomaly directly in terms of the mean anomaly,
without going through the eccentric anomaly.
To find the equatorial rectangular components,
rotate about the x-axis through an angle of i. This
gives y' = y cos i - z sin i and z' = y sin i + z cos i. If
x', y' and z' are the components of a vector in
equatorial rectangular components, then the right

Pure Mathematical Physics

 28

ascension α is given by cos α = x'/√(x'2 + y'2), and

the declination δ by sin δ = z'/√(x'2 + y'2 + z'2). The
proper quadrants have to be determined. In this way,
the direction in which a planet is seen from the earth
can be found.

Pure Mathematical Physics

 29

Cometary Orbits

Comets can be divided into two classes, the short-

period comets, like Halley's Comet, with a period of
76 years, and long-period comets, like Hale-Bopp
(1995 O1), which do not return for thousands of
years, if at all. A period of 200 years is the
conventional dividing-line. Short-period comets
have elliptical orbits like the planets, except that the
eccentricity is larger, the inclination can take any
value, and the comets can move in a direct or
retrograde direction. Retrograde motion is usually
expressed by an inclination greater than 90°. Comets
were the first real test for Newton's theory, which
finally showed that they were normal members of
the solar system, not mysterious atmospheric
happenings, as had always been supposed. Positions
of short-period comets are calculated in the same
way as planetary positions, and the orbital elements
are presented the same way.

Pure Mathematical Physics

 30

Long-period comets come from the periphery of the

solar system, where they wander in the hypothetical
Oort Cloud of cometary debris, normally water and
carbon dioxide ice and dust. Their total energy is
about zero, so the eccentricity of their orbits when
they make an excursion towards the sun is about 1.
That is, their orbits are parabolic. Short-period
comets are those that have suffered an energy-losing
collison with a planet (usually Jupiter) and have
dropped into an elliptical orbit. There must also have
been comets that have gained energy, entered
hyperbolic orbits, and were ejected from the solar
system. Few, if any, comets have eccentricities
significantly greater than 1, which would indicate
that they were encountered by the solar system in its
path through space, and are not part of the family.
The parabola is a considerably simpler orbit than the
ellipse. Its polar formula is r = p/(1 + cos θ) = (p/2)
sec2 (θ/2). The areal velocity A = √(GMp/4). Setting
x = tan (θ/2), direct integration gives Kepler's
Equation as x + x3/3 = √(GM/2)(2/p)3/2(t - T). The
one real root of this equation gives the position in
orbit. There are many ways to solve a cubic
equation: refer to mathematics handbooks (see
References), algebra texts or mathematics programs.

Pure Mathematical Physics

 31

The time to move from perihelion to the end of the

latus rectum is t1 = (1/2)√(p3/GM). 2t1 is a good
measure of the time that the comet will spend near
the sun.
The orbital elements include the inclination i,
longitude of the ascending node Ω, and argument of
perihelion ω. The perihelion distance p/2 is given,
and the JD of perihelion passage. For Hale-Bopp, the
perihelion distance was 0.91399384 and the time of
perihelion passage was JD 245 0539.60742. The
inclination was 89.42064850°, longitude of
ascending node 282.47215310°, and argument of
perihelion 130.59561740. From these elements, it
was possible to find the position of the comet as
seen from the earth. Cometary orbital elements are
available on the internet soon after the discovery of a
comet. The Jet Propulsion Laboratory supplied the
elements for Hale-Bopp shortly after it was
discovered.
Halley's comet was the first periodic comet to be
recognized. Edmund Halley (1656-1742), friend of
Newton's and later Astronomer Royal, suspected that
the comets of 1531 and 1607 were the same as the
comet of 1682, since they had similar orbits. He

Pure Mathematical Physics

 32

predicted the return of the comet in 1758, which

duly occurred. Halley had published Newton's
Principia at his own expense in 1678, and made an
extensive study of comets using the new theory.
Halley's comet is the brightest of the short-period
comets; most are quite dim, and can even have orbits
of small eccentricity like asteroids.
To show how to use the orbital elements given by
JPL, let's look at them for Halley's comet. The time
of perihelion passage, Tp is shown as
19860209.45895, which I interpret to mean 1986
February 9, at 0.45895 part of the day, or 11h 0m
53s UT. I should have preferred the Julian Day,
which would be unambiguous. The "epoch" is the
date to which the elements apply, which is given as
46480. This is a modified JD. The JD can be
obtained by adding 24400000.5, or JD 2446480.5,
which is 1986 February 18. The size and shape of
the orbit is specified by the perihelion distance q =
0.58710374, in AU, and e = 0.96727724. The
eccentricity is nearly, but not quite, unity. Near
perihelion, the orbit will be indistinguishable from a
parabola, so the parabolic formulas can be used for
predicting its position while it is near the sun and
visible. The orientation of the orbital plane is

Pure Mathematical Physics

 33

specified by the longitude of the ascending node,

58.86004°, and the inclination i = 162.24220°. When
the numbers are located, look at the symbols used in
the website, which may not be the usual ones
because of font limitations. Since i > 90°, the motion
of the comet in its orbit appears to be retrograde,
opposite to the direction of motion of the planets, if
the inclination is taken to be 180° - 162.24220° =
17.75780°. Finally, the orientation of the orbit in its
plane is given by the argument of perihelion,
111.8656°, measured from the ascending node in the
direction of movement. It is not easy to comprehend
the orbital position from the bare numbers; a
drawing, made with some care, will show exactly
what is going on.
As the perihelion distance is q = a(1 - e), it is a
simple matter to determine a. In fact, a = q/(1 - e) =
0.58710374/0.03272276 = 17.941755 AU. From
this, the orbital period P can be found. Since P2/a3 =
1 if P is in years (earth = 1) and a is in AU (earth =
1), P = 75.99716 years, very close to 76.0 years. The
aphelion distance is 2(17.941755) - 0.58710374 =
35.296 AU. This is outside the orbit of Neptune (a =
30) but not as far as Pluto (a = 39.5). In accordance
with the law of areas, Halley's comet spends most of

Pure Mathematical Physics

 34

its time drifting in this dark region, periodically

darting in to the sun to have a little more of itself
boiled off into space.
Although orbital elements are given to high
precision, the accuracy is not necessarily as high.
Especially for comets, they can be changed by
planetary perturbations, sometimes by large
amounts. The ejection of gases from comets also
leads to reaction forces that can affect the orbit.
Also, the reference coordinates may change with
time, because of precession of the equinoxes and
other effects, and this must be taken into account in
accurate calculations. These changes are not changes
in the orbit, of course.

Pure Mathematical Physics

 35

Determination of Orbits

Finding the orbital elements from observations, and

predicting the changes in orbital elements due to
perturbations, are two of the most important
problems in celestial mechanics, and have received
close attention from Newton's time onwards. We
cannot give any reasonable account of this work
here, but we can show how orbital elements come
from observed motions by a graphical analysis that
is very instructive, though of little practical use
where high precision is required. The reader with
drawing supplies is encouraged to follow along.
We suppose known the position r of a body M from
the sun and its velocity v relative to the sun. These
are six parameters that will serve to determine the
six orbial elements. The orbital plane is defined by
the plane of r and v, and can be represented in two
views by means of orthographic projection
(descriptive geometry; see Monge's Procedure). The
line of intersection of the orbital plane with the

Pure Mathematical Physics

 36

ecliptic plane can then be found, and the dihedral

angle between the two planes. Since the line of
intersection will be the line of nodes, we have now
found two of the orbital elements, Ω and i, which
determine the plane of the orbit.
The next step is to calculate
the mean distance a, using the
energy relation v2 = GM(1/r -
1/2a). From a we can find the
orbital period by P2/a3 =
4π2/GM. The radius vector
SM and the line from M to the
other focus of the ellipse make equal angles with the
tangent to the ellipse at M, which is in the direction
of v. In a view showing the orbital plane in true size,
Lay off SM and MQ, where MQ is an arbitrary line
in the direction of the empty focus. The empty focus
S' is located by laying off a distance 2a - r along
MQ', since the sum of the distances from the foci to
M is 2a. Now we can draw the line SS', find its
center at O, and then the locations of perihelion and
aphelion. The eccentricity of the orbit is e = OS/a =
OS'/a. The orbit can now be drawn, and the angle
from the line of nodes to the perihelion, the
argument of perihelion, ω, can be measured, as well

Pure Mathematical Physics

 37

as the true anomaly θ, which is the angle PSM. We

now have the five elements Ω, i, ω, e, and a.
The remaining element is the time of perihelion
passage, which can be found from the true anomaly.
First, we find the eccentric anomaly by tan(E/2) =
[(1-e)/(1+e)]1/2 tan(θ/2), and from it the mean
anomaly M = E - e sin E. The mean anomaly is M =
2π(t - T)/P, where T is the time of perihelion
passage, and t is the time of observation. That is, the
body passed through perihelion at a time t - T =
MP/2π earlier. We have now determined all six
orbital elements from the six components of the
initial position and velocity. This was easy because
the distance r was one of the given quantities. In
general, orbits must be determined by observations
of the angular position from earth, not directly in
terms of distance and velocity. The position at three
different times is sufficient to determine the orbit in
this case, but the analysis is more difficult than what
we have done.
Orbital elements vary only slowly in the solar
system, except when there are near encounters,
usually of light bodies with massive planets or
satellites, that affect the orbits of the light bodies,

Pure Mathematical Physics

 38

but not of the heavier ones. If we know the orbit at

one instant, then we can predict the velocity and
position at a later time from the orbit. We can also
integrate the actual changes in position and velocity
at the later time, taking into account forces exerted
by bodies other than the sun. The difference will be
reflected in the orbital elements, which will slowly
change. This is a useful way to take perturbations
into account, and is widely used. The maximum
errors in using the Keplerian orbits are on the order
of 5" to 30" for the inner planets, somewhat more for
Jupiter and Saturn. Neptune was discovered in 1846
as a result of its perturbations of the motion of
Uranus. This was a remarkable demonstration of the
accuracy of Newtonian mechanics.
The earth's mean distance is decreasing by 5 x 10-8
AU per century, about 75 m per year. The
eccentricity is decreasing by 3.804 x 10-5 per
century. The inclination is decreasing by 46.94" per
century. These are average rates, and are simply the
current values of changes that may be periodic.
Nevertheless, they show how gradual any changes
are, and allow the calculation of the earth's position
with reasonable accuracy for thousands of years
either side of the present. The longitude of the node,

Pure Mathematical Physics

 39

however, moves at -3.038" per year because of the

precession of the equinoxes, which moves the
reference point. The argument of perihelion is
increasing at 11.9828" per year. It is typical for an
orbit to show a relatively rapid change in longitude
of the node and motion of the perihelion; the moon
is an excellent example, and also earth satellites,
which are perturbed by the equatorial bulge of the
earth and tidal forces exerted by the sun and moon,
as well as by atmospheric drag if they are in low
orbits.

Pure Mathematical Physics

 40

Earth Satellites

The mass of an earth satellite is infinitesimally small

compared to the mass of the earth, so the centre of
the earth may be considered as a fixed point. A
reference frame fixed to the earth's center is not an
inertial frame, due to revolution about the sun and
the moon's motion. The earth-moon barycentre is a
"more inertial" reference point, but in any case the
difference from an inertial frame is negligible. There
is a change in terminology, as perihelion becomes
perigee and aphelion becomes apogee. This useless
distinction is regrettable, and can be carried to
excess (perijove, etc.). It would be good if there
were terms usable in any orbit.
GM for the earth is 3.98600440 x 1014 m3/s2.
Therefore, a3/P2 = 10097 if a is in km and P is in
seconds. For P = 86164.0989 s (sidereal rotation
period of the earth), a = 42,164,172 m, or an altitude
above the surface of 22,232 miles. A satellite with
this mean distance will return to the same point in
the sky each day. If the inclination and the
eccentricity are zero, the satellite will appear to hang

Pure Mathematical Physics

 41

motionless in the sky. Such an orbit is called

geostationary, obviously of great value as a
communications relay location. An area of very
choice space real estate is thereby created, for which
there is considerable demand. Satellites occupy fixed
stations around the equator. They remain in a
constant direction from an observing point on the
earth's surface, so an antenna can be permanently
pointed at them without having to be continually
redirected, which would not be easy nor cheap. Their
positions may be adjusted slightly to keep the
satellites on station by on-board thrusters during
their effective life, and the thrusters are used to
remove them from orbit when their fuel runs out,
since space is at a premium in the geostationary belt.
From a geostationary orbit, the earth subtends an
angle of 17.4°, the same as a 12"-diameter globe at a
viewing distance of 40", so it is quite reasonable to
take photos of the whole visible hemisphere at once.
Examples are available at the Dundee University
website in the References. The poles are not visible,
since the view extends only to latitude 81.3°.
The first artificial earth satellite was Sputnik I,
launched on 4 October 1957 and weighing 83.6 kg.

Pure Mathematical Physics

 42

Its orbital period was 96.15 min, so its mean

distance was 6953 km. The eccentricity was 0.0517,
giving a perigee altitude of 228 km and an apogee
altitude of 947 km. The inclination of the orbit was
64.26°, so it could easily be seen in every part of the
earth when the sun shone on it. Sputnik remained in
orbit until January 1958, making 1350 revolutions.
Because of the low perigee, it was considerably
affected by atmospheric drag, its final period being
near 90 min. The period of a satellite that just grazes
the surface of the earth would be 84.48 min, and its
velocity 7906 m/s.
The Global Positioning System (GPS) uses satellites
that continually radiate very accurate timing
information, corrected for relativistic time dilatation
due to the satellite's speed, so that their distances
from an observer can be found to within a few
centimetres from the time differences. The orbits
must be known to a similar accuracy, so the
effectiveness of this system is evidence of the
correctness of the dynamics. Orbit information is
broadcast along with timing information. The carrier
frequencies are 1.57542 and 1.22760 GHz. The
complete GPS system was deployed in 1993,
consisting of 21 operational satellites and three

Pure Mathematical Physics

 43

spares, on circular orbits inclined at 55° and with a

period of one-half a sidereal day (a sidereal day is
about 23h 56m), so the satellites are seen to rise and
set about 4 minutes earlier each day, and appear
twice a day. Four satellites in good postions are
intended to be visible from any point at any time,
giving a good intersection, and eliminating the
necessity for accurate calibration of the receiver
clock (the offset of the receiver clock is one of the
four unknowns that can be determined). The
altitudes, about 20,200 km (a = 26,560 km), are high
enough to make atmospheric drag negligible. There
are six orbital planes, with four satellites spaced
equally in each. The GPS is by no means the only
satellite navigation system that has been developed,
but it has become the most used, and is replacing the
others.
Earth satellites are affected by many small forces in
addition to the main inverse-square force directed
toward the centre of the earth, and these
perturbations cause the orbital elements to change
slowly with time. In the case of earth satellites, the
mean distance a, the eccentricity e and the
inclination i do not change on the average (they may
fluctuate slightly over the short term). The longitude

Pure Mathematical Physics

 44

of the line of nodes, Ω, the argument of perigee, ω,

and the rate of change of mean anomaly change
steadily with time. The drag of the atmosphere is
negligible for satellites that do not come lower than
about 1000 km at perigee, which includes most
practical satellites. For lower satellites, atmospheric
drag causes a loss of energy at perigee that pulls in
the apogee position, making the orbit less eccentric
and decreasing the orbital period. Paradoxically, the
drag speeds up the satellite! The pulls of the moon
and sun, tidal effects cause orbital changes. The
direct tidal effect is due to the force exerted directly
on the satellite, while the indirect tidal effect is due
to the changes in mass distribution of the earth
caused by tidal motions. Pressure of solar radiation
and solar wind is another disturbing effect. The most
important perturbing force, however, is that exerted
by the equatorial bulge of the earth.
The gravitational potential of the earth can be
expressed approximately as V = GM/r [1 -
(aE/r)2J2 P2 (sin φ)], where aE is the equatorial radius
of the earth, 6 378 137 m, J2 = 0.001082630, and
P2 (x) is the Legendre polynomial of second degree,
(1/2)(3x2 - 1). The argument of this polynomial is
usually seen as cos θ, but here we use the latitude φ

Pure Mathematical Physics

 45

= 90° - θ. This is the first part of an expansion in

spherical harmonics, which can even allow for a
lumpy earth that is not symmetric about the polar
axis. The correction to the 1/r field shown here is a
zonal harmonic of order 2, related to the flattening f
= a/(a - b) = 1/298 of the earth. However, the
gravitational potential depends on the distribution of
mass in the earth, so it cannot be expressed simply in
terms of f. The radial component of the force is -
∂V/∂r, and the tangential part is (1/r)∂V/∂θ. The
motion of earth satellites has led to a much better
knowledge of the earth's gravitational potential.
If a satellite has a mean motion n = 2π/P, a mean
distance a and an inclination i, the movement of the
line of nodes is dΩ/dt = -(3n/2)(aE/a)2[cos i/(1 -
e2)2]J2 , and the change of the argument of perigee is
dω/dt = +(3n/4)(aE/a)2[(5cos2i - 1) /(1 - e2)2]J2 . The
quantities are in MKS units, and the results are in
radians per second. For the moon, the line of nodes
rotates once in 18.6 years. For a typical GPS
satellite, the line of nodes moves about -0.03° per
year, the perigee about 0.01° per year. All the other
perturbations are much smaller in amount, but must
be taken into consideration in accurate work.
Relativistic perturbing accelerations are inversely

Pure Mathematical Physics

 46

proportional to the fourth power of the distance r.

Their magnitudes for GPS satellites are about 3 x 10-
10
m/s2, less than other small perturbations, so they
have no practical effect, although they have been
considered.
The International Space Station
(ISS), shown at the right, is a
famous earth satellite. The image
is from the NASA website. The first two modules,
Zarya and Unity, were assembled in orbit in 1998,
and now there are 14. It is 73 m across the solar
arrays, 44.5 m long and 27.5 m high, with 425 m3 of
habitable space. Power comes mainly from the 892
m2 of solar arrays. There are thrusters to adjust the
orbit when necessary. The total mass is about 179
metric tons. Although this is called "space travel" by
NASA, it is not even outside the earth's atmosphere!
Approximate elements of the ISS orbit at 13.49Z, 6
April 2003, are i = 51.6°, Ω 19.99°, e = 0.00083, ω =
52.91°, M = 307.29°, and mean motion 15.59579861
rev/day, or 0.0649825°/s. ("Z" is another designation
for UT.) Exact orbits, predicted about 10 days in
advance, can be found at the NASA link in the
References. The period is P = 5539.95 s (1.54 hr),

Pure Mathematical Physics

 47

and the mean distance is 6767009 m. This is a low,

almost circular orbit at an altitude of about 400 km.
Because the orbit is low, the ISS is seldom
illuminated by the sun at night, so it is not frequently
seen, though a frequent, large and prominent object.
The longitude (RA) of the ascending node changes
rapidly, by about -5.00° per day, as does the
argument of perihelion, by -3.64° per day.
Perturbation by the oblateness of the earth is the
reason for most of this change. The inclination,
eccentricity and mean distance do not change
rapidly. The orbit is just above the maximum
ionization in the F layer of the ionosphere, and so is
affected by atmospheric drag. The density of the
atmosphere at this height is about ρ = 9 x 10-12
kg/m3. The drag is given by F = CρV2A/2, where C
is the drag coefficient, V the satellite velocity, and A
the effective projected area. NASA gives A = 344
m2 and C = 2.36. V = 2πa/P = 6541 m/s (assuming a
circular orbit), so F = 0.16 N, which will produce an
acceleration of 8.9 x 10-7 m/s2. This drag causes the
mean distance and period to decrease. NASA notes
that the "decay" is 4.11 x 10-4 rev/day2. In 100 days,
the mean motion will increase by 0.0411 rev/day, or
P = 5527 s, a decrease of 13 s, or 0.23%. If

Pure Mathematical Physics

 48

uncorrected, the satellite would spiral inward at an

increasing rate, eventually burning up
catastrophically.
It may be interesting to find out where you should
look for a satellite in a known orbit. The procedure
will be illustrated here without taking all the
refinements into account that are necessary for, say,
GPS positioning. The idea is to find the rectangular
coordinates of the satellite in an approximately
inertial system with its origin at the centre of the
earth, and then to find the rectangular coordinates of
the point of observation in the same system. The
differences of the coordinates then give a vector
from the point of observation to the satellite. The
motion of a satellite as seen from a fixed location on
earth may be very complex, because of the
interaction of the two motions involved, the
revolution of the satellite and the rotation of the
earth. Only when these are approximately equal is
the situation more or less simple.
Let's suppose the satellite is in a circular orbit with a
= 26 500 000 m and i = 55°, like a GPS satellite. Let
the longitude of the ascending node be 40°. Since
the orbit is circular, there is no perigee point, so we

Pure Mathematical Physics

 49

measure the true anomaly from the line of nodes in

the orbital plane. Let the true anomaly at the time we
are considering be 70°. The z-axis is taken along the
axis of rotation of the earth, and the x-axis in the
direction of the vernal equinox. The y-axis then
makes a right-handed coordinate system. This is just
like the case of planetary motion, except that the role
of the ecliptic plane is played by the equatorial
plane. The orientation of these axes remains fixed in
space as the earth revolves about the sun.
First, resolve the radius vector along and
perpendicular to the line of nodes. The components
are a cos θ = 9 084 055.0 m and a sin θ = 24 958
236.0. Then, z = a cos θ sin i = 20 444 590.0 m, x =
a cos θ cos Ω - a sin θ sin Ω = -9 084 055.0 m, and y
= a cos θ sin Ω + a sin θ cos Ω = 24 958 236.0 m.
This is the instantaneous position of the satellite in
the inertial system, and it could obviously be found
for any time equally easily.
Now for the point of
observation. Let's use my own
location, which is longitude λ =
-104.92583°, latitude φ =
39.72694 and height above the

Pure Mathematical Physics

 50

ellipsoid (or mean sea level) h = 1633.7 m. Your

own coordinates, if you do not know them already,
can be found from a USGS topographic map. A
meridional section of the earth is shown at the right,
where φ is the geographic or spheroidal latitude. The
dimensions of the earth are a = 6378136 m and b =
6356752 m. There is a number of reference
spheroids, any one of which will work satisfactorily.
The eccentricity of the elliptical cross-section of the
spheroid can then be found to be e = 0.006694167.
The meridional radius of curvature at latitude φ is N
= a[cos2φ + (1 - e2)sin2φ]-1/2, or 6 378 194.38 m at
my location. The length of N below the equatorial
plane is e2N = 285.8 m.
We now choose an earth-fixed rectangular
coordinate system with the z'-axis along the
rotational axis of the earth, the x'-axis in the
meridian of Greenwich, and the y'-axis making a
right-handed system. The rectangular coordinates of
my location are then x' = (N + h)cos φ cos λ = -1 263
816.21 m, y' = (N + h)cos φ sin λ = -4 741 167.76 m,
z' = (N + h - e2N)sin φ = 4 077 353.73 m. These
coordinates do not change as the earth rotates. The
distance from the centre is 6 379 711.32 m.

Pure Mathematical Physics

 51

The relation between the inertial

and earth-fixed coordinate
systems is shown at the left. The
earth-fixed system rotates
steadily in the positive direction
about the common z and z' axes,
and the angle of rotation Θ is the local apparent
sidereal time at Greenwich, which is easy to find
from the Astronomical Almanac. The rotational
period of the earth is 23h 56m 04.09890369732s.
When the time is stated to such precision, there are
many small corrections to be considered, but we
shall ignore them here so the principle is not lost.
We shall assume that the inertial and earth-fixed
coordinates are related simply by the rotation
through an angle Θ, not correcting for polar motion,
precession of the equinoxes and nutation. These are
simply rotation matrices applied in addition to the
one for the main rotation. For accurate work,
however, they cannot be ignored. Since sideral time
goes from zero to 24 hours in one rotation, we must
multiply an ordinary time interval by 24/23.9344719
or 1.00273781195 to find the equivalent sidereal
time interval. A sidereal hour is simply equivalent to

Pure Mathematical Physics

 52

15° of angle; it is not an "hour" in the usual sense of

time.
The transformation equations are the familiar x = x'
cos Θ - y' sin Θ and y = x' sin Θ + y' cos Θ. All you
have to do is check that the signs of the sines are
correct. If we put in coordinates in the earth-fixed
system (x',y',z') we will get out coordinates in the
inertial system (x,y,z) which we can compare with
the satellite coordinates. Let's pick 10.00 am MST
on 2003 January 1. Since my time zone is +7, the
UT is 17.00h January 1. From the Astronomical
Almanac, page B8, the sidereal time at 0h UT was
6.61651428 hours. The sidereal time elapsed is then
1.00273781195 x 17.00 = 17.04654280 hours.
Adding the two, the Greenwich LAST at my 10.00
am will be 23.66305708 hours, or 354.945856°. If it
seems easier, this can also be expressed as -
5.0541438°. Now I can find my inertial coordinates
at this instant. The results are x = -1 676 585.46 m, y
= -4 611 395.05, z = 4 077 353.73 m. The square
root of the sum of the squares is 6 379 711.32, so the
arithmetic checks.
Now we can subtract the coordinates of the
observation point from the coordinates of the

Pure Mathematical Physics

 53

satellite to find the relative vector (X,Y,Z). I find X

= -7 407 469.54 m, Y = 29 569 631.05 m, Z = 16
367 236.27 m. The direction of this vector can now
be expressed in terms of right ascension and
declination, since it is a vector in the inertial system
referred to the vernal equinox. The distance from the
observer to the satellite is 34 599 423.53 m. The
projection on the equatorial plane is 30 483 334.55
m, so that the declination δ = 28.2324° and the right
ascension is 75.9363° or 5h 3m 45s. At the date and
time specified, this direction is beneath the earth. In
fact, the antipodeal point of right ascension 18h and
declination -29° is in the southern sky at an altitude
of about 20° above the horizon.
Although these calculations are tedious and subject
to error, even when done with an electronic
calculator, a computer program can do them with
ease, speed and correctness, and even the more
precise calculations are not much bother. It would
not be difficult to write a program to determine the
visibility of all of the 24 GPS satellites, and the
directions in which they are to be seen.

Pure Mathematical Physics

 54

Other Orbit Lore

The total energy E of a body of

unit mass moving under the
influence of a potential energy
per unit mass V(r) is the kinetic
energy T = v2/2 plus V(r). In
polar coordinates, the
components of velocity are
dr/dt and r dθ/dt, so the total energy becomes E =
(dr/dt)2/2 + r2(dθ/dt)2/2 + V(r). The constancy of
angular momentum for a central force gives r(dθ/dt)2
= h = constant. Substituting for dθ/dt, we find E =
v2/2 + h2/2r2 + V(r). This is the energy equation for a
particle of unit mass moving in one dimension with
the effective potential h2/2r2 + V(r). The additional
term is called the centrifugal potential. In the case of
orbital motion under gravity, V(r) = -GM/r. The
form of the resulting potential is shown in the figure.
For E = 0, a particle coming in from infinity
accelerates steadily until r = h2/GM, then is rapidly
decelerated, coming to rest at r = h2/2GM, and

Pure Mathematical Physics

 55

afterwards retracing its path in reverse. Of course,

the particle swings round the centre of force in this
motion. For E = -(GM)2/2h2, r has the constant value
that gives the minimum of the curve, and the orbit is
circular. For intermediate energies, r oscillates
between a minimum and a maximum value,
corresponding to perhelion and aphelion.
Consider the vector B = v x h + (GM/r)r. The
angular momentum per unit mass, h, is here taken to
be a vector perpendicular to the orbital plane. The
cross product then lies in the orbital plane, and so
then does B. Taking its time derivative, we have
(dv/dt) x h - (GM/r3)r(v·r) + (GM/r)v. Now dv/dt =
-(GM/r3)r, and h = r x v, so the first term becomes -
(GM/r3)[r(r·v) - r2v]. All the terms cancel, so dB/dt
= 0, and B is a constant of the motion. This is
extraordinary, since we already have one constant of
the motion, the total energy E, and this is usually the
only one that can be found. For an inverse-square
force, we have found a second one. It is known that
this happens when we have an additional symmetry
in the problem--here it happens to be symmetry
under rotation in a four-dimensional space, but there
is no room to explain this here. We'll just be happy
with the result.

Pure Mathematical Physics

 56

Since B is a constant, it can be evaluated at any point

of the motion, and it is easiest to choose perihelion.
At this point, v x h = vp2rp , so B is a vector parallel
to the line between the focus and perihelion. We will
find the same vector if we evaluate B at any point of
the motion. Its constancy is a good proof of an exact
inverse-square force.
We have seen that with an attractive inverse-square
force, closed orbits exist for negative energies. For
zero or positive energies, the orbits are open. We
have considered briefly the zero-energy parabolic
orbits that divide the two classes, but not the
hyperbolic orbits, since they are not of much use in
the solar system. For a repulsive inverse-square
force, there are no closed orbits of negative energy at
all, only the hyperbolic orbits of positive energy. On
these orbits, the body comes in from infinity along a
path asymptotic to a line of a certain direction. It
then interacts with the other body, in what is called a
"collision," though the bodies do not come into
contact, and leaves to go to infinity asymptotic to a
different direction. The angle between the
asymptotes is a function of the impact parameter,
which is the distance of closest approach if the
bodies did not exert forces on each other. This is the

Pure Mathematical Physics

 57

problem of Coulomb scattering in atomic and

nuclear physics. Rutherford and his students used
classical mechanics to analyze the scattering of
alpha particles (Z = 2, A = 4) by atomic nuclei to
verify the existence of nuclei. This is an area
different enough from celestial mechanics to leave it
for another article.

Pure Mathematical Physics

 58

References

P. van de Kamp, Elements of Astromechanics (San

Francisco: W. H. Freeman, 1964).
F. R. Moulton, An Introduction to Celestial
Mechanics, 2nd ed. (New York: MacMillan, 1914).
Y. Ryabov, An Elementary Survey of Celestial
Mechanics (New York: Dover, 1961).
P. K. Seidelmann, ed., Explanatory Supplement to
the Astronomical Almanac (Mill Valley, CA:
University Science Books, 1992). Chapter 1 is
especially recommended.
The Astronomical Almanac for the Year 2003
(Washington, DC: U.S. Government Printing Office,
2001). pp. B8, C8 and E3-E5. Orbital elements for
the minor planets (asteroids) are also included.
Bronshtein and Semendyayev, A Guide Book to
Mathematics (New York: Springer, 1973). p. 161
(solution of cubic equation).

Pure Mathematical Physics

 59

I. Ridpath, ed., Norton's 2000.0 Star Atlas and

Reference Handbook, 18th ed. (New York: John
Wiley & Sons, 1989). pp. 114-117. There is a table
of cometary orbital elements.
L. D. Landau and E. M. Lifshitz, Mechanics
(Oxford: Pergamon Press, 1960). This is an excellent
concise text in classical mechanics, looking towards
applications in quantum mechanics and atomic
physics.
Orbital elements for planets, asteroids and comets
can be found at JPL.
Orbital elements for the ISS can be found at ISS
Orbit. The website contains information on NASA
missions.
An excellent source of earth satellite orbits
is Dundee University Satellite Station. Free
registration opens up excellent things, like IR and
VIS images from GEOS and other satellites, as well
as orbital elements of a number of satellites.
B. Hofmann-Wellenhof, H. Lichtenegger and J.
Collins, GPS Theory and Practice (Wien: Springer-
Verlag, 1993).

Pure Mathematical Physics

Pure Mathematical Physics