# How to find eccentric anomaly for orbit propagation

Given a set of orbit elements (assume inclination = 0, arg perigee = 0, RAAN = 0, e = known, and a = known), combined with a radius position and a velocity in orbit I should be able to forward propagate the position of the satellite by a given time t.

My current method:

$n = \sqrt{\frac{\mu}{a^3}}$

$E_0 = \mathrm{acos}((1-\frac{R}{a})\cdot(\frac{1}{e}))$

$M_0 = E_0 - e \cdot \mathrm{sin}(E_0)$

$M_1 = M_0 + n \cdot t$

$E_1 = M_1 + e \cdot \mathrm{sin}(E_1)$ <- for this I take an initial value of $E_1 = M_1$ and iterate

future positions:

$R = a \cdot (1-e \cdot \mathrm{cos}(E_1))$

$TrueA = 2 \cdot \mathrm{atan}(\sqrt{\frac{1+e}{1-e}} \cdot \mathrm{tan}(\frac{E_1}{2}))$

$V = \sqrt{\mu \cdot (2/R-1/a)}$

$fpa = \mathrm{atan}\frac{e \cdot \mathrm{sin}(TrueA)}{1+e \cdot \mathrm{cos}(TrueA)}$

I'm reasonably confident on every step except the finding the initial eccentric anomaly ($E_0$) and as it turns out this returns an error. The root of the error is trying to find the $\mathrm{acos}$ of a number less than -1 (first value I calculate is -1.0000000023241). I'm unsure of why this is occurring and any advice/help/alternative equations would be a great help.

• To any reasonable precision, -1 is the same as -1.0000000023241, so you could simply clamp the input to acos(). What are your starting R, a, e? – Russell Borogove Mar 4 '15 at 16:45

$E_0 = \mathrm{acos}((1-\frac{R}{a})\cdot(\frac{1}{e}))$
The other is the problem you discovered. This expression is subject to precision loss. Suppose you want to start at apofocus, so you set $R=a(1+e)$. The expression $(1-\frac{R}{a})\cdot(\frac{1}{e})$ should be -1, but using double precision numbers on my computer, this expression evaluates to 0 with $e=1\times10^{-16}$ and -1.97372982 with $e=1.125\times10^{-16}$.
This is a result of precision loss. (My cherry-picked values represent 100% precision loss.) You are always going to suffer some amount of precision loss when you subtract two numbers close to one another. For nearly circular orbits, $R/a$ will always be close to one.
Both problems can be solved by using a different expression. If you know the eccentricity $e$ and the initial true anomaly $f_0$, the initial eccentric anomaly $E_0$ is given by $E_0 = 2\arctan \left(\sqrt{\frac{1-e}{1+e}} \tan \frac f 2\right)$.
At the furthest from the primary, $R=a(1+e)$, so the argument of your $\mathrm{acos}$ should be exactly $-1$. It should never get larger. Your computation is $(1-(1+e))\cdot \frac 1e$, which has numeric cancellation if $e \ll 1$. Slight numeric inaccuracy between your values of $R,a,e$ can also be the problem. As suggested, you might just use if statements to bring the value inside $[-1,1]$