Going to Jupiter with Python using Jupyter and poliastro

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

import astropy.units as u
from astropy.time import Time
from astropy.coordinates import solar_system_ephemeris

from poliastro.bodies import Sun
from poliastro.twobody import Orbit
from poliastro.maneuver import Maneuver

from poliastro.iod import izzo
from poliastro.ephem import get_body_ephem

from poliastro.plotting import plot, OrbitPlotter

from poliastro.util import norm
In [2]:
solar_system_ephemeris.set("jpl")
Out[2]:
<ScienceState solar_system_ephemeris: 'jpl'>
In [3]:
## Initial data
# Links and sources: https://github.com/poliastro/poliastro/wiki/EuroPython:-Per-Python-ad-Astra
date_launch = Time("2011-08-05 16:25", scale='utc')
C_3 = 31.1 * u.km**2 / u.s**2
date_flyby = Time("2013-10-09 19:21", scale='utc')
date_arrival = Time("2016-07-05 03:18", scale='utc')
In [4]:
r_e0, v_e0 = get_body_ephem("earth", date_launch)
In [5]:
r_e0
Out[5]:
$[1.0246553 \times 10^{8},~-1.023135 \times 10^{8},~-44353346] \; \mathrm{km}$
In [6]:
v_e0
Out[6]:
$[1847708.5,~1594323.4,~691089.12] \; \mathrm{\frac{km}{d}}$
In [7]:
# Initial state of the Earth
ss_e0 = Orbit.from_vectors(Sun, r_e0, v_e0, epoch=date_launch)
In [8]:
# State of the Earth the day of the flyby
r_efly, v_efly = get_body_ephem("earth", date_flyby)
ss_efly = Orbit.from_vectors(Sun, r_efly, v_efly, epoch=date_flyby)
In [9]:
# Assume that the insertion velocity is tangential to that of the Earth
dv = C_3**.5 * v_e0 / norm(v_e0)
man = Maneuver.impulse(dv)
In [10]:
# Inner Cruise 1
ic1 = ss_e0.apply_maneuver(man)
ic1.rv()
Out[10]:
(<Quantity [  1.02465527e+08, -1.02313505e+08, -4.43533465e+07] km>,
 <Quantity [ 2198705.82621214, 1897186.74383867,  822370.88977492] km / d>)
In [11]:
ic1.period.to(u.year)
Out[11]:
$2.1515475 \; \mathrm{yr}$
In [12]:
op = OrbitPlotter()
op.plot(ss_e0)
op.plot(ic1)

plt.grid(linestyle='--')
plt.gca().set_axisbelow(True)  # Please tone down that grid
../_images/examples_Going_to_Jupiter_with_Python_using_Jupyter_and_poliastro_12_0.png
In [13]:
# We propagate until the aphelion
ss_aph = ic1.propagate(ic1.period / 2)
ss_aph.epoch
Out[13]:
<Time object: scale='utc' format='iso' value=2012-09-01 14:38:56.604>
In [14]:
# Let's compute the Lambert solution to do the flyby of the Earth
time_of_flight = date_flyby - ss_aph.epoch
time_of_flight
Out[14]:
<TimeDelta object: scale='tai' format='jd' value=403.1958726432301>
In [15]:
(v_aph, v_fly), = izzo.lambert(Sun.k, ss_aph.r, ss_efly.r, time_of_flight)
In [16]:
# Check the delta-V
norm(v_aph - ss_aph.v)  # Too high!
Out[16]:
$1.0798665 \; \mathrm{\frac{km}{s}}$
In [17]:
ss_aph_post = Orbit.from_vectors(Sun, ss_aph.r, v_aph, epoch=ss_aph.epoch)
ss_junofly = Orbit.from_vectors(Sun, r_efly, v_fly, epoch=date_flyby)
In [18]:
op = OrbitPlotter()
op.plot(ss_e0, label="Earth")
op.plot(ic1, label="Inner Cruise 1")
op.plot(ss_efly, osculating=False)
op.plot(ss_aph_post, label="Back to Earth")

plt.grid(linestyle='--')
plt.gca().set_axisbelow(True)
../_images/examples_Going_to_Jupiter_with_Python_using_Jupyter_and_poliastro_18_0.png
In [19]:
# And now, go to Jupiter!
r_j, v_j = get_body_ephem("jupiter", date_arrival)
ss_j = Orbit.from_vectors(Sun, r_j, v_j, epoch=date_arrival)
In [20]:
(v_flypre, v_oip), = izzo.lambert(Sun.k, r_efly, r_j, date_arrival - date_flyby)
In [21]:
ss_oip = Orbit.from_vectors(Sun, r_j, v_oip, epoch=date_flyby)
In [22]:
fig, ax = plt.subplots(figsize=(9, 12))

op = OrbitPlotter(ax)
op.plot(ss_e0, label="Earth")
op.plot(ic1, label="Inner Cruise 1")
op.plot(ss_efly, osculating=False)
op.plot(ss_aph_post, label="Back to Earth")
op.plot(ss_oip, label="Jupiter Orbit Insertion Phase")
op.plot(ss_j, label="Jupiter")

ax.grid(linestyle='--')
ax.set_axisbelow(True)  # Please tone down that grid
../_images/examples_Going_to_Jupiter_with_Python_using_Jupyter_and_poliastro_22_0.png