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    m1 is the mass of the sun,

             m2 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 r1 and r2 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


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.3x1024 kg and the sun’s mass is 1.9855x1030 kg.  For this system, the reduced mass is 1893.2x1024 kg: a reduction of about 0.3%, and M = 1.9874x1030 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 t0, 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, az, 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 t0, which is zero at t = t0 because at that moment the velocity vector is in the x-y plane; we are therefore left with vz(t) = 0.  Now we integrate a second time,

Here, the constant of integration is z(t0) = 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 t0, 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


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 rp be the perihelion distance and vp the velocity at perihelion, and let ra be the aphelion distance and va 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 t0.  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


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.




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:





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 r0 = rp and v0 = vp.  We then have

This integral is of the form




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 θ.


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.



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.


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


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 m1 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)

Function written in the form for Graph program:




Mercury; θ versus r

Function written in the form for Graph program:




Mars; r versus θ, Equation (16)

Function written in the form for Graph program:




Mars; θ versus r

Function written in the form for Graph program:




Pluto; r versus θ, Equation (16)

Function written in the form for Graph program:




Pluto; θ versus r

Function written in the form for Graph program:



webpages and Eos image copyright 2015 M  Nealon

Orbital Types

Planetary Orbits


Let’s look up the perihelion distances and speeds of the planets and calculate the eccentricity.  The perihelion and aphelion values found on NASA webpages are listed in the table below.  NASA’s eccentricities are not measured values, but are calculated from orbital parameters.


Mass of the sun, m1= 1.988500 x 1030 kg



Universal Gravitational Constant, G = 6.67384 x 10-11 m3 / kg s2

(  [see also the note at the top of this NASA webpage for the sources of their information])

Graphs of Planetary Orbits, Equation (16)


On the scales of these graphs only Mercury, Pluto and perhaps Mars have readily obvious elliptical orbits.  The elliptical nature of the other orbits is a bit more subtle.

Top   Elliptical Orbits

Page 2  Time Dependency

Page 3  Kepler’s Equation

Using the Work-Energy Theorem


Equation (12) is the starting point in finding the shape of the orbit, and so we need an expression for the velocity, v.  The v2 term makes us think of the kinetic energy of the planet, so we might look at the work-energy theorem to aid us in finding the solution.  Let us say that the planet has some velocity v0 at some distance r0 and moves to some other (variable) position along its orbit and has a speed v and variable position r.  


The work done by the gravitational force in moving the planet from one position to another is equal to the change in the planet’s kinetic energy.  If the planet moves closer to the sun then the work done by the gravitational force is positive and the change in kinetic energy is also positive.  If the planet moves further away from the sun then the work done by the gravitational force is negative and the change in kinetic energy is also negative.  Either case gives us the same expression for v.

where M = m1 + m2; so the integral becomes

Using this and solving the work-energy theorem for v2 we have

Note:  The initial conditions r0 and v0 may occur at any place one may choose along the orbit, but typically these are evaluated at either the aphelion or perihelion.  At those places the position vector and velocity vectors are at right angles to each other.  Since aphelion, ra is the furthest distance from the sun, it is always greater than r, and the difference (1/r – 1/ra) is always positive.  However, since perihelion is the point of closest approach, the difference (1/r – 1/rp) is always negative.  This is of no consequence now, but later we must choose perihelion or aphelion carefully in order that the argument of an arcsine function is never greater than 1 or less than -1.

recalling that                and                   are parallel, so there is no need, really, for the dot product.  Tackle the integral first.  Write the magnitude of the force as

Notes on Eccentricity

This eccentricity is an interesting result:


1.    Notice that if we had chosen the aphelion as our initial condition, then the result would be, well, not exactly the same, but similar:

2.    Let us look at the results of the work energy theorem again.  Using the perihelion as initial condition we have

Now then, if the planet just happens to be at aphelion, we would let r = ra, and, correspondlingly set v = va, then

We also know, based on the idea of conservation of angular momentum, that (vp rp)2 = (va ra)2.  We could substitute at will either of these into the expression for eccentricity if we like.  Instead, however, let’s rewrite that expression and substitute for va 2.

Further manipulation gives us

and now,

Simplifying, we have

Subtract one from both sides and simplify to obtain the expression for ε.

Naturally, the aphelion expression for ε produces the same result.

The eccentricity is zero if ra = rp, or, if

This is the orbital speed of a planet in uniform circular motion.  Since

the orbital speed is constant, we can write the speed as vp = v.  In other

words, the orbit is a circle with radius

The orbital period, τ, is given by the circumference of the orbit

divided by the linear velocity.





If  ε = 0.





If  ε = 1.


The eccentricity is 1 if


and the orbit is parabolic.   The velocity is exactly that

of “escape”, so the body in orbit should never return;

however, a slight perterbation of the orbit may cause

the body to be captured.  [The escape velocity refered

to here is not the velocity of escape from the surface

of the sun, but the escape velocity at the point of

closest approach.]




If  ε > 1.


The eccentricity is greater than 1 if the velocity of

closest approach is greater than the corresponding

escape velocity:

In this case the orbit is hyperbolic.  The body will not return.

and the orbit is an ellipse with the sun at one

of the focal points.




If  0 < ε < 1.


If the eccentricity is between zero and 1, then


Perihelion distance

(x109 m)*

Perihelion velocity

(x103 m/s)*

Aphelion distance

(x109 m)*

Aphelion velocity

(x103 m/s)*


Calculated from (18)


(x1024 kg)*