Skip to content

Conversion Functions

Jump to bottom
Evan Bauer edited this page Apr 23, 2024 · 7 revisions

Below are the various conversions functions currently defined in Astrodynamics, including conversion between Keplerian orbital elements and Cartesian state vectors (and vice versa) as well as more straightforward conversions:

Trajectory determination

Conversions between Cartesian the state vectors (position $\vec{r}$ and velocity $\vec{v}$) and the six Keplerian orbital elements (semi-major axis $a$, eccentricity $e$, inclination $i$, longitude of the ascending node $\Omega$, argument of periapsis $\omega$, and true anomaly $\nu$) are a fair bit more complicated conversions than those discussed later on this page, as they contain many more degrees of freedom.

The relevant sources for these functions will, of course, be provided, as neither of these functions are totally written from scratch by myself, though I have made fairly significant alterations to the original code to best suit this project.

Keplerian to Cartesian

The keplerian_2_cartesian function takes a two-body system's standard gravitational parameter and a dictionary of Keplerian orbital elements and computes the equivalent Cartesian state vectors.

The first thing keplerian_2_cartesian does is verify that the signs of the semi-major axis and eccentricity are consistent with their implications:

$$\begin{align} e > 1 &\iff a < 0 \tag{1.1} \\ e < 1 &\iff a > 0 \tag{1.2} \end{align}$$

If the sign of the semi-major axis is inconsistent with that implied by the value of the eccentricity, the semi-major axis will be changed accordingly. Then, keplerian_2_cartesian computes the magnitude of the specific angular momentum $h$ depending on the value of $e$:

$$\begin{align} \tag{2} h^2 &= \begin{cases} \mu a (1 - e^2) & \text{if } e < 1 \\ -\mu a (e^2 - 1) & \text{if } e > 1 \end{cases} \end{align}$$

Next, we need to evaluate the length of the semi-latus rectum $\ell$:

$$\begin{align} \tag{3} \ell &= \begin{cases} a & \text{if } e = 1 \\ a (1 - e^2) & \text{if } e < 1 \\ -a (e^2 - 1) & \text{if } e > 1 \end{cases} \end{align}$$

Which is then used to determine the magnitude of the position vector $r$ at the corresponding true anomaly:

$$\begin{align} r &= \frac{\ell}{1 + e \cos{\nu}} \tag{4} \end{align}$$

Now, both the position and velocity vectors must be defined with respect to the perifocal frame, denoted by the subscript $p$::

$$\begin{align} \vec{r_{p}} &= r \langle \cos \nu, \sin \nu, 0 \rangle \tag{5.1} \\ \vec{v_{p}} &= \frac{\mu}{h} \langle - \sin \nu, e + \cos \nu, 0 \rangle \tag{5.2} \end{align}$$

With the position and velocity vectors in the perifocal frame defined, they can be applied to a rotation matrix given by the three Euler angles ($\omega$, $i$, and $\Omega$). The rotation matrix is formed by the proper Euler angle $Z\text{-}X\text{-}Z$ sequence of rotations. To accomplish the rotation from the perifocal frame to the equatorial frame, Astrodynamics utilizes SciPy's Rotation algorithm. The rotation is initialized by indicating the rotation sequence and supplying the opposite of each angle. Finally, SciPy provides a method for applying the given rotation to vectors, so we can apply the inverse rotation to the vectors in the perifocal frame and the algorithm returns the relevant position and velocity vectors in the equatorial frame, thus completing the conversion from Keplerian elements to Cartesian state vectors.

Code originally adapted from: https://orbital-mechanics.space/classical-orbital-elements/orbital-elements-and-the-state-vector.html#orbital-elements-state-vector

Cartesian to Keplerian

The cartesian_2_keplerian function takes a two-body system's standard gravitational parameter and a dictionary of Cartesian state vectors and computes the corresponding Keplerian orbital elements.

The conversion from state vectors to Keplerian elements is fairly straightforward and doesn't require any fancy-pants rotations. To start, we need to evaluate the specific angular momentum vector $\vec{h}$ and specific orbital energy $\epsilon$:

$$\begin{align} \vec{h} &= \vec{r} \times \vec{v} \tag{6.1} \\ \epsilon &= \frac{\vec{v}^2}{2} - \frac{\mu}{\vec{r}} \tag{6.2} \end{align}$$

With the specific angular momentum vector, the eccentricity vector $\vec{e}$ can be calculated and the value of $e$ determined:

$$\begin{align} \vec{e} &= \frac{\vec{v} \times \vec{h}}{\mu} - \frac{\vec{r}}{\lvert \vec{r} \rvert} \tag{7.1} \\ e &= \lvert \vec{e} \rvert \tag{7.2} \end{align}$$

The length of the semi-major axis is related to the specific orbital energy:

$$\begin{align} a &= - \frac{\mu}{\epsilon} \tag{8} \\ \end{align}$$

The inclination can be computed by comparing the magnitudes of the z-component of $\vec{h}$ to the magnitude of $\vec{h}$:

$$\begin{align} i &= \arccos{\frac{h_z}{\lvert \vec{h} \rvert}} \tag{9} \\ \end{align}$$

Next, the longitude of the ascending node can be calculated by defining the convenience vector $\vec{N}$ by computing the cross-product of a reference unit vector (chosen to be $+\hat{z}$) and the specific angular momentum:

$$\begin{align} \vec{N} &= \hat{z} \times \vec{h} \tag{10} \\ \end{align}$$

The longitude of the ascending node $\Omega$ is dependent on the x- and y-components of $\vec{N}$:

$$\begin{align} \tag{11} \Omega &= \begin{cases} \arccos \frac{N_x}{\lvert \vec{N} \rvert} & \text{if } N_y \ge 0 \\ 2\pi - \arccos \frac{N_x}{\lvert \vec{N} \rvert} & \text{if } N_y < 0 \end{cases} \end{align}$$

Next, the argument of periapsis is dependent on the sign of the eccentricity vector's z-component as follows:

$$\begin{align} \tag{12} \omega &= \begin{cases} \arccos \left(\frac{\vec{e} \cdot \vec{N}}{\lvert \vec{e} \rvert \lvert \vec{N}\rvert}\right) & \text{if } e_z \ge 0 \\ 2\pi - \arccos \left(\frac{\vec{e} \cdot \vec{N}}{\lvert \vec{e} \rvert \lvert \vec{N}\rvert}\right) & \text{if } e_z < 0 \end{cases} \end{align}$$

Lastly, the true anomaly, which will first be denoted by $\theta$, can be determined directly from the eccentricity and state vectors:

$$\begin{align} \tag{13} \theta &= \arccos \left( \frac{\vec{e} \cdot \vec{r}}{\lvert \vec{e} \rvert \lvert \vec{r} \rvert} \right) \end{align}$$

For the "true" true anomaly $\nu$, we must interpret $\theta$ based on the context of the orbit. In the case of hyperbolic orbits, $e &gt; 1$, the sign of $\nu$ is also the sign of the radial velocity of the satellite:

$$\begin{align} \tag{14} \nu &= \begin{cases} \theta & \text{if } \vec{r} \cdot \vec{v} > 0 \\ -\theta & \text{if } \vec{r} \cdot \vec{v} < 0 \end{cases} \end{align}$$

And in the case of bound orbits:

$$\begin{align} \tag{15} \nu &= \begin{cases} \theta & \text{if } \vec{r} \cdot \vec{v} > 0 \\ 2 \pi - \theta & \text{if } \vec{r} \cdot \vec{v} < 0 \end{cases} \end{align}$$

Code adapted from: https://orbital-mechanics.space/classical-orbital-elements/orbital-elements-and-the-state-vector.html#state-vector-orbital-elements

Anomaly conversion

The following functions are each defined in Astrodynamics and are specifically for converting between various anomaly types.

Mean to true

Caution

This section is pending review.

The mean_2_true_anom function converts a mean anomaly to the equivalent true anomaly via:

$$\begin{align} \text{Let } \beta &= \frac{1 - \sqrt{1 - e^2}}{e} \tag{16.1} \\ \text{and } f(j) &= \sum^{\infty}_{k=1}{\beta^j (J_v(j-k, je) + J_v(j+k, je))} \tag{16.2} \\ \end{align}$$

$$\begin{align} \nu = M + 2 \sum^{\infty}_{j=1} \frac{\sin (Mj)}{j} \left(J_v(j, je) + f(j) \right) \tag{16.3} \\ \end{align}$$

Where $J_v$ is the Bessel function.

Mean to eccentric

Converting the mean anomaly to the eccentric anomaly $E$ we have the mean_2_eccentric_anom function which computes the eccentric anomaly using the relation:

$$\begin{align} E = M + 2 \sum^{\infty}_{n=1} \frac{J_n(n, ne)}{n} \sin(Mn) \tag{17} \\ \end{align}$$

Where $J_n$ is the Bessel function of the first kind.

Mean to hyperbolic

The mean_2_hyperbolic_anom function converts the mean anomaly to the hyperbolic anomaly $F$ using the Newton-Raphson method of finding the root of the function $f(F)$:

$$\begin{align} f(F) &= e \sinh{F} - F - M \tag{18.1} \\ f'(F) &= e \cosh{F} - 1 \tag{18.2} \\ f''(F) &= e \sinh{F} \tag{18.3} \\ \end{align}$$

True to eccentric

Converting the true anomaly to the eccentric anomaly is handled by the true_2_eccentric_anom function following the relation:

$$\begin{align} E &= \arctan \left( \sqrt{1 - e^2} \frac{\sin \nu}{e + \cos \nu} \right) \tag{19} \end{align}$$

True to hyperbolic

Converting the true anomaly to the hyperbolic anomaly is handled by the true_2_hyperbolic_anom function following the relation:

$$\begin{align} F &= 2 \text{arctanh} \left( \sqrt{\frac{e-1}{e+1}} \tan \frac{\nu}{2} \right) \tag{20} \end{align}$$

Eccentric to true

Converting the eccentric anomaly to the true anomaly is handled by the eccentric_2_true_anom function following the same relation as in Eq. 19:

$$\begin{align} \text{Let } \beta &= \frac{e}{1 + \sqrt{1 - e^2}} \tag{21.1} \\ \nu &= E + 2 \arctan \left( \frac{\beta \sin E}{1 - \beta \cos E} \right) \tag{21.2} \end{align}$$

Hyperbolic to true

Converting the hyperbolic anomaly to the true anomaly is handled by the hyperbolic_2_true_anom function following the same relation as in Eq. 20:

$$\begin{align} \nu &= 2 \arctan \left( \sqrt{\frac{e + 1}{e - 1}} \text{tanh} \frac{F}{2} \right) \tag{22} \end{align}$$

Hyperbolic to true

Caution

This function is not defined in the current version of Astrodynamics, but is anticipated to be included in a future release.

Converting the hyperbolic anomaly to the mean anomaly is handled by the hyperbolic_2_mean_anom function and is equivalent to evaluating the hyperbolic Kepler equation at the given hyperbolic anomaly:

$$\begin{align} M &= e \text{sinh} F - F \tag{23} \end{align}$$

Convenience conversion functions

This last section covers the few conversion functions included for convenience.

Degrees to radians

The radians function takes an angle measured in degrees and returns the equivalent angle in radians within the domain $0 \le \theta &lt; 2 \pi \text{ rad}$:

$$\begin{align} \theta_{rad} &= \mathrm{mod} \left( 2 \pi \frac{\theta_{deg}}{360}, 2 \pi \right) \tag{24} \end{align}$$

Radians to degrees

The degrees function takes an angle measured in radians and returns the equivalent angle in degrees within the domain $0 \le \theta &lt; 360\degree$:

$$\begin{align} \theta_{deg} &= \mathrm{mod} \left( 360\frac{\theta_{rad}}{2\pi}, 360 \right) \tag{25} \end{align}$$

Vector to unit vector

The unit function takes a vector $\vec{v}$ and returns a parallel vector with a magnitude of $1$ (the unit vector $\hat{v}$):

$$\begin{align} \hat{v} &= \frac{\vec{v}}{\lvert \vec{v} \rvert} \tag{26} \end{align}$$

Clone this wiki locally