# Going to Mars with Python using poliastro¶

This is an example on how to use poliastro, a little library I've been working on to use in my Astrodynamics lessons. It features conversion between classical orbital elements and position vectors, propagation of Keplerian orbits, perturbation equations, initial orbit determination using the solution of the Lambert's problem and medium-precision ephemerides.

In this example we're going to draw the trajectory of the mission Mars Science Laboratory (MSL), which carried the rover Curiosity to the surface of Mars in a period of something less than 9 months.

Note: This is a very simplistic analysis which doesn't take into account many important factors of the mission, but can serve as an starting point for more serious computations (and as a side effect produces a beautiful plot at the end).

In [1]:
%pylab inline


Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].



First of all, we are importing the necessary modules. We used the latest version of poliastro from master on GitHub (88dff731d02) because version 0.1.0 has several flaws not yet solved in a 0.2.0.

In [2]:
from datetime import datetime

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from poliastro import ephem, iod, twobody, util
from poliastro import plotting
from poliastro.constants import k_Sun


The initial data was gathered from Wikipedia: the date of the launch was on November 26, 2011 at 15:02 UTC and landing was on August 6, 2012 at 05:17 UTC. We compute then the time of flight, which is exactly what it sounds. It is a crucial parameter of the mission.

In [3]:
# Initial data
N = 50

date_launch = datetime(2011, 11, 26, 15, 2)
date_arrival = datetime(2012, 8, 6, 5, 17)
tof = (date_arrival - date_launch).total_seconds()


With the date of launch and the date of landing we can compute the Julian days. The Julian day is an integer assigned to a date, and it's useful for not having to deal with leap years, changes of calendar and other messy stuff. It is measured from around 4713 BC so it is a pretty big number, as we'll see:

In [4]:
# Calculate Julian days
jd_launch = ephem.jd(date_launch)
jd_arrival = ephem.jd(date_arrival)
print(jd_launch, jd_arrival)

# Build a vector of N Julian days between the launch and the landing
jd_vec = np.linspace(jd_launch, jd_arrival, num=N)

2455892.12639 2456145.72014



Once we have the Julian day, we can compute the orbital elements of the Earth and Mars. And what the heck are those? They are a set of six parameters which geometrically define the orbit in space, as you can see in this excellent picture from Wikipedia:

https://en.wikipedia.org/wiki/File:Orbit1.svg

With the function mean_elements we can compute them. As we're giving it a vector of Julian days, it will compute a N x 6 matrix with the evolution of the six elements in time.

In [5]:
# Compute position of the planets in time
coe_earth = ephem.mean_elements(jd_vec, ephem.EARTH)
coe_mars = ephem.mean_elements(jd_vec, ephem.MARS)

# These are the elements of the Earth at the time of launch
print("""Semi-major axis (a): {:.0f} km
Eccentricity (e): {:.5f}
Longitude of ascending node (Ω): {:.3f} rad
Argument of the pericenter (ω): {:.3f} rad
""".format(*coe_earth[0]))

Semi-major axis (a): 149598023 km
Eccentricity (e): 0.01670
Longitude of ascending node (Ω): 0.000 rad
Argument of the pericenter (ω): 1.800 rad



And what's the point of having these six arcane names here? With them we can compute the actual position and velocity of the planets, using the function coe2rv (classical orbital elements to r and v). We have to provide the standard gravitational parameter of the Sun, which is a measure of the force with which it attracts the planets.

In [6]:
# Compute position and velocity in time
rr_earth, vv_earth = twobody.coe2rv(k_Sun, *coe_earth.T)
rr_mars, vv_mars = twobody.coe2rv(k_Sun, *coe_mars.T)


Now we have the history of position and velocity of the Earth and Mars. To compute the transfer orbit, we have the useful function lambert: according to a theorem with the same name, the transfer orbit between two points in space only depends on those two points and the time it takes to go from one to the other. We have the starting and final position and we have the time of flight: there we go!

In [7]:
# Compute the transfer orbit!
r0 = rr_earth[0]
rf = rr_mars[-1]

va, vb = iod.lambert(k_Sun, r0, rf, tof)  # Lots of math happening here

coe0_trans = twobody.rv2coe(k_Sun, r0, va)
coef_trans = twobody.rv2coe(k_Sun, rf, vb)


What just happened is that lambert computed the velocity in both the starting and the final point, completely defining the orbit, and we converted it back to classical orbital elements using the function coe2rv. Once we have them, we are going to use a little trick with the coe2rv functions to retrieve the part of the orbits not traversed by the planets. For that, we pass coe2rv all of the elements except from the true anomaly, and we reverse the values:

In [8]:
# Extract whole orbit of Earth, Mars and transfer (for plotting)
# New syntax on Python 3: extended iterable unpacking. Very cool!
# http://www.python.org/dev/peps/pep-3132/
*_, nu0_trans = coe0_trans
*_, nuf_trans = coef_trans
rr_trans, _ = twobody.coe2rv(k_Sun, *coe0_trans[:-1], nu=np.linspace(*util.direct_angles(nu0_trans, nuf_trans), num=N))

*_, nu0_earth = coe_earth[0]
*_, nuf_earth = coe_earth[-1]
*_, nu0_mars = coe_mars[0]
*_, nuf_mars = coe_mars[-1]
rr_earth_rest, _ = twobody.coe2rv(k_Sun, *coe_earth.T[:-1], nu=np.linspace(*util.direct_angles(nuf_earth, nu0_earth), num=N))
rr_mars_rest, _ = twobody.coe2rv(k_Sun, *coe_mars.T[:-1], nu=np.linspace(*util.direct_angles(nuf_mars, nu0_mars), num=N))


And finally, we are in conditions to plot the figure! There is no more magic here, just passing the position vectors to matplotlib plot function and adding some style to the plot. Two views are offered.

In [9]:
# Plot figure
# https://github.com/matplotlib/matplotlib/blob/master/lib/matplotlib/streamplot.py#L140
fig = plt.figure(figsize=(10, 10))

def plot_body(ax, r, color, size, border=False, **kwargs):
"""Plots body in axes object.

"""
return ax.plot(*r[:, None], marker='o', color=color, ms=size, mew=int(border), **kwargs)

# I like color
color_earth0 = '#3d4cd5'
color_earthf = '#525fd5'
color_mars0 = '#ec3941'
color_marsf = '#ec1f28'
color_sun = '#ffcc00'
color_orbit = '#888888'
color_trans = '#444444'

# Plotting orbits is easy!
ax.plot(*rr_earth.T, c=color_earth0)
ax.plot(*rr_mars.T, c=color_mars0)
ax.plot(*rr_trans.T, c=color_trans)

ax.plot(*rr_earth_rest.T, ls='--', c=color_orbit)
ax.plot(*rr_mars_rest.T, ls='--', c=color_orbit)

# But plotting planets feels even magical!
plot_body(ax, np.zeros(3), color_sun, 16)

plot_body(ax, r0, color_earth0, 8)
plot_body(ax, rr_earth[-1], color_earthf, 8)

plot_body(ax, rr_mars[0], color_mars0, 8)
plot_body(ax, rf, color_marsf, 8)

ax.text(-0.75e8, -3.5e8, -1e8, "MSL mission:\nfrom Earth to Mars", size=20, ha='center', va='center', bbox={"pad": 30, "lw": 0, "fc": "w"})
ax.text(r0[0] * 1.15, r0[1] * 1.25, r0[2] * 10, "Earth at launch\n(26 Nov)", ha="left", va="bottom", backgroundcolor='#ffffff')
ax.text(rf[0] * 0.7, rf[1] * 1.1, rf[2], "Mars at arrival\n(6 Ago)", ha="left", va="top")
ax.text(-1.9e8, 8e7, 0, "Transfer\norbit", ha="right", va="center")

# Tune axes
ax.set_xlim(-3e8, 3e8)
ax.set_ylim(-3e8, 3e8)
ax.set_zlim(-3e8, 3e8)

# And finally!
ax.view_init(30, 260)
fig.savefig("trans_30_260.png", bbox_inches='tight')

In [10]:
ax.view_init(80, 260)
display(fig)


Not bad! Hope you found it interesting. In case you didn't but are still reading, here is some music that you may enjoy:

In [11]:
from IPython.display import YouTubeVideo