Page 1Elliptical Orbits

Page 2Time Dependency

Page 3Keplerís Equation

Keplerís Equation

This formula, derived at the bottom of the previous page, is called Keplerís equation.† The planetís time-in-orbit, t, is given as a function of its angular position, θ.† Keplerís idea to solve the equation for the angular position as a function of time (despite the transcendental nature of the problem) was to use the area swept out by the planet as a proxy for a clock.† This is possible since, as indeed Kepler was first to observe, the areal velocity is a constant throughout the orbit.


In order to accomplish this same task let us begin with the x- and y-coordinates of the planetís position, P, along the ellipse.† The origin of our x-y coordinate system is at the sun, and the distance from the sun to center of the ellipse is aε.† In Cartesian coordinates an ellipse with this off-set and with semi-major axis a and semi-minor axis b is given by

Next, if we escribe a circle of radius b inside the ellipse and project the planetís position parallel to the x-axis to a point onto this new circle, as shown below, then the line drawn from to the center of the ellipse makes the same angle ψ with respect to the x-axis.† (We leave this assertion as an exercise for the reader to prove.)

We can see these relationships geometrically, as follows: Kepler inscribed a circle around the ellipse with the radius a, and projected the position of the planet perpendicular to the x-axis onto the circle, point Q in this figure.† The parameter ψ is a measure of the angle between the x-axis and the line from the center of the ellipse to Q.† From the figure we see that x is related to the cosine of the angle ψ by

so we know right away that we could characterize the ellipse with a parameter ψ, where

From this figure, we have the y-coordinate,

What we want is the relationship between ψ and θ.† Recall that (somewhere or other) we had

These two equations give us r as a function of ψ,

so we can write the magnitude of the position vector as

or, as we prefer for now

Add εr to both sides of equation (26), and recalling that (1-ε2) = (1-ε)(1+ε), we have

Now we can substitute in the ψ dependence of r from equation (25) into the right-hand-side only and rearrange:

Subtract εr from both sides of equation (26), and simplify in a similar manner to obtain

Divide equation (27) by (26); this leads to

At this point we can use the half-angle formula for the tangent function, which in general is

With this, then, we now have

Solve this last for ψ:

This equation gives us the angular position of the planet as measured from the center of the ellipse, given any angular position θ as measured from the sun.† This happens to be the first term on the right-hand-side of Keplerís equation (24).


At this juncture we wish to put Keplerís equation into the standard form found in all the texts.† In order to do so, we could either re-derive Keplerís equation with this variable change from θ to ψ using integral calculus, or we could simply transform the rest of Keplerís equation (24) as already written using algebra only.† We shall do both, just for the exercise.

The fraction in the second term of Keplerís equation is

A)†††††††† Transform Keplerís Equation in Variable θ into Variable ψ


In order to transform the second term on the right-hand-side of Keplerís equation (24) [we have transformed the first term already], we begin by solving equation (29) for θ and taking the sine and cosine of θ .† Let us introduce the momentarily useful symbol γ.

We now have

and consequently,

where we have used the double-angle formulś for the sine and cosine functions, which are, in general,

Now substitute in the expression for γ

Therefore, Keplerís equation becomes

B)†††††††† Derive Keplerís Equation by Integration Using Variable ψ


Since we are going to use this equation, (23) from the previous page (substituting in for b),


we will need the differential dθ.† Begin by solving equation (29) for θ and using, temporarily, the symbol δ.† If

and we have this one, equation (27), rearranged and having had the half-angle formula applied to the cosines:

since the secant is the inverse of the cosine.† Substituting the original variable, θ, back into the solution may seem odd, but it removes the δ in the denominator and in just a moment we will cancel it anyway:

The integral weíre after solving includes an r2 term, so we will find two separate expressions for r and multiply them together. We have this one, equation (25):




A)† Series Expansion


For Keplerís equation in this form

How to Use Keplerís Equation


Solutions to Keplerís equation give the planetís angular position, either θ or ψ, as a function of time.† From that, then, the radial distance, r, as a function of time may be calculated with equation (25) or (25í).† For a parametric representation in (r, θ) of the orbit in polar coordinates, time may be used as the parameter.† The degree to which any of these approximate solutions are either appropriate or sufficiently accurate for any particular application is left as a decision for the reader.† As always, an appeal to nature is warranted.††


This is a look at several techniques for finding a solution to the equation:

and the integral is

which is, after simplifying,

Thus, upon integrating and rearranging, we have again Keplerís equation in the standard form,

This function is plotted below in red for Mercury, Pluto, the asteroid Eros, and the dwarf planet Ceres.† The blue curves are plots of θ(t) with all terms except the linear first term (2πt/τ); these blue curves are included in the graphs simply to illustrate the asymmetric oscillatory nature of the function a bit more clearly.

this series expansion in θ works for planetary orbits if terms in ε3 are included, but still-and-all it works best for small ε.† Without derivation, we have

ďReconstructedĒ image from NASA


Eros Fact Sheet
Size: 33 km x 13 km x 13 km
Approximate mass: 7.2 x 10^15 kg
Rotation Period: 5.270 hrs
Orbital Period: 1.76 yrs
Spectral Class: S
Semi-major Axis: 1.458 AU
Perihelion Distance: 1.13 AU
Aphelion Distance: 1.78 AU
Orbital Eccentricity: 0.223
Orbital Inclination: 10.8 deg
Geometric Albedo: 0.16






tau = 87.969 earth day†††††††† †††††††††††††††† epsilon = .2056


ďGraphĒ program.


Graph of θ versus t over a few orbits for Mercury (in green, left axis) and a graph of r versus t (violet, right axis) using

where† a = 5.787x1010 m.


ďLogger-ProĒ program.

Parametric representation of Mercuryís orbit in r and θ

as calculated above using Series Expansion.





Streaked Conical Mountain on Ceres Captured by Dawn Spacecraft, NASA

Ceres, Keck Observatory

Ceres, NASA, Hubble Telescope


The red curve is the series expansion in θ for Pluto.† The blue curve is a plot of θ(t) with all terms except the linear first term, 2πt/τ.



Series expansion in θ for Eros the asteroid.†

B)† Besselís function



According to Encyclopśda Britannica, Besselís functions were ď. . . derived around 1817 by the German astronomer Friedrich Wilhelm Bessel during an investigation of solutions of one of Keplerís equations of planetary motion.Ē

Letís begin with Keplerís equation in this form:

to answer the following questions:

The solution can be found using the sum of Besselís functions of the first kind, Jk, multiplied by a sine function.†† Note that although the index k begins at 1 and goes to , the series may be sufficiently accurate if truncated after only a few terms.

With x = kε, we shall write the first few terms explicitly.


For k = 1,

A close encounter with Besselís functions of the first kind reveals that they are each themselves a series.† When k is a positive integer, and for some argument x, they are given, in general, by

For k = 4,

where the factorial is

For k = 3,

For k = 2,

etc., etc.† So the solution becomes a tedious accounting of the sum of all the sums.† To make the calculation easier we should find a mathematics program, Mathematica for instance, which will call up Besselís functions as subroutines.† For that particular program, here is the code:


Psi[t_] is the function to be calculated, namely, ψ(t);

BesselJ[k,kε]is the Bessel function subroutine for Jk(kε);

this is multiplied by the sine function, Sin[k2πt/τ];

and the bracket {k,1,12} sets the limits of index for the sum from k = 1 to k = 12.


In Mathematica, if we input this code:

this output will be returned:

Specifying the value for ε will cause the program to calculate and return a result which, if then added to 2πt/τ, will give the function ψ(t).


This Mathematica code produces graphs for ψ(t), θ(t) and r(t):



ε = .2056

m = (2 * 3.1415926 t)/87.969

a = 57871000000


psi[t_] = m + Sum[(2/k)*BesselJ[k,k*ε ]Sin[k*m],{k,1,12}]



theta[t_] = 2*ArcTan[Sqrt[(1+ε)/(1-ε)]*Tan[psi[t]/2]]



r[t_] = a*(1-ε*ε)/(1+ε*Cos[theta[t]])


Notes on Mathematica code:


††† a)††††† In using the notation psi[t_], the blank space underline is required when defining the function on the left-hand-side of the equal sign.† Drop the blank space underline when using the function in calculations.

††† b)††††† All subroutines and embedded functions in Mathematica must have upper case first letter, for example, Sin, Tan, ArcTan, Sqrt, etc., and their arguments must be placed in square brackets.

††† c)†††† The bracket {t,0,200} in the† Plot command places the variable, t, on the x-axis, and sets the minimum (0) and maximum (200) limits of the x-axis.† The y-axis is auto-scaled with this command.

††† d)†††† A right-arrow is generated by a minus sign followed by a greater-than sign:† ->

††† e)†††† There are many ways to go wrong in Mathematica.† If for some (unknown) reason these short programs do not work as exactly written here, youíll just have to figure out Mathematica for yourself ó as you should anyway, of course.


††† Those commands produce these results for:



ψ(t) versus t

θ(t) versus t

r(t) versus t


Mathematica code for Mars:



ε = .0934

m = 2 * 3.1415926 t /1.8809


a = 1.5237*149500000000

"AU * meter/AU"


psi[ε_,m_] = m + Sum[(2/n)*BesselJ[n,nε]Sin[n*m],{n,1,12}]



theta[t_] = 2*ArcTan[Sqrt[(1+ε)/(1-ε)]*Tan[psi[ε,m]/2]]



r[t_] = a*(1-ε*Cos[psi[ε,m]])


ψ(t) versus t

θ(t) versus t

Halleyís Comet Fact Sheet


Aphelion 35.1 AU * 1.495x1011m/AU = 5.247x1012 m

Perihelion .586 AU* 1.495x1011m/AU = 8.761x1010 m

eccentricity† ε = .967

Period† τ = 75.3 year

Mass† =† 2.2x1014 kg

For a smooth curve this orbit requires many terms in the sum of Besselís functions

n = number of terms


Note that the eccentricity is quite high.

r(t) versus t

Eros the asteroid

Mathematica code for Eros:



ε = .2230

m = 2*3.1415926 t / 1.761


a = 1.4583*149500000000

"AU * meter/AU"


psi[ε_,m_] = m + Sum[(2/n)*BesselJ[n,n*ε]Sin[n*m],{n,1,12}]



theta[t_] = 2*ArcTan[Sqrt[(1+ε)/(1-ε)]*Tan[psi[ε,m]/2]]



r[t_] = a*(1-ε*ε)/(1+ε*Cos[theta[t]])


ψ(t) versus t

θ(t) versus t

r(t) versus t


ψ(t) versus t

θ(t) versus t

r(t) versus t

Halleyís Comet

Mathematica code for Halleyís Comet:


"Halley's Comet"

"ε eccentricity"

ε = .967

"tau† Period"

tau = 75.3

m = 2*π*t /tau

"a†† semi-major axis"

a = 17.8*149500000000

"requires many terms for smooth curve"

"n = number of terms in Bessel's Function"

n = 48



"Sum[2/kBesselJ[k,k ε]Sin[k m],{k,1,n}]"

psi[t_] = m+Sum[2/k BesselJ[k,k ε]Sin[k m],{k,1,n}]



theta[t_] = 2*ArcTan[Sqrt[(1+ε)/(1-ε)]*Tan[psi[t]/2]]



r[t_] = a*( 1-ε*ε)/(1+εCos[theta[t]])


ψ(t) versus t

θ(t) versus t

r(t) versus t

Halleyís Comet, the nucleus

Haleyís Comet, with Milky Way background

Text Box: The Comet, in 1986

C)† Trial and error


Consider a hypothetical planet with period τ = 100 days and eccentricity ε = .25, with our clock set such that t = 0 occurs when the planet is at perihelion. Let the semi-major axis be a = 5x109 m.† Use Keplerís equation in this form:

1)†††† How far along in its orbit is the planet after a time of t = 50 days?

2)†††† How far along is this same planet after a time of t = 57.335 days?

3)†††† How far along is the planet at i)† t = 17.126 days, and at ii)† t = 25 days?† We leave this as an exercise, for the reader will be greatly edified by calculating ψ, θ, and r for each of the specified times.

Oops!, ran out of room.† Click here to continue to Part 2:


Keplerís Equation, Part 2

Keplerís Equation, Part 2

webpages and Eos image copyright 2015 M  Nealon

Top†† Keplerís Equation

Page 1†† Elliptical Orbits

Page 2†† Time in Orbit

Keplerís Equation, Part 2