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.

  • 1
    $\begingroup$ 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? $\endgroup$ – Russell Borogove Mar 4 '15 at 16:45

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


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.