What is this algorithm to calculate orbital state vectors?

Question

I came across the following algorithm to calculate the cartezian vectors from the orbital elements:

First calculate some coefficients that will be used in determining the position: $$ r_x' = a (\cos E - \varepsilon)\\ r_y' = b \sin E\\ $$

Then calculate coefficients for the velocity: $$ \dot{\varepsilon} = \sqrt{\frac{\mu}{ar}}\\ v_x' = -a\dot{\varepsilon}\sin E\\ v_y' = b\dot{\varepsilon}\cos E $$

Finally, what appears to be some common coordinates for position and velocity: $$ p_x = \cos\omega \cos\Omega - \sin\Omega \sin\omega \cos i\\ p_y = \cos\omega \sin\Omega + \cos\Omega \sin\omega \cos i\\ p_z = \sin\omega \sin i\\ q_x = -\sin\omega \cos\Omega - \sin\Omega \cos\omega \cos i\\ q_y = -\sin\omega \sin\Omega + \cos\Omega \cos\omega \cos i\\ q_z = \cos\omega \sin i $$

Finally, position and velocity are calculated: $$ r_x = r_x'p_x + r_y'q_x\\ r_y = r_x'p_y + r_y'q_y\\ r_z = r_x'p_z + r_y'q_z\\ v_x = v_x'p_x + v_y'q_x\\ v_y = v_x'p_y + v_y'q_y\\ v_z = v_x'p_z + v_y'q_z $$

Where:

• $a$ and $b$ are the semi-major and semi-minor axes,
• $E$ is the eccentric anomaly,
• $\varepsilon$ is the eccentricity,
• $\mu$ is the gravitational parameter,
• $i$ is the inclination,
• $\omega$ is the argument of periapsis,
• $\Omega$ is the longitude of the ascending node, and
• $r$ is the distance from the barycenter to the orbiter.

Does anybody recognize this algorithm? I'd like to understand what $\dot{\varepsilon}$ represents, and how those coefficients are combined to produce position and velocity.

Answer 1

Then calculate coefficients for the velocity:$$\dot{\varepsilon} = \sqrt{\frac{\mu}{ar}}$$

The eccentricity ($\varepsilon$) is constant. That should be $\dot E$ rather than $\dot\varepsilon$, and that expression is incorrect. The correct expression is $$\dot E = \sqrt{\frac{\mu}{ar^2}}$$

---

I'd like to understand what $\dot{\varepsilon}$ represents, and how those coefficients are combined to produce position and velocity.

First, I'll look at what the eccentric anomaly represents. In terms of the eccentric anomaly, the orbital plane cartesian coordinates of an orbiting body are (as noted in the question) given by $$\begin{aligned}r_x &= a(\cos E - \varepsilon)\\r_y &= a\sqrt{1-e^2}\sin E\end{aligned}$$ The radial distance in terms of the eccentric anomaly is $$r = a(1- \varepsilon\cos E)$$

The time derivative of the eccentric anomaly can be calculated by differentiating Kepler's equation, $M = E - \varepsilon \sin E$, with respect to time: $\dot M = \dot E\,(1 - \varepsilon \cos E)$. Note the common term $1-\varepsilon \cos E$ in the expressions for the radial distance and the time derivative of Kepler's equation.

The time derivative of the mean anomaly is the mean motion: $\dot M = n = \sqrt{\frac \mu {a^3}}$. Thus the time derivative of the eccentric anomaly is $$\dot E = \frac {\dot M}{1-\varepsilon\cos E} = \sqrt{\frac \mu {a^3}} \frac 1 {1-\varepsilon\cos E} = \sqrt{\frac \mu {a^3}} \frac a r = \sqrt{\frac \mu {a r^2}}$$

Differentiating the coordinates $r_x$ and $r_y$ with respect to time yields $$\begin{aligned}\dot r_x = v_x &= \dot E a\sin E\\ \dot r_y = v_y &= \dot E a\sqrt{1-e^2}\cos E\end{aligned}$$