Home

Page 1 Elliptical Orbits

Page 2 Time Dependency

Page 3 Kepler’s Equation

Elliptical Planetary Orbits Derived from Newtonian Mechanics

The notion that planets revolve about the sun in elliptical orbits was found empirically by Johannes Kepler during his little “war with Mars.” Later, the fundamental principles elucidated by Isaac Newton provided a firm theoretical basis for understanding such orbits. Newton’s principles can produce quite general solutions to this idealized two-body problem of course, showing that planets may follow any of a number of trajectories described by conic sections. Of interest here, however, is only in showing that an expression for one particular orbit - the ellipse - falls right out of the mathematics.

In general, a straightforward technique for deriving an object’s trajectory is to begin with Newton’s second law of motion, F = ma. Substituting an expression for the force into the left-hand-side and writing the acceleration as the second time derivative of position results in a set of equations which can be integrated in order to find the object’s position as a function of time. Integrating twice may be a bit difficult for some trajectories.

Since the time dependence is not really the desired solution at the moment, the problem can be simplified by removing any time dependency and solving for the planet’s path in ordinary two-dimensional space. This solution displays the elliptical nature of the orbit. Properly speaking, though, for a planet with a moon or moons, it is the center-of-mass, or, barycenter of the planet/moon system which traces the elliptical orbit about the sun.

Newton’s expression for the gravitational force, F, is written in terms of the planet’s radial position, r:

where m₁ is the mass of the sun,

m₂ is the mass of the planet,

r is the radial distance between the sun and planet,

G is the universal gravitational constant,

is the unit vector in the radial direction.

The radial vector is

and the force can be written explicitly as a function of r,

Note that the unit vector in the radial direction is

Note on Center of Mass of Sun/Planet System

If the origin of our coordinate system is as shown in this figure:

then the kinetic energy of the system, KE, is the sum of the kinetic energy of motion of the center of mass plus the kinetic energy about the center of mass, given as

where R is the position vector of the center of mass and r₁ and r₂ describe the positions of the bodies relative to the center of mass. We have employed the dot notation for the first derivative with respect to time. These two position vectors are related as

If the origin of our coordinate system is at the center of the sun as described above in the first figure, then solving for the motion of the planet presumes that the sun is fixed in place; that is, as though there were no rotation of both bodies about a common center of mass. This implies that the sun has no kinetic energy within the solar system. This may be approximately true in our solar system, where the largest planet is small (but not necessarily insignificant considering modern astronomical measurement techniques) compared to the sun.

If instead we wish to place the origin at the center of mass of the Sun/Planet system, we would write the position vectors as

and

With these, the kinetic energy about the center of mass can be reduced to

Since we have no interest here in the motion of the solar system within the galaxy we may ignore altogether R and its time derivative. This process of ignoration is simply a matter of dropping the first term in the expression for kinetic energy.

Thus the system is reduced to an equivalent one-body problem where the planet acts as though it has an effective, or “reduced” mass

and the sun acts as though it is a body of mass

Newton’s universal law of gravitational attraction for a body of mass M and a body of mass μ results in the same force as given above. This equivalent one body problem is as though there were a single mass, M, at the center of rotation and a body of mass μ in orbit.

In our solar system, we may wish to ignore the reduced mass if we consider the planetary masses small compared to the sun's mass. For example, the mass of Jupiter is 1898.3x10²⁴kg and the sun’s mass is 1.9855x10³⁰kg. For this system, the reduced mass is 1893.2x10²⁴kg: a reduction of about 0.3%, and M = 1.9874x10³⁰ kg; so we can see here that using the reduced mass in the solar system is not really necessary.

There are other systems in which the two masses are comparable, such as binary star or planet/moon orbits. The analysis presented here is general and can be applied to these systems. For our purpose we will use terminology of the sun/planet system: aphelion and perihelion, but the general terms are apoapsis and periapsis.

Two-dimensional nature of the orbit

That the orbital path of a planet about the sun is two-dimensional has been asserted without any proof. This assertion should be shown to be true, perhaps, before moving on to the derivation of elliptical orbits. This two-dimensional nature of the orbits can be understood with a physical argument, as follows: the gravitational force vector is always anti-parallel to the position vector; at any given moment the planet’s position vector and velocity vector form a plane (except in the particular case of free-fall in which they are co-linear). Since the force vector is also in this plane there is no component of force perpendicular to the plane, and thus, there is no force to drive the planet outside the plane. This is a consequence of the centripetal (“center-directed”) nature of the gravitational force.

The idea of planar orbits can easily be shown mathematically. Define the moment which was just described as occurring at time t₀, and arrange the coordinate system such that the x-y plane aligns with the r-v plane. Since the component of force in the z-direction is zero, we can write this component in terms of acceleration, a_z, as

Dividing through by the mass of the planet and writing an indefinite integral over time, we obtain

This constant of integration is the z-component of speed at time t₀, which is zero at t = t₀ because at that moment the velocity vector is in the x-y plane; we are therefore left with v_z(t) = 0. Now we integrate a second time,

Here, the constant of integration is z(t₀) = 0. This result means that the time dependency of the z-component of position is zero; in other words, since the planet is in the plane z = 0 at some moment t₀, then it must always be in that plane.

Yet another way to prove the assertion of two-dimensional motion is to take a look at the time derivative of the cross product of position with momentum. The momentum, p, of the planet is defined as the planet’s mass times its velocity. The reason for taking the derivative will be clear in a moment. Using vector notation,

Notice that the first term on the right hand side is

since the cross product of a vector with itself is zero. Also, recall that the gravitational force vector and the position vector are anti-parallel. The cross product of any two co-linear vectors is zero; therefore

so,

Substituting these into (4), we have

Now, taking an indefinite time integral of both sides results in a constant of integration,

This means that if r and v are in a plane at some moment and then at some later time one of them (or the other or both) deviates from that plane – then their cross product would point in a different direction and thus not be a constant. Therefore r and v are restricted to one plane, and the two-dimensional nature of the orbit is proven.

Perhaps it seems as though equation four was written as if the end of the proof and the direction towards it were already known. They were, of course, but it does not make the proof any less rigorous. This seemingly post res nature of some mathematical proofs will appear again later.

The linear momentum of the planet is

Angular momentum, , is the cross product of position with linear momentum, so from equation (5) the angular momentum must be constant

Incidentally, since the angular momentum is a constant throughout the orbit, we could evaluate it at any point. At both perihelion, the point of closest approach to the sun, and aphelion, the furthest point in the orbit, the planet’s position vector is perpendicular to the velocity vector. Let r_p be the perihelion distance and v_p the velocity at perihelion, and let r_a be the aphelion distance and v_a the velocity at aphelion; the magnitude of angular momentum, then, is

Attempt at a perspective rendering of the x-y plane drawn co-planer with the r-v plane at time t₀. It is not necessary that either the x- or y-axis align with F, v, or indeed, r.

Polar coordinate diversion

Since the problem has been reduced to two-dimensions, and the force is centripetal, it would seem to be ideal to operate in polar coordinates (cylindrical coordinates, really, if ever we invoke the z-direction).

At this point it may be a useful diversion to step explicitly through the conversion between coordinate systems. Let θ be the angle between the position vector and the x-axis. We now have

The position vector is

so the unit vector in the direction of the position of the planet is

where the usual unit vector notations i and j in the x-y coordinate system are used. Therefore, the radial vector is written as

The dot notation for time derivative is employed here. Likewise, we have

and

Since the unit vector in the direction of increasing angle θ is in the x-y plane and is perpendicular to , it can be found from the negative (to give it the correct direction) of the cross product between and the unit vector in the z-direction, .

If required, these last two equations could be rearranged to give i and j in terms of

We may consider the two unit vectors in polar coordinates, , as moving along with the planet. This means that their direction is time dependent. Let’s take a look at the time derivative of the position unit vector.

Velocity

Now let’s use these polar coordinates to discuss the velocity of the planet before we return to the subject of angular momentum. The velocity of the planet is the time derivative of the position vector, namely,

[Or, equally, the time derivative in rectangular coordinates may be taken first, and then the complete conversion to polar coordinates is made:

Now substituting in the unit vectors we have

and with a little algebra the same result obtains.]

Angular Momentum

We know from equation (6) that the angular momentum is constant, so why not calculate the cross product? We have the position vector, equation (7), and we can use the velocity vector as expressed in equation (9).

Remove the dot notation and write the magnitude of the angular momentum as

By the way, dθ/dt = ω, is the planet’s angular velocity. Rearrange the expression for angular momentum to find

Now use equation (8) to write an expression for the magnitude of the velocity, again removing the dot notation:

Solve this for dr/dt, and substitute equation (10) in for dθ/dt:

Now, solve for dt

In order to show that the orbit is an ellipse, the time dependency of this equation needs to be removed and the equation should be written in terms of the angle. Use equation (10) once again; solve it for dt and substitute that into (11), thusly:

or,

Integration

Substitute that last expression into equation (12).

Now to find the position as a function of angle it is necessary to integrate both sides of equation (13). An indefinite integral will give the general result.

At this juncture it might be helpful to simplify the appearance of this equation by using the expression given above for angular momentum evaluated at either perihelion or aphelion. If we choose perihelion, say, then let r₀ = r_p and v₀ = v_p. We then have

This integral is of the form

and

where

and

Calculus diversion

A table of integrals would be handy in finding the solution to that integral. However, seeing this integration performed might prove instructive. Start with the trigonometric function cosine. Let x = cos(θ). A similar solution could be found if we were to begin with the sine function; at any rate, let's take the derivative of cosine with respect to θ.

or,

Now, notice that

and take the arcsine of both sides;

Here we take the derivative of both sides with respect to θ

Now, substitute x for cos(θ) and substitute equation (15) for dθ.

Using a trigonometric identity and again reverting to the variable x, we have

That’s the derivative of the arcsine. Let’s try a different argument in the arcsine — an inverse function — and include an arbitrary constant, - β; the chain rule for derivatives gives us

This is not quite the result we are after – but close. We can find what we need by using a slightly more complex argument in the arcsine function, and again using the chain rule.

Written in the form we recognize, this becomes

Integrate to obtain the final result

There we have it.

where

and

The Solution

The solution to equation (14) is now clear.

The condition which makes this solution possible is that the argument in the arcsine function must be between -1 and 1. This occurs if the numerator is less than the denominator:

Another solution to equation (14) is

This requirement is that

Now let us write the solution in terms of the orbit parameters. We will use the perihelion.

Now take the sine of both sides and solve for r.

Therefore,

Simplifying, this becomes

This is the equation which describes in polar coordinates the orbit as a conic section. The other solution to equation (14) produces a plus sign instead of a negative sign in the denominator (the one in front of the sine function). This puts the sun at the other focal point in the case of elliptical or hyperbolic orbits and opens the parabolic orbit in the opposite direction, but the overall shape of the curve would be unaltered. Both solutions, however, place the major axis of the curve along the y-axis. For the elliptical solution, this figure illustrates the general shape as written in equation (16):

The expression inside the square bracket is the eccentricity, ε, of the orbit. Rewriting, we have

where

Phase of the Ellipse

An inspection of the figure shows that the phase of the ellipse may at first seem a bit odd: perihelion occurs at θ = - π/2. Notice that if θ = 0 then the distance along the x-axis from the origin to the curve is given by

As mentioned above, there is a solution which might be found when starting with the cosine function. This can produce two expressions for r; one with (1 + ε cos(θ)), and one with (1 - ε cos(θ)) in the denominator. These solutions have the major axis aligned with the x-axis.

This cosine solution places either the aphelion or perihelion at θ = 0, which might be preferable, for exampe, when graphing the solution. However, a simple phase shift in the trig functions after finding the solution will accomplish the same rotation of the ellipse.

We might as well graph the expression for the distance to the sun versus angle just to see what it looks like. Here is the graph for Mercury, Mars, and Pluto over two cycles. Also, while we are at it, let’s look at θ versus r from perihelion to aphelion – over an angular difference of π. Note the asymmetry in both curves. This asymmetry indicates that the planet’s speed is not constant throughout its orbit. These graphs were made using the mass m₁ rather than M since the masses of these planets are so small compared with the sun. Also, the graphs of Mars and Pluto have been phase shifted to place the perihelion at θ = 0.

--------------------------------------------------------------------------------------------------------------------------------

Mercury; r versus θ, Equation (16)