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 twobody 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 lefthandside 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 twodimensional space. This solution displays the elliptical nature of the orbit. Properly speaking, though, for a planet with a moon or moons, it is the centerofmass, 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_{1} is the mass of the sun, m_{2} 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_{1} and r_{2} 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 onebody 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^{24 }kg and the sun’s mass is 1.9855x10^{30 }kg. For this system, the reduced mass is 1893.2x10^{24 }kg: a reduction of about 0.3%, and M = 1.9874x10^{30} 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. 
Twodimensional nature of the orbit
That the orbital path of a planet about the sun is twodimensional 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 twodimensional nature of the orbits can be understood with a physical argument, as follows: the gravitational force vector is always antiparallel 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 freefall in which they are colinear). 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 (“centerdirected”) 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_{0}, and arrange the coordinate system such that the xy plane aligns with the rv plane. Since the component of force in the zdirection 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 zcomponent of speed at time t_{0}, which is zero at t = t_{0} because at that moment the velocity vector is in the xy 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}) = 0. This result means that the time dependency of the zcomponent of position is zero; in other words, since the planet is in the plane z = 0 at some moment t_{0}, then it must always be in that plane.
Yet another way to prove the assertion of twodimensional 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 antiparallel. The cross product of any two colinear 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 twodimensional 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 xy plane drawn coplaner with the rv plane at time t_{0}. It is not necessary that either the x or yaxis align with F, v, or indeed, r. 
Polar coordinate diversion
Since the problem has been reduced to twodimensions, and the force is centripetal, it would seem to be ideal to operate in polar coordinates (cylindrical coordinates, really, if ever we invoke the zdirection).
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 xaxis. 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 xy 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 xy 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 zdirection, . 
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, 
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_{0} = r_{p} and v_{0} = 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 yaxis. 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 xaxis 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 xaxis.
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_{1} 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: 5.5466E10/(1.2056sin(x))  
 Mercury; θ versus r Function written in the form for Graph program: asin((3.6058E11(2/x))/7.4203E12)  
 Mars; r versus θ, Equation (16) Function written in the form for Graph program: 2.2591E11/(1.09336sin(x+4.7135))  
 Mars; θ versus r Function written in the form for Graph program: asin((8.8531E12(2/x))/8.26516E13)+1.5708  
 Pluto; r versus θ, Equation (16) Function written in the form for Graph program: 5.519E12/(1.244sin(x+4.7135))  
 Pluto; θ versus r Function written in the form for Graph program: asin((3.6235E13(2/x))/8.8423E14)+1.5708  
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, m_{1}= 1.988500 x 10^{30} kg (http://nssdc.gsfc.nasa.gov/planetary/factsheet/sunfact.html)
Universal Gravitational Constant, G = 6.67384 x 10^{11} m^{3 }/ kg s^{2} (http://nssdc.gsfc.nasa.gov/planetary/factsheet/fact_notes.html [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 WorkEnergy 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 v^{2} term makes us think of the kinetic energy of the planet, so we might look at the workenergy theorem to aid us in finding the solution. Let us say that the planet has some velocity v_{0} at some distance r_{0} 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 = m_{1} + m_{2}; so the integral becomes 
Using this and solving the workenergy theorem for v^{2} we have 
Note: The initial conditions r_{0} and v_{0} 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, r_{a} is the furthest distance from the sun, it is always greater than r, and the difference (1/r – 1/r_{a}) is always positive. However, since perihelion is the point of closest approach, the difference (1/r – 1/r_{p}) 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 = r_{a}, and, correspondlingly set v = v_{a}, then 
We also know, based on the idea of conservation of angular momentum, that (v_{p }r_{p})^{2} = (v_{a }r_{a})^{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 v_{a }^{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 r_{a} = r_{p}, 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 v_{p} = 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. 
or, 
Circular
If ε = 0. 
Kepler. 
Parabolic
If ε = 1.
The eccentricity is 1 if 
or, 
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.] 
Hyperbolic
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. 
Elliptical
If 0 < ε < 1.
If the eccentricity is between zero and 1, then 

Perihelion distance (x10^{9} m)* 
Perihelion velocity (x10^{3} m/s)* 
Aphelion distance (x10^{9} m)* 
Aphelion velocity (x10^{3} m/s)* 
Eccentricity Calculated from (18) 
Mass (x10^{24} kg)* 
Mercury 
46.00 
58.98 
69.82 
38.86 
0.2058 
0.3301 
Venus 
107.48 
35.26 
108.94 
34.79 
0.006908 
4.8676 
Earth 
147.09 
30.29 
152.10 
29.29 
0.01690 
5.9726 
Mars 
206.62 
26.50 
249.23 
21.97 
0.09336 
0.64174 
Jupiter 
740.52 
13.72 
816.62 
12.44 
0.04937 
1898.3 
Saturn 
1352.55 
10.18 
1514.50 
9.09 
0.05590 
568.36 
Uranus 
2741.30 
7.11 
3003.62 
6.49 
0.0458 
86.816 
Neptune 
4444.45 
5.50 
4545.67 
5.37 
0.0130 
102.42 
Pluto** 
4436.82 
6.10 
7375.93 
3.71 
0.244 
0.0131 






