<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step1:  Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 M.Z. Jorisch
<h1 align="center">Orbital</h1>
<h1 align="center">Perturbations</h1>
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!)
<h2 align="center">Encke's Method</h2>
<h4 align="center">Figure 1. Visualization of Encke's Method ([Analytical Mechanics of Aerospace Systems Pg 342](http
Step3:  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
Step5:  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
Step6:  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.
<h3 align="center">Initial Conditions</h3>
Step7:  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!
Step8:  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
Step9:  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!
<h3 align="center">Dig Deeper</h3>
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?
<h2 align="center">References</h2>
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
<ASSISTANT_TASK:>
Python Code:
from matplotlib import pyplot
import numpy
from numpy import linalg
%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
def Kepler_eqn(e, M):
    Takes the eccentricity and mean anomaly of an orbit to solve Kepler's equation
    
    Parameters:
    ----------
        e : float
            eccentricity of orbit
        M : float
            Mean anomaly of orbit
        
    Returns:
    -------
        E : float
            Eccentric anomaly
    
    
    E = M + e * numpy.sin(M) # eccentric anomoaly
    fofE = E - e * numpy.sin(E) - M #eccentric anomaly as a function of E
    fdotE = 1 - e * numpy.cos(E) #derivative with respect to E of fofE
    dE = - fofE / fdotE # change in E
    Enew = E + dE
    tolerance = 1e-2
    
    while abs(fofE) > tolerance:
        E = M + e * numpy.sin(Enew)
        fofE = E - numpy.sin(E) - M
        fdotE = 1 - e * numpy.cos(E)
        dE = - fofE / fdotE
        Enew = E + dE
    
    return E
    
    #Based on code from Ashish Tewari
Explanation: Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 M.Z. Jorisch
<h1 align="center">Orbital</h1>
<h1 align="center">Perturbations</h1>
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!)
<h2 align="center">Encke's Method</h2>
<h4 align="center">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))</h4>
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.
<h4 align="center">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!
End of explanation
def ellip_orb(a, Period, mu, e, t0, r0, v0, t):
    
    Calculates the orbital position for an elliptical orbit
    
    Parameters:
    ----------
    
    a  : float
        Semi-major axis
    Period : float
        Period of planetary orbit
    mu : float
        Gravitational parameter
    t0 : float
        Initial time t = 0
    r0 : array of float
        Initial positional array
    v0 : array of float
        Initial velocity array
    t : float
        time
    
    Returns:
    -------
    
    r : array of float
        Array of radius at each time t
    v : array of float
        Array of velocity at each time t
    
    
    r0_norm = numpy.linalg.norm(r0) # Normalized initial radius
    
    v0_norm = numpy.linalg.norm(v0) # Normalized initial velocity
    
    alpha = r0 * v0 # Constant used for Legrangian coefficients
    
    theta0 = numpy.pi # Initial true anomaly
            
    n = 2 * numpy.pi / (Period) # Mean motion, given the period
    E0 = 2 * numpy.arctan(numpy.sqrt((1 - e) / (1 + e)) * numpy.tan(0.5 * theta0)) # Initial eccentric anomaly
    
    tau = t0 + (- E0 + e * numpy.sin(E0)) / n # t - tau is time since Perigee
    M = n * (t - tau) # Mean anomaly
    E = Kepler_eqn(e, M) # Eccentric anomaly found through Kepler's equation
    
    r_leg = a * (1 - e * numpy.cos(E)) # Radius used for legrangian coefficients
    
    #Legrangian Coefficients
    
    F = 1 + a * (numpy.cos(E - E0) - 1) * r0_norm
    
    G = a * alpha * (1 - numpy.cos(E - E0)) / mu + r0_norm * numpy.sqrt(a / mu) * numpy.sin(E - E0)
    
    F_dot = - numpy.sqrt(mu * a) * (numpy.sin(E - E0)) / (r_leg * r0_norm)
    
    G_dot = 1 + a * (numpy.cos(E - E0) - 1) / r_leg
    
    r = numpy.zeros_like(r0)
    
    v = numpy.zeros_like(v0)
    
    r = F * r0 + G * v0 # Radial value of orbit for specified time
    
    v = F_dot * r0 + G_dot * v0 # Velocity value of orbit for specified time
    
    return r, v
Explanation: 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
End of explanation
def acceleration_d(m1, m2, m3, r, r3):
    
    Calculates the acceleration due to the disturbing orbit
    
    Parameters:
    ----------
    m1 : float
        Mass of central body
    m2 : float
        Mass of second body
    m3 : float
        Mass of third (disturbing) body
    r : array of float
        Radial distance between body two and one
    r3: array of float
        Radial distance between body three and one
        
    Returns:
    -------
    a_d : array of float
        Acceleration due to the disturbing orbit
    
    a_d = numpy.zeros((2, 1))
    
    G = 6.674e-11 # Gravitational constant
    
    r13 = r3 # Radial distance between Jupiter and the Sun
    
    r23 = r - r3 # Radial distance between Jupiter and Mars
    
    r23_norm = numpy.linalg.norm(r23) # Normalized radius between Jupiter and Mars
    
    r13_norm = numpy.linalg.norm(r13) # Normalized radius between Jupiter and the Sun
    
    a_d = (((1 / m2) * ((G* m2 * m3)/ (r23_norm ** 3))) * r23) - (((1 / m1) * ((G * m1 * m3) / (r13_norm ** 3))) * r13)
    
    return a_d
Explanation: 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.
End of explanation
mu3 = 1.2669e17                  # Standard gravitational parameter of Jupiter in m^3 / s^2
m3 = 1.8983e27                   # Mass of Jupiter in kg
e3 = .0489                       # Eccentricity of Jupiter
a3 = 778000000.                  # Semi-major Axis of Jupiter in km
Period3 = 4332.589 * 3600 * 24   # Period of Jupiter Orbit in seconds
mu = 4.2828e13                   # Standard gravitational parameter of Mars in m^3 / s^2
m2 = 6.4174e23                   # Mass of Mars in kg
e = .0934                        # Eccentricity of Mars
a = 228000000.                   # Semi-major Axis of Mars in km
Period = 686.980 * 3600 * 24     # Period of Mars Orbit in seconds
mu1 = 1.3271e20                  # Standard gravitational parameters of the Sun in m^3 / s^2
m1 = 1.989e30                    # Mass of the Sun in kg
dt = 24 * 3600                   # Time step
tfinal = 4000 * dt               # Final time
N = int(tfinal / dt) + 1         # Number of time steps
t = numpy.linspace(0,tfinal,N)       # Time array
r0 = numpy.array([228000000., 0.])   # Initial radius of Mars
v0 = numpy.array([-21.84, -10.27])   # Initial velocity
r3_0 = numpy.array([778000000., 0.]) # Initial radius of Jupiter
v3_0 = numpy.array([-13.04, -.713])  # Initial velocity of Jupiter
# Set arrays for radial and velocity components
r = numpy.empty((N, 2))
v = numpy.empty((N, 2))
gamma =  numpy.empty((N, 2))
delta = numpy.empty((N, 2))
r_n = numpy.empty((N, 2))
v_n = numpy.empty((N, 2))
a_d = numpy.empty((N, 2))
r_osc = numpy.empty((N, 2))
r_osc_n = numpy.empty((N, 2))
v_osc = numpy.empty((N, 2))
v_osc_n = numpy.empty((N, 2))
r3_n = numpy.empty((N, 2))
Explanation: 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.
<h3 align="center">Initial Conditions</h3>
End of explanation
for i,ts in enumerate(t):
    delta = numpy.zeros_like(r0)
    gamma = numpy.zeros_like(r0)
    r_osc, v_osc = ellip_orb(a, Period, mu1, e, t[0], r0, v0, ts) # Trajectory of the osculating orbit of Mars
    r_osc_norm = numpy.linalg.norm(r_osc) # Normalized osculating orbit of Mars
    r0_norm = numpy.linalg.norm(r0) # Normalized initial orbit of Mars
    r3, v3 = ellip_orb(a3, Period3, mu3, e3, t[0], r3_0, v3_0, ts) # Trajectory of Jupiter
    a_d = acceleration_d(m1, m2, m3, r_osc, r3) # Acceleration due to Jupiter
    gamma = mu3 * (dt) * ((1 - (r_osc_norm / r0_norm) ** 3) / r_osc_norm ** 3) + a_d * (dt) # Difference in velocity between osculating orbit and perturbed
    delta = gamma * (dt) # Difference between osculating orbit and perturbed orbit radius
    r = r_osc + delta # Perturbed orbital radius
    v = v_osc + gamma # Perturbed orbital velocity
    r_osc_n[i,:] = r_osc # Value of osculating orbital radius for every time step
    v_osc_n[i,:] = v_osc # Value of osculating orbital velocity for every time step
    r3_n[i,:] = r3 # Value of Jupiter's radius for every time step
    r_n[i,:] = r # Value of the perturbed orbital radius for every time step
    v_n[i,:] = v # Value of the perturbed orbital velocity for every time step
    
Explanation: 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!
End of explanation
x = numpy.linspace(t[0], t[-1], N)
pyplot.figure(figsize = (10,10))
pyplot.grid(True)
pyplot.xlabel(r'X Distance (km)', fontsize = 18)
pyplot.ylabel(r'Y Distance (km)', fontsize = 18)
pyplot.title('Trajectory of Osc vs Perturbed Orbit, Flight Time = %.2f days' %(tfinal / dt), fontsize=14)
pyplot.plot(r_n[:,0], r_n[:,1])
pyplot.plot(r_osc_n[:,0], r_osc_n[:,1])
pyplot.legend(['Perturbed Orbit', 'Osculating Orbit']);
pyplot.figure(figsize = (10,10))
pyplot.grid(True)
pyplot.xlabel(r'X Distance (km)', fontsize = 18)
pyplot.ylabel(r'Y Distance (km)', fontsize = 18)
pyplot.title('Trajectory of Osc, Perturbed and Jupiter Orbit, Flight Time = %.2f days' %(tfinal / dt), fontsize=14)
pyplot.plot(r_n[:,0], r_n[:,1])
pyplot.plot(r_osc_n[:,0], r_osc_n[:,1])
pyplot.plot(r3_n[:,0], r3_n[:,1])
pyplot.legend(['Perturbed Orbit', 'Osculating Orbit', 'Jupiter\'s Orbit']);
Explanation: 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.
End of explanation
from IPython.core.display import HTML
css_file = '../numericalmoocstyle.css'
HTML(open(css_file, "r").read())
Explanation: 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!
<h3 align="center">Dig Deeper</h3>
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?
<h2 align="center">References</h2>
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
End of explanation
<END_TASK>