Trajectory Simulation

This page covers the procedures used when simulating satellites' motion about primary objects and the formulae used to arrive at these individual solutions.

Below are the methods for determination of a satellite's motion based on the type of trajectory (circular, elliptical, parabolic, hyperbolic):

• ODE integration (any trajectory)
• Analytical (circular)
• Bessel summation (elliptical)
• Newton-Raphson (hyperbolic)

ODE Integration - Any Trajectory

Simulating a satellite's orbit using ordinary differential equation integration is perhaps the most physically intuitive method for estimating the motion of the satellite, as it relies on Newtonian gravitational physics. However, ODE integration is also subject to fluctuation as each subsequent data point is dependent on the solution to the previous points.

The only data necessary to compute the motion of a satellite is the system's standard gravitational parameter $\mu$, the Cartesian state vectors for position $\vec{r}$ and velocity $\vec{v}$ and a continuous array of times to perform integration for.

From the initial position $\vec{r_0}$ and velocity $\vec{v_0}$ of the satellite, the change in velocity in each direction is given by:

$$\begin{align} \vec{r_0} &= x_0 \hat{x} + y_0 \hat{y} + z_0 \hat{z} \tag{1.1} \\ \ddot{x} &= - \frac{\mu}{\lvert \vec{r}\rvert ^3} x_0 \hat{x} \tag{1.2} \\ \ddot{y} &= - \frac{\mu}{\lvert \vec{r}\rvert ^3} y_0 \hat{y} \tag{1.3} \\ \ddot{z} &= - \frac{\mu}{\lvert \vec{r}\rvert ^3} z_0 \hat{z} \tag{1.4} \\ \end{align}$$

Which then is fed back into the SciPy odeint function along with the initial velocity $\vec{v_0}$ for integration to determine the position and velocity vectors after time $\Delta t$. The "new" state vectors are then fed back into the equations for acceleration to continue the integration for every increment in the range of times being simulated.

Lastly, though not necessary for plotting, the Cartesian state vectors are used to determine the eccentricity vector $\vec{e}$:

$$\begin{align} \tag{2} \vec{e} = \frac{1}{\mu} \vec{v} \times (\vec{r} \times \vec{v}) - \frac{\vec{r}}{\lvert \vec{r}\rvert} \end{align}$$

Which can be used to determine the true anomaly $\nu$ of the satellite:

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

$$\begin{align} \tag{4} \nu &= \begin{cases} \theta & \text{if } \vec{e} \cdot \vec{r} < 0 \text{ or } \vec{e} \ge 1 \\ 360 - \theta & \text{otherwise} \end{cases} \end{align}$$

And the rest of the Keplerian elements can be computed using the cartesian_2_keplerian function.

Analytical - Circular Trajectories

Circular trajectories are the simplest orbits to predict as it is a special case of an elliptical orbit.

First, circular trajectories are unique in that the eccentricity $e=0$, and therefore the semi-major axis $a$ is identical to the satellite's radial distance $r$ from the system's center of mass. This necessarily implies that the magnitude of the orbital velocity $\lvert \vec{v} \rvert$ is constant.

From Kepler's Third Law, the orbital period of a satellite in a bound orbit ($e&lt;1$) is proportional to $a^3$:

$$\begin{equation} \tag{5} T^2 \propto a^3 \end{equation}$$

The mean angular motion of the satellite in its circular orbit is then closely related to the above, and is therefore useful in determining the true anomaly $\nu$ of the satellite along its orbit after $\Delta t$ span of time from $\nu_0$. Given that the mean anomaly $M$ is equivalent to the true anomaly, the mean angular motion $n$ of the satellite is proportional to the system's standard gravitational parameter $\mu$:

$$\begin{equation} \tag{6} n = \sqrt{\frac{\mu}{a^3}} \end{equation}$$

Then, Kepler's equation can be used in conjunction with the equation for mean angular motion to determine the true anomaly of the satellite from its initial mean anomaly $M_0$ (which in the case of circular orbits is identical to $\nu_0$) after length of time $\Delta t$:

$$\begin{align} M - M_0 &= n (t - t_0) \tag{7.1} \\ \nu - \nu_0 &= n \Delta t \tag{7.2} \\ \nu &= n \Delta t + \nu_0 \tag{7.3} \end{align}$$

Finally, the Cartesian state vectors representing the position and velocity vectors is calculated by supplying the keplerian_2_cartesian function with $\mu$ and the dictionary of Keplerian orbital elements for each increment of time.

Bessel Summation - Elliptical Trajectories

Elliptical trajectories are most accurately calculated by determining the initial mean eccentric anomaly $M_{E_0} = E_0$ and using the mean angular motion to estimate the eccentric anomaly. The mean anomaly $M$ in terms of the eccentric anomaly $E$ is given by:

$$\begin{align} M &= E - e \sin{E} \tag{8} \end{align}$$

And thus there is not an analytical solution to solve for the eccentric anomaly as a function of time and must be solved numerically. If it is assumed the epoch $t_0 = 0$ and the satellite is not at periapsis, the initial eccentric anomaly $E_0$ can be computed from the initial true anomaly $\nu_0$, which in turn can be expressed as the initial mean anomaly $M_0$:

$$\begin{align} M_0 &= \arctan2 \left( -\sqrt{1 - e^2} \sin{\nu_0}, -e-\cos{\nu_0}\right) + 180 - e \frac{\sqrt{1-e^2} \sin{\nu_0}}{1 + e \cos{\nu_0}} \tag{9} \end{align}$$

Which can be substituted into Kepler's equation:

$$\begin{align} M &= E - e \sin{E} + M_0 \tag{10.1} \\ M(t) &= n \Delta t + M_0 \tag{10.2} \end{align}$$

Due to the lack of a closed-form solution, the relations can be solved using the Newton-Raphson method and finding the root of the function:

$$\begin{align} f(E) &= 0 = E - e \sin{E} - M(t) \tag{11} \end{align}$$

Whose solution can be expressed as a Fourier series including a Bessel function of the first kind $J_n$ as:

$$\begin{align} E(t) &= M(t) + 2 \sum^{\infty}_{x=1}{\frac{J_n(xe)}{x} \sin(xM(t))} \tag{12} \end{align}$$

Once the eccentric anomaly at each time $t$ is computed, the range of values can then be converted back into the true anomaly with:

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

Finally, the Cartesian state vectors then can be determined with the keplerian_2_cartesian function.

Newton-Raphson Method - Hyperbolic Trajectories

Simulating hyperbolic trajectories by determination of the true anomaly as a function of time is more accurate than ODE integration albeit more computationally expensive. Hyperbolic trajectories also have certain unique caveats setting them apart from bound orbits: the semi-major axis is negative, the eccentricity is greater than 1, and the orbit is aperiodic, etc.

As usual, before simulation, it is necessary to determine the initial mean anomaly $M_0$ and hyperbolic anomaly $F_0$:

$$\begin{align} F_0 &= \text{sign} (\nu) \text{arccosh} \left( \frac{e + \cos \nu}{1 + e \cos \nu} \right) \tag{14} \end{align}$$

The hyperbolic anomaly also must have the same sign as the true anomaly, provided the true anomaly is negative and the satellite is approaching periapsis. With the initial hyperbolic anomaly, the corresponding mean anomaly is obtained from the Kepler equation for hyperbolic trajectories:

$$\begin{align} M_0 &= e \sinh F_0 - F_0 \tag{15} \end{align}$$

Next, the mean anomaly can be expressed as a function of time:

$$\begin{align} M(t) &= n \Delta t + M_0 \tag{16} \end{align}$$

The relation between the mean anomaly and the hyperbolic anomaly $F$, (sometimes $H$), must then be solved with the Newton-Raphson method and SciPy's newton algorithm. The function $f(F)$ and its derivatives then are passed to the algorithm:

$$\begin{align} f(F) &= e \sinh{F} - F - M(t) \tag{17.1} \\ f'(F) &= e \cosh{F} - 1 \tag{17.2} \\ f''(F) &= e \sinh{F} \tag{17.3} \\ \end{align}$$

After solving for the hyperbolic anomaly at each time being simulated, the next steps are slightly more complicated, as Astrodynamics handles the hyperbolic anomalies differently. When the results of the Newton-Raphson solver are returned, depending on where the satellite is along its trajectory, the returned anomaly may be outside the allowed domain:

$$\begin{align} |F| &\le \arccos\left(-\frac{1}{e}\right) \tag{18} \end{align}$$

For the indices in the solution array where the hyperbolic anomaly is not within the allowed domain, $360\degree$ are subtracted from the angle, making the hyperbolic anomaly at those indices negative, ensuring they are in the allowed domain. The solutions are then converted to their equivalent true anomaly angles via applying the hyperbolic_2_true_anom function to the absolute values of the estimated hyperbolic anomalies.

The array of true anomalies then must be made negative at the same indices as are negative in the array of hyperbolic anomalies, corresponding to a negative radial velocity. If the radial velocities are positive (and the satellite is moving away from the primary body), the hyperbolic anomalies can be converted to their equivalent true anomaly values without any special treatment. Why? Well:

1. It ensures that the angles are within the allowed domain
2. It is virtually identical to the plot of true anomalies derived from an equivalent simulation using ODE integration
3. I think the resulting true anomalies to make more sense in the context of gravitational encounters
4. It emphasizes the aperiodicity of hyperbolic trajectories
5. It works.

To elaborate, the true anomaly is defined as the angle between a satellite's argument of periapsis and the position of the satellite, meaning the satellite's true anomaly at periapsis is $0$. So, when the satellite enters a massive planet's gravitational sphere of influence, the true anomaly (when omitting the above outlined steps) will be positive, implying that the satellite has already surpassed periapsis. However, if the true anomaly when the satellite enters an object's sphere of influence is negative, the implication is that, as the mean anomaly monotonically increases as a function of time, the satellite is instead approaching periapsis. Additionally, with respect to the application of this convention to a hyperbolic trajectory, the sign of the true anomaly is also the sign of the radial velocity of the satellite.

Anyway, with the hyperbolic anomalies converted to the corresponding true anomalies, the state vectors can again be computed with the keplerian_2_cartesian function.

For more information on simulating hyperbolic trajectories, please see the source and inspiration the Astrodynamics algorithm was adapted from.