# Going to Jupiter with Python using Jupyter and poliastro¶

In [1]:

import numpy as np
import matplotlib.pyplot as plt
plt.ion()

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

from poliastro.bodies import Sun, Earth, Jupiter
from poliastro.twobody import Orbit
from poliastro.maneuver import Maneuver
from poliastro.iod import izzo
from poliastro.plotting import plot, OrbitPlotter
from poliastro.util import norm

solar_system_ephemeris.set("jpl")

Out[1]:

<ScienceState solar_system_ephemeris: 'jpl'>

In [2]:

## Initial data
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 [3]:

# Initial state of the Earth
ss_e0 = Orbit.from_body_ephem(Earth, date_launch)
r_e0, v_e0 = ss_e0.rv()

/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:193: TimeScaleWarning:

Input time was converted to scale='tdb' with value 2011-08-05 16:26:06.183. Use Time(..., scale='tdb') instead.


In [4]:

r_e0

Out[4]:

$[1.0246553 \times 10^{8},~-1.023135 \times 10^{8},~-44353346] \; \mathrm{km}$
In [5]:

v_e0

Out[5]:

$[1847708.5,~1594323.4,~691089.12] \; \mathrm{\frac{km}{d}}$
In [6]:

# State of the Earth the day of the flyby
ss_efly = Orbit.from_body_ephem(Earth, date_flyby)
r_efly, v_efly = ss_efly.rv()

/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:193: TimeScaleWarning:

Input time was converted to scale='tdb' with value 2013-10-09 19:22:07.182. Use Time(..., scale='tdb') instead.


In [7]:

# 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 [8]:

# Inner Cruise 1
ic1 = ss_e0.apply_maneuver(man)
ic1.rv()

Out[8]:

(<Quantity [ 1.02465527e+08, -1.02313505e+08, -4.43533465e+07] km>,
<Quantity [2198705.82621214, 1897186.74383867,  822370.88977492] km / d>)

In [9]:

ic1.period.to(u.year)

Out[9]:

$2.1515474 \; \mathrm{yr}$
In [10]:

op = OrbitPlotter()

op.plot(ss_e0)
op.plot(ic1);

/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:434: UserWarning:

Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned


In [11]:

# We propagate until the aphelion
ss_aph = ic1.propagate(ic1.period / 2)
ss_aph.epoch

Out[11]:

<Time object: scale='tdb' format='iso' value=2012-09-01 14:40:01.690>

In [12]:

# 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[12]:

<TimeDelta object: scale='tai' format='jd' value=403.1958969055241>

In [13]:

(v_aph, v_fly), = izzo.lambert(Sun.k, ss_aph.r, ss_efly.r, time_of_flight)

In [14]:

# Check the delta-V
norm(v_aph - ss_aph.v)  # Too high!

Out[14]:

$1.079866 \; \mathrm{\frac{km}{s}}$
In [15]:

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 [16]:

op = OrbitPlotter()

op.plot(ss_e0, label="Earth")
op.plot_trajectory(ic1.sample(np.linspace(0 * u.deg, 180 * u.deg, 50)), label="Inner Cruise 1")
#op.plot(ss_efly)
op.plot(ss_aph_post, label="Back to Earth");

/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:434: UserWarning:

Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned


In [17]:

# And now, go to Jupiter!
ss_j = Orbit.from_body_ephem(Jupiter, date_arrival)
r_j, v_j = ss_j.rv()

/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:193: TimeScaleWarning:

Input time was converted to scale='tdb' with value 2016-07-05 03:19:08.184. Use Time(..., scale='tdb') instead.


In [18]:

(v_flypre, v_oip), = izzo.lambert(Sun.k, r_efly, r_j, date_arrival - date_flyby)

In [19]:

ss_oip = Orbit.from_vectors(Sun, r_j, v_oip, epoch=date_flyby)

In [20]:

ss_oip.nu

Out[20]:

$3.1126494 \; \mathrm{rad}$
In [21]:

fig, ax = plt.subplots(figsize=(9, 12))

op = OrbitPlotter(ax)

op.plot(ss_e0, label="Earth")
op.plot_trajectory(ic1.sample(np.linspace(0 * u.deg, 180 * u.deg, 50)), label="Inner Cruise 1")
#op.plot(ss_efly)
op.plot_trajectory(ss_aph_post.sample(np.linspace(180 * u.deg, 360 * u.deg + ss_junofly.nu, 50)),
label="Back to Earth");
# Plotting approximation, suggestions welcome
op.plot_trajectory(ss_oip.sample(np.linspace(ss_junofly.nu, 180 * u.deg, 50)),
label="Jupiter Orbit Insertion Phase")
op.plot(ss_j, label="Jupiter")

fig.savefig("jupiter.png");

/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:434: UserWarning:

Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned