# Kepler

In orbital mechanics, Kepler's equation relates various geometric properties of the orbit of a body subject to a central force. This can be used to calculate the position of a body (relative to the central body) at any time, without having to iteratively integrate acceleration and velocity over time.

In Ephemeris all gravity sources (planets, moon etc) are kept to Kepler orbits, only free floating items (spaceships, satellites, missiles etc) are plotted step-by-step.

## Orbital Elements

A keplerian orbit is defined by six **Orbital Elements**:

- $e$ = Eccentricity
- $a$ = Semi Major Axis
- $i$ = Inclination
- $Ω$ = Longitude of the Ascending Node
- $ω$ = Argument of periapsis
- $M_0$ = Mean Anomaly At Epoch

## Calculating Position (Overview)

To calculate position the process is:

- Choose a
**time**since epoch. - Calculate
**Mean Anomaly**at this time. - Calculate
**Eccentric anomaly**from**Mean anomaly**. - Calculate Position (X/Y) from
**Eccentric Anomaly**. - Tilt X/Y Position into orbital plane.

## Choose A Time

The **True Anomaly** element defines the position of the body along the orbit at an arbitrary time, called the "epoch". To calculate a position a time, relative to the epoch, must be chosen.

## Calculate Mean Anomaly

The mean anomaly is a measure of how far along the orbit the body is at a particular time. Since the body completes a full orbit every $orbital\_period$ (by definition) then the **Average Angular Velocity** ($n$) is:

$n = \frac{360\degree}{orbital\_period}$

The **Mean Anomaly** is simply:

$M = M_0 + n * time\_since\_epoch$

## Calculate Eccentric Anomaly

**Mean Anomaly** is related to **Eccentric Anomaly**:

$M = E - e \sin E$

However, there is no closed form solution to this equation. Instead an iterative solution such as Newton's Method to refine an initial guess.

$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$

In our case:

$f(E) = E - e * sin(E) - M$

$f'(E) = 1 - e * cos(E)$

Implementing this in code:

`## Calculate Eccentric anomaly from M (Mean Anomaly) and (e) eccentricity`

def solve_kepler(M, e):

f = lambda E: E - e * sin(E) - M

E = initial_guess(M, e) # Make an initial guess

delta = 10000

while (delta > 0.001): # exit once refinement process isn't doing much

e1 = E - f(E) / (1 - e * cos(E)) # refine the guess

delta = abs(e1 - E) # calculate how much the guess has improved

E = e1

return E

Now we just need an initial guess. The simplest option is simply to use the mean anomaly ($M$) as the initial guess.

For the original source of this code and also a better initial approximation, see this blog post.

## Calculate Position

Finally **Eccentric Anomaly** can be used to calculate position:

$x = a * (cos(E) - e)$

$y=b*sin(E)$

## Orbital Plane

We now have the position in 2D, on the orbital plane. This simply needs to be tilted into the orbital plane to achieve a 3D position:

`_quaternion = Quaternion.AngleAxis(LongitudeOfAscendingNode, float3(0, 0, 1))`

* Quaternion.AngleAxis(Inclination, float3(0, 1, 0))

* Quaternion.AngleAxis(ArgumentOfPeriapsis, float3(0, 0, 1));

var pos3 = _quaternion * new Vector3(x, y, 0);

This is where the final orbital parameters are used, they define a plane which tilts the $x$,$y$ position into 3d.