This document presents my attempt to solve Kepler's Equation of Elliptical Motion due to Gravity.
Introduction
The geometric equation for an ellipse is quite simple; most high-school students are exposed to conic sections and their features.
However, when a body moves in an elliptical orbit, tracking its motion as a function of time is not as straightforward. The position of a body in elliptical motion is usually found by way of an intermediate property: the eccentric anomaly, E (refer to the diagram below).
Given an ellipse with semi-minor axis b and semi-major axis a, a circle of radius a is superposed upon the ellipse. Assume that at time, t, the body in elliptical motion is at point P. A vertical line is then drawn through P, up to where it intersects the superposed circle. Call this point Pc. The angle formed by drawing a line from Pc to the origin, to the body's position at periapse defines the eccentric anomaly, E. Once E is found, calculating the true position of the body (point P) is quite simple.
Determining the eccentric anomaly requires solving Kepler's Equation of Elliptical Motion:
E - e sin(E) = M
where E is the eccentric anomaly;
e is the eccentricity of the motion; and
M is the mean anomaly, 2πt/T (t is the time since periapse, and T is the period of the motion).
To date, an explicit solution for E has not been found. All solutions are found either numerically (by using computers to calculate a result to a desired accuracy) or by assuming that e is small and finding an approximation based on this assumption.
An Example
Planetary motion is elliptical. Plotting the path of a planet therefore requires solving Kepler's Equation of Elliptical Motion.
Consider Earth's motion. The Earth revolves in an elliptical orbit around the Sun, which is at one focus of the ellipse. The eccentricity of this ellipse is about 0.0167. In addition, the time it takes for the Earth to make one complete orbit of the Sun is approximately 365 days. We also know that the Earth approaches closest to the Sun approximately January 4 of each year; in other words, perihelion occurs approximately January 4. Now say we want to know the Earth's position February 1, or 27 days since perihelion. Finding the Earth's position at this time requires solving Kepler's Equation of Elliptical Motion. In this case,
e = 0.0167, T = 365 days, and t = 27 days.
Therefore, M = 2*3.14159*27/365 = 0.46478 radians.
Finding the Earth's position at February 1 requires solving Kepler's Equation of Elliptical Motion for E:
E - 0.0167 sin(E) = 0.46478
References
Kepler's Equation is presented in many books about Celestial Physics, Celestial Mechanics, and Spacecraft Dynamics. These books give a short introduction and then present the equation in the form presented above; this form is the most concise form into which Kepler's Equation has so far been resolved.
For more detailed coverage of attempts that have been made to solve Kepler's Equation, the following is an excellent book:
"Solving Kepler's Equation over Three Centuries"
by Peter Colwell
Published by Willmann-Bell, Inc.
Richmond, Virginia 1993
This book presents attempts that have been made to solve the equation by eminent scientists and mathematicians since Kepler's Equation was first established. In addition, it
presents approximations that have been proposed and some numerical methods used today that, with the aid of a computer, solve Kepler's Equation to any desired degree of accuracy. (If you should need one, the page at the following link is a numerical utility that computes solutions of Kepler's Equation:
Calculator for Kepler's Equation of Elliptical Motion.)
Yet, to date, an exact, theoretical, solution has not been found. My attempt follows.
My Approach
Using data generated by a numerical method, my approach is to analyze the form of the curve and reduce Kepler's Equation to a composition of simple factors.
A plot of sin(E) goes to zero each time t takes a value of nT/2, where n is an integer. Because I do not want to introduce any "magic factors" into my solution to take this feature into account, I assume that the solution to sin(E) includes a factor of sin(M):
sin(E) = Y(e, t) sin(M)
Plotting Y(e, t) = sin(E)/sin(M) indicates certain features.
Specifically, the plot implies that Y(e, t) = n(e, t)/[1 - e cos(M)]. This proposal satisfies conditions for t = nT/2.
In turn, plotting n(e, t) indicates that it takes the form n(e, t) = √[1 - r(e, t)].
At this point, Kepler's Equation has been reduced to the form
sin(E) = {√[1 - r(e, t)]/[1 - e cos(M)]} sin(M)
The task now is to find r(e, t).
Continuing in this manner, I have worked toward a solution for sin(E). I have not yet reached a final solution, but I feel I have made progress.
In what follows, I will use P(e, t) as my proposed solution.
This work is ongoing. If you know of any published work in which an attempt has been made to solve Kepler's Equation for an exact solution, I would like to hear
from you. Or if you have any suggestions, please send me an e-mail.
Does anybody know of an exact, explict expression for sin(E T/4)? In other words, sin(E), or E, when the time is one-quarter of the period?
I believe that knowing an exact expression for sin(E T/4) would provide a skeleton on which to build the rest of the solution.
In fact, sin(E T/4) is given by the following equation, which must be solved numerically to determine sin(E T/4) for a particular value of e:
sin(E T/4) = √(1 - sin ²(e sin(E T/4)))
Rearranging, this equation can be restated as follows:
Define θ = sin(E T/4).
Then sin(E T/4) for a particular e is found as the solution to
cos(e θ) = θ
Either of these expressions is verified by numerical data.
Here are a few data points that can be used to check proposed expressions for sin(E T/4):
At e = 0, sin(E T/4) = 1
For e = 0.001, sin(E T/4) = 0.999999500001
For e = 0.1, sin(E T/4) = 0.995053426874
For e = 0.3, sin(E T/4) = 0.958906950862
For e = 0.5, sin(E T/4) = 0.90036722259
For e = 0.542199961348, sin(E T/4) = 0.886651779307 = √ 0.786151377748
For e = 0.7, sin(E T/4) = 0.834269951868
For e = 0.773634829186, sin(E T/4) = 0.81
For e = 0.786151377748, sin(E T/4) = 0.805917329563, cos(E T/4) = -0.592028088785
For e = 0.9, sin(E T/4) = 0.769576421734
For e = 0.963645663751, sin(E T/4) = 0.75
For e = 0.999, sin(E T/4) = 0.739382684206
As e approaches 1, sin(E T/4) approaches 0.739085133215
I have been posting information about much of my current work on my blog, for example, Keplers Equation: Eccentric Anomaly Values at the Quarter-period