How to find eccentric anomaly for orbit propagation

Question

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.

Answer 1

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

This has two problems. One problem is that using the inverse cosine only gives you half of the orbit.

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)$.

Answer 2

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]$