Given orbital parameters and position vector along that orbit, how can I determine the eccentric anomaly at that point?

Question

For example:

Consider a prograde Earth orbit with $\mu = 3.98588738 \times 10^{14} m^3s^{-2}$ semi-major axis $a= 4 \times 10^6 m$ and an eccentricity $\epsilon = 0.3$. Let's ignore the argument of periapsis $\omega$ since the we can just rotate the input coordinates by it.

For a given point along that orbit, say $(2.418542, 1.626223) \times 10^6 m$ (about 126 seconds after periapsis), the eccentric anomaly is about 0.440273 radians.

How can I compute the eccentric anomaly?

To get position from eccentric anomaly, I'm using the first two equations that I found in this answer:

$$x=a\left(\cos\tau-e\right)$$

$$y=a\sqrt{1-e^2}\sin\tau$$

But rearranging either of those equations to isolate $\tau$ gives me a value totally off from where it should be.

I'm working in 2D if that matters. Greatly appreciate any help.

Answer 1

From the first equation, we get $\cos\tau = \frac x a + e$. This on its own is not quite enough to find $\tau$, since a given value of $\cos\tau$ usually corresponds to two values of $\tau$ in $[0,2\pi)$. To find the correct value, we can use the fact that for a prograde orbit with the argument of periapsis $0$, $0\leq \tau\leq \pi$ if $y\geq 0$, and $\pi<\tau<2\pi$ if $y<0$. So we have $$ \tau = \begin{cases} \arccos\left(\frac x a + e\right)&\text{if $y\geq 0$};\\ 2\pi - \arccos\left(\frac x a + e\right)&\text{if $y < 0$}. \end{cases} $$ Substituting the values from your example, we have $\tau = \arccos(2.418542/4 +0.3) = 0.440273$, just as it should be. Alternatively, you can express $\tau$ in terms of $y$ from the second equation. It's better to use this approach for positions very close to the periapsis or the apoapsis as it give better precision in these cases: $$ \tau = \begin{cases} \arcsin\left(\frac y {a\sqrt{1-e^2}}\right)& \text{if $y\geq 0$ and $x\geq -ae$};\\ \pi - \arcsin\left(\frac y {a\sqrt{1-e^2}}\right)& \text{if $x < -ae$};\\ 2\pi + \arcsin\left(\frac y {a\sqrt{1-e^2}}\right)& \text{if $y < 0$ and $x\geq -ae$}. \end{cases} $$ Again, substituting your values, we get $\tau = \arcsin(1.626223/(4\sqrt{0.91})) = 0.440273$.

Answer 2

This may be a little thin for an answer but maybe I or someone else can expand it from a start.

Take note of the first figure in the Wikipedia article on Eccentric Anomaly. You have r (position to point P in the diagram) as a vector. You have a and e as orbital parameters, from which you can locate C and F in the figure (the center of the orbital ellipse and the focus for the central body).

Once you know all of those you can compute the eccentric anomaly with standard trigonometric methods.

Answer 3

Use the eccentricity vector, $$ \vec e = \frac{\vec v \times (\vec r \times \vec v)}{\mu} - \frac{\vec r}{||\vec r||} = \left( \frac1r - \frac1a \right) \vec r - \frac{\vec r \cdot \vec v}{\mu}\,\vec v$$

This vector is a constant of motion. Its magnitude is the eccentricity of the orbit and its direction is along the line from the center of the central body to the periapsis point. The eccentricity vector is thus useful for determining both the eccentricity and the argument of periapsis.

You need the argument of periapsis to compute the true anomaly given the position of some orbiting body, and you need the eccentricity to compute the eccentric anomaly given the true anomaly.

To compute the true anomaly, first compute the argument of latitude, which is the sum of the true anomaly and the argument of periapsis. In 2D, the relationship between argument of latitude and the position of an orbiting body in a coordinate system that has the central body at the origin is given by $\tan u = r_y/r_x$, where $r_x$ and $r_y$ are the $x$ and $y$ coordinates of the orbiting body. Use a two argument inverse tangent (e.g., $u = \mathrm{atan2}(y, x)$) to compute $u$. Then subtract the argument of periapsis to yield the true anomaly: $\nu = u - \omega$. Finally, to compute the eccentric anomaly, use $$\tan \frac\nu2 = \sqrt\frac{1+e}{1-e} \tan\frac E2$$