Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 M.Z. Jorisch

Orbital

Perturbations

In this lesson, we will discuss the orbits of bodies in space, and how those bodies can be affected by others as they fly by. We will look at Encke's method, which was created by Johann Franz Encke in 1851, and uses a second order ODE to describe the true orbit of a body when affected by the pull of an additional body flying by.

In traditional orbital dynamics, the standard two-body problem is used to describe two bodies in motion with one orbiting the other. This fails to take into account the affect of outside bodies on the orbital trajectory and can produce an orbit very different than the "true" orbit.

These orbits play a large role in our daily lives. There are numerous satellites currently orbiting Earth, which are used for communications, GPS, as well as other data grabbers. These satellites can have slight changes to their orbits around Earth caused by other satellites, planets, moons, or comets that need to be taken into consideration when designing orbital parameters for them.

We will compare the traditional two-body motion and Encke's method to see how much the orbits vary over time.

For our example, we will use Mars and Jupiter orbiting the Sun, with Jupiter being the disturbing body. These two planets were chosen for their proximity to one another and for Jupiter's large mass (Over 300 times greater than Earth!)

Encke's Method

Figure 1. Visualization of Encke's Method ([Analytical Mechanics of Aerospace Systems Pg 342](http://www.control.aau.dk/~jan/undervisning/MechanicsI/mechsys/chapter10.pdf))

Encke's method uses the difference between the standard Keplerian orbit and the true orbit affected by a third body, and can be represented by the following equations.

The Keplerian orbit or the osculating orbit is represented by the equation:

$$\frac{d^2 \vec{r}_{osc}}{dt^2} = - \frac{\mu}{r^3 _{osc}} \vec{r}_{osc}$$

The perturbed orbit is represented by a similar equation:

$$\frac{d^2 \vec{r}}{dt^2} = - \frac{\mu}{r^3} \vec{r} + \vec{a}_d$$

When looking at the two equations, the only difference between the osculating orbit and the perturbed orbit is the term $\vec{a}_d$, which is the acceleration vector caused by the third body flyby.

The acceleration vector can be found using the following:

$$\vec{a}_d = \frac{1}{m_2} \vec{f}_{d_2} - \frac{1}{m_1} \vec{f}_{d_1}$$

The two accelerations cancel out many times. (Schaub pg 257)

This leads to:

$$\vec{a}_d = \frac{1}{m_2} \frac{G m_2 m_3}{|\vec{r}_23|^3} \vec{r}_23 - \frac{1}{m_1} \frac{G m_1 m_3}{|\vec{r}_13|^3} \vec{r}_13$$

Where $m_1$ is the mass of the central body, $m_2$ is the mass of the body orbiting around $m_1$, and $m_3$ is the mass of the disturbing body.

Initially, at time $t_0 = 0$, the osculating and perturbed orbits are equal. The change occurs at a time $t = t_0 + \Delta t$.

Let's define the difference between the radius of the osculating and perturbed orbits as $\delta$ and the difference between the velocities of the two orbits as $\gamma$

Therefore at time $t$, which we just defined, the radial and velocity components are:

$$\vec{\delta}(t) = \vec{r}(t) - \vec{r}_{osc} (t)$$$$\vec{\gamma}(t) = \vec{v}(t) - \vec{v}_{osc} (t)$$

We have some initial conditions as well. As mentioned before the obrbits are equal at $t_0$ which gives us $\vec{\delta} (t_0) = 0$. The velocity difference at $t_0$ is also zero. $\frac{d \vec{\delta} (t_0)}{dt} = \vec{\gamma} (t_0) = 0$

If we subtract our two initial equations we get:

$$\frac{d^2 \vec{\delta}}{dt^2} = \vec{a}_d + \mu \left( \frac{\vec{r}_{osc}}{r_{osc} ^3} + \frac{\vec{r}}{r^3} \right)$$

This can be simiplied to:

$$\frac{d^2 \vec{\delta}}{dt^2} + \frac{\mu}{r_{osc} ^3} \vec{\delta} = \frac{\mu}{r_{osc} ^3} \left( 1 - \frac{r_{osc} ^3}{r^3} \right) \vec{r} + \vec{a}_d$$

Our term $1 - \frac{r_{osc} ^3}{r^3}$ can be an issue because at the beginning of flight $r_{osc}$ and $r$ are approximately equal. That can't be too good, can it? We'll take a look at that a little later on.

We can find the radial and velocity components using the initial values of the radius and velocity along with the Legrangian coefficients in terms of the eccentric anomaly $E$, where the eccentric anomaly is the angle between the major axis and any point on the orbit.

$$\vec{r} = F \vec{r}_0 + G \vec{v}_0$$$$\vec{v} = \dot{F} \vec{r}_0 + \dot{G} \vec{v}_0$$

Where

$$F = 1 + \frac{a}{r_0} \left[ cos(E - E_0) - 1 \right]$$$$G = \frac{a \alpha _0}{\mu} \left[ 1 - cos(E - E_0) \right] + r_0 \sqrt{\frac{a}{\mu}} sin(E - E_0)$$$$\dot{F} = - \frac{\sqrt{\mu a}}{r r_0} sin(E - E_0)$$$$\dot{G} = 1 + \frac{a}{r} \left[ cos(E - E_0) - 1 \right]$$

(Tewari pg 104)

The eccentric anomaly, $E$, is equal to $M + e sin(E)$. $M$ is the mean anomaly and $e$ is the eccentricity.

Figure 2. Orbital anomalies for elliptic motion ([AIAA pg158])

As you can see from the equation, the eccentric anomaly is a function of itself. In order to calculate it we will have to start with a guess and then iterate until the difference between our new value of E and the guess is within a certain tolerance. This is based on Newton's approximation using a Taylor series expansion of $f(E) = E -e sin E - M$

The expansion is: $f(E + \Delta E) = \sum_{k = 0} ^{\infty} f^{(k)} (E) \frac{(\Delta E) ^k}{k!}$, in which $f(k) \approx \frac{d^k f(E)}{dE^k}$ The first two terms of the Taylor series can be used for Newton's approximation:

$$f(E + \Delta E) \approx f(E) + f^{(1)} (E) (\Delta E)$$

The following method and code for the eccentric anomaly approximation, are based upon those from Ashish Tewari's book, Atmospheric and Space Flight Dynamics. The initial guess we will start with uses the mean anomaly $M$ to give $E$:

$$E + M + e sin M$$

$\Delta E$ will be calculated using $-\frac{f(E)}{f^{(1)} (E)} = \frac{-E + e sin E + M}{1 - e cos E}$ so $f(E + \Delta E)$ is equal to $0$. After that, the value of E is updated, where $E = E + \Delta E$. This operation is repeated until a small enough difference is found, and our $E$ value is found.

Let's get to the coding!

The Kepler_eqn function uses the eccentricity, $e$, and mean anomaly, $M$, and then uses the iterative method to spit out our value of $E$.

The eccentric anomaly $E$ is the basis for the elliptical trajectory and is the value that changes at each time step, in turn changing the radial and velocity vectors.

Now let's define a function, that given our orbital parameters, will give us our trajectory at a time $t$.

The orbital elements that we will need to input into the function are:

• $a$, the semi-major axis, the distance from the center of the orbit to either of the foci
• $P$, the period, the time to complete one orbit
• $\mu$, the gravitational parameter, the gravitational constant $G$ times the mass of the body,
• $e$, the eccentricity
• $t_0$, an initial time
• $r_0$, an initial radial vector
• $v_0$, an initial velocity vector
• $t$, a time, where the new radial and velocity vectors will be found.

Using these initial values, we have equations that are set up within the function to give us our radial and velocity vectors. We also have to calculate our Legrangian coefficients, described earlier, $F$, $G$, $\dot{F}$, and $\dot{G}$

We will use the following equations:

• $\alpha _0 = r_0 \circ v_0$ : a constant that is used in the coefficient equations

• $\theta _0 = \pi$: our initial true anomaly (We are starting at the semi-major axis where it is equal to $\pi$)

• $n = \frac{2 \pi}{P}$: the mean motion, which is $2 \pi$ over the period

• $E_0 = 2 \tan ^{-1} (\sqrt{\frac{1 - e}{1 + e}} \tan (0.5 \theta _0))$: initial eccentric anomaly

• $\tau = t_0 + \frac{- E_0 + e \sin (E_0)}{n}$: where $t - \tau$ is the time since the closest point of orbit

• $M = n (t - \tau)$ : the mean anomaly which is used for the eccentric anomaly iteration in Kepler's equation

Now that we have a function to solve for the orbital trajectory of an elliptical orbit (which Mars and Jupiter both have) we have to create a function for the disturbing acceleration caused by the third body, Jupiter.

The equation for $a_d$ as mentioned earlier, is:

$$\vec{a}_d = \frac{1}{m_2} \frac{G m_2 m_3}{|\vec{r}_23|^3} \vec{r}_23 - \frac{1}{m_1} \frac{G m_1 m_3}{|\vec{r}_13|^3} \vec{r}_13$$

Our inputs for this function will be the mass of the Sun ($m_1$), the mass of Mars ($m_2$), the mass of Jupiter ($m_3$), and the radial vectors of Mars and Jupiter with respect to the Sun.

Now that we have our Kepler function, our elliptical orbit function, and our acceleration function, we can set up our initial conditions. After that we will be able to get our final radial and velocity vectors for the osculating and perturbed orbits in addition to Jupiter's orbit.

Initial Conditions

Next we will have a for-loop that will give us our solution at every time step for each of our orbits. In this case, we will use a final time of 4000 days. Each time step will be one day (in seconds), and we will be able to see the trajectories over that time period.

The loop uses a crude method of integration (Tewari pg 166) to solve for $\gamma$ and $\delta$, the difference between the osculating and perturbed radial and velocity vectors respectively.

The osculating orbit is used here as our base orbit. The orbit for Jupiter is calculated as well and used with the osculating orbit to get our disturbing acceleration. The acceleration is then used to find $\gamma$, which in turn is used to calculate $\delta$.

$\delta$ is added to the radial vector of the osculating orbit at every time step to give us our perturbed orbit. The same is done with $\gamma$ for the velocity vector for the perturbed orbit. The solutions are then entered into arrays so that we can plot them!

We mentioned earlier that sometimes the term $1 - \frac{r_osc ^3}{r^3}$ can cause issues because towards the beginning of orbit the two terms are approximately equal and can cause the solution to blow up, but this can be solved as follows:

$$1 - \frac{r_{osc} ^3}{r^3} = -B \frac{3 + 3B + B^2}{1 + (1 + B) ^{\frac{3}{2}}}$$

Where $$B = \frac{\vec{\delta} (\vec{\delta} - 2 \vec{r})}{r^2}$$

Now that we have our solutions arrays, we can plot our answers. Each array contains an x and y component for each time step. We will look at a plot of the osculating orbit and the perturbed orbit of Mars, and in a separate plot, we will look at those two orbits along with Jupiter's orbit.

Looking at our first plot we can see that there is a change in Mars' orbit due to the gravitational pull from Jupiter as it flies by!

Dig Deeper

See what happens when you change the orbital parameters! What happens with different planets or with a satellite orbiting earth? What about when both planets don't start at the zero point of the y axis?

References

• NASA Mars Fact Sheet

• NASA Jupiter Fact Sheet

• Standard Gravitational Parameter Wikipedia

• Perturbation (Astronomy)

• Battin, Richard H., AIAA Education Series, An Introduction to the Mathematics and Methods of Astrodynamics, AIAA, 1999

• Schaub, Hanspeter & John. L Junkins, AIAA Education Series, Analytical Mechanics of Aerospace Systems, AIAA, 2000

• Tewari, Ashish, Atmospheric and Space Flight Dynamics: Modeling and Simulation with MATLAB and Simulink, Birkhauser, 2007