 Page 1  Elliptical Orbits
 Page 2  Time Dependency
 Page 3  Kepler’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 Q’ onto this new circle, as shown below, then the line drawn from Q’ 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),
 or,
 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):
 Therefore,
 then
 Now,        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 http://nssdc.gsfc.nasa.gov/planetary/text/eros.txt   ///////////////////////////////////////////////////////////////////////////////\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
 -----------------------------------------------------------------------------------------------------------------------------  (2*pi*x/tau)+2*epsilon*sin(2*pi*x/tau)+(5*(epsilon^2)/4)*sin(4*pi*x/tau) +((epsilon^3)/12)*(13*sin(6*pi*x/tau)-3*sin(2*pi*x/tau)) 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.   //////////////////////////////////////////////////////////////////////\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
 Mercury Mercury  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:  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:
 where 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).  Mercury This Mathematica code produces graphs for ψ(t), θ(t) and r(t):   "Mercury" ε = .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}] Plot[psi[t],{t,0,200}]   theta[t_] = 2*ArcTan[Sqrt[(1+ε)/(1-ε)]*Tan[psi[t]/2]] Plot[theta[t],{t,-87.969/2,87.969/2}]   r[t_] = a*(1-ε*ε)/(1+ε*Cos[theta[t]]) Plot[r[t],{t,0,200}]
 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:
 Mercury \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////////////////////////////////////
 ψ(t) versus t θ(t) versus t r(t) versus t
 Mars
 Mathematica code for Mars:   "Mars"  ε = .0934 m = 2 * 3.1415926 t /1.8809 "years" a = 1.5237*149500000000 "AU * meter/AU"   psi[ε_,m_] = m + Sum[(2/n)*BesselJ[n,nε]Sin[n*m],{n,1,12}] Plot[psi[ε,t],{t,0,2}]   theta[t_] = 2*ArcTan[Sqrt[(1+ε)/(1-ε)]*Tan[psi[ε,m]/2]] Plot[theta[t],{t,-1.88/2,1.88/2}]   r[t_] = a*(1-ε*Cos[psi[ε,m]]) Plot[r[t],{t,0,4}]
 ψ(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:   "Eros" ε = .2230 m = 2*3.1415926 t / 1.761 "years" a = 1.4583*149500000000 "AU * meter/AU"   psi[ε_,m_] = m + Sum[(2/n)*BesselJ[n,n*ε]Sin[n*m],{n,1,12}] Plot[psi[ε,t],{t,0,3.6}]   theta[t_] = 2*ArcTan[Sqrt[(1+ε)/(1-ε)]*Tan[psi[ε,m]/2]] Plot[theta[t],{t,-1.761/2,1.761/2}]   r[t_] = a*(1-ε*ε)/(1+ε*Cos[theta[t]]) Plot[r[t],{t,0,4}] ψ(t) versus t θ(t) versus t r(t) versus t
 Pluto   ψ(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}] Plot[psi[t],{t,0,3*tau},AxesLabel->{t,ψ}]   theta[t_] = 2*ArcTan[Sqrt[(1+ε)/(1-ε)]*Tan[psi[t]/2]] Plot[theta[t],{t,-tau/2,tau/2},AxesLabel->{t,θ}]   r[t_] = a*( 1-ε*ε)/(1+εCos[theta[t]]) Plot[r[t],{t,0,3*tau},AxesLabel->{t,r}]   ψ(t) versus t
 θ(t) versus t
 r(t) versus t Halley’s Comet, the nucleus Haley’s Comet, with Milky Way background  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   