View source for Astrodynamics/The Kepler Problem

{{astrodynamics}}

== The Kepler Problem ==
The Kepler Problem is a Two-Body Problem to find the position and velocity of an orbiting body over a period of time. Given the masses, position, and velocity, or alternatively the six Keplerian orbit elements. 

== Using Kepler's Orbital Elements ==
Given that the orbit is defined by the Semi-Major axis, <math>a</math>, the Eccentricity, <math>e</math>, the inclination, <math>i</math>, the Longitude of the Ascending Node, <math>\Omega</math>, and the Argument of Periapsis, <math>\omega</math>. It is possible to find the the position of an orbiting body at time, <math>t</math>. The procedure is as follows:

'''1. Compute the mean anomaly'''

<math>M = \sqrt{\frac{\mu}{a^3}} (t - \tau)</math>

Where, <math> \mu</math> is the gravitational parameter, <math> \tau</math> is the time at which the orbiting body passes the periapsis. For convenience the mean anomaly can be expressed in the interval <math>[0,2\pi)</math>.

'''2. Compute the eccentric anomaly using Kepler's equation'''

<math>M = E - e \sin{(E)}</math>

While Kepler's equation is easy to solve for time, there is no general solution for the reverse problem. To determine eccentric anomaly (and thus spacecraft position) at a given time, generally an iterative numerical method is used, such as Newton's method:
:<math> x_{k+1}=x_{k}-\frac{F(x_k)}{F'(x_k)}</math>
Where
:<math>F(E) = E - e \sin(E) - M(t)</math>
The iteration takes the form
:<math>E_{k+1}=E_{k} - \frac{ E_{k} - e \sin(E_{k}) - M(t) }{ 1 - e \cos(E_{k})}</math>
For most elliptical orbits an initial guess of ''E<sub>0</sub>'' = ''M'' is sufficient; for orbits with eccentricities greater than 0.8, ''E<sub>0</sub>'' = &pi; may be used. Better initial guesses are possible, but are not generally required. The iterative process is repeated until required accuracy conditions are achieved, for example:
:<math> |E_{k+1}-E_k|<\delta</math>
Where <math> \delta</math> is the desired accuracy.

Further, like with the mean anomaly the eccentric anomaly can be limited to the interval <math>[0,2\pi)</math>.

'''3. Compute the true anomaly'''

<math>\tan\frac{\nu}{2}=\sqrt{\frac{1+e}{1-e}}\tan\frac{E}{2}</math>

Note, that due to the range and domain of the tangent function and it's inverse a quadrant check is needed. If the eccentric anomaly has been limited to the interval <math>[0,2\pi)</math> then there is no need for a quadrant check when <math>E < \pi</math>.

== Using Cartesian Cooridnates ==
<math>\ddot{\boldsymbol{r}} = - \frac{\mu}{r^{2}} \hat{\boldsymbol{r}}</math>

From the special case of the Two-Body Problem the governing equation of the motion is shown above where <math>\mu</math> is the gravitational parameter. Given an initial position vector, <math>\boldsymbol{r} = (x, y, z)^\dotplus</math>, and corresponding velocity vector, <math>\dot{\boldsymbol{r}} = (u, v, w)^\dotplus</math> it is possible to solve for the future positions using numerical integration. To use ODE Runge-Kutta method the governing equation must be expressed as a system of first order ODEs as shown below. Note, that the system shown below is a scaler representation of the vectorised governing equation.  

<math>\begin{cases}
\frac{dx}{dt} = u \\ 
\frac{dy}{dt} = v \\ 
\frac{dz}{dt} = w \\
\frac{du}{dt} = - \frac{\mu x}{r^{3}} \\ 
\frac{dv}{dt} = - \frac{\mu y}{r^{3}} \\ 
\frac{dw}{dt} = - \frac{\mu z}{r^{3}} 
\end{cases}</math>

The Runge-Kutta method is an iterative semi-implicit numerical method, Many software's include this method such as MATLAB and in Python's SciPy integration method. These methods use an adaptive step where they estimate a local truncation error by using two orders of the method. For example RK45 would use the 4th and 5th orders, RK89 the 8th and 9th order. The error would be compared and checked against a user defined error and adjust the step size making it smaller or bigger depending if the error is too small or too large. A simple constant positive step size, <math>h</math>, the 4th order Runge-Kutta method is shown below.

Let a system of ODEs be represented as shown with an initial condition.

<math>\frac{d\boldsymbol{y}}{dt} = \boldsymbol{f}(t,\boldsymbol{y}),\qquad \boldsymbol{y}(t_{0}) = \boldsymbol{y}_{0}</math>

Then the 4th order Runge-Kutta method is formulate as

<math>\boldsymbol{y}_{n+1} = \boldsymbol{y}_{n} + \frac{h}{6} (\boldsymbol{k}_{1} + 2\boldsymbol{k}_{2} + 2\boldsymbol{k}_{3} + \boldsymbol{k}_{4})</math>

Where,

<math>\boldsymbol{k}_{1} = \boldsymbol{f}(t_{n}, \boldsymbol{y}_{n})</math>

<math>\boldsymbol{k}_{2} = \boldsymbol{f}(t_{n} + \frac{h}{2}, \boldsymbol{y}_{n} + h\frac{\boldsymbol{k}_{1}}{2})</math>

<math>\boldsymbol{k}_{3} = \boldsymbol{f}(t_{n} + \frac{h}{2}, \boldsymbol{y}_{n} + h\frac{\boldsymbol{k}_{2}}{2})</math>

<math>\boldsymbol{k}_{4} = \boldsymbol{f}(t_{n} + h, \boldsymbol{y}_{n} + h\boldsymbol{k}_{3})</math>

== Numerical Solution ==