# poliastro - Astrodynamics in PythonÂ¶

**poliastro** is an open source (MIT) collection of Python functions useful
in Astrodynamics and Orbital Mechanics, focusing on interplanetary applications.
It provides a simple and intuitive API and handles physical quantities with
units.

View source code of poliastro!

Some of its awesome features are:

- Analytical and numerical orbit propagation
- Conversion between position and velocity vectors and classical orbital elements
- Coordinate frame transformations
- Hohmann and bielliptic maneuvers computation
- Trajectory plotting
- Initial orbit determination (Lambert problem)
- Planetary ephemerides (using SPICE kernels via Astropy)
- Computation of Near-Earth Objects (NEOs)

And more to come!

poliastro is developed by an open, international community. Release announcements and general discussion take place on our mailing list and chat.

The source code, issue tracker and wiki are hosted on GitHub, and all contributions and feedback are more than welcome. You can test poliastro in your browser using binder, a cloud Jupyter notebook server:

See benchmarks for the performance analysis of poliastro.

poliastro works on recent versions of Python and is released under the MIT license, hence allowing commercial use of the library.

```
from poliastro.examples import molniya
molniya.plot()
```

## Success storiesÂ¶

âMy team and I used Poliastro for our final project in our Summer App Space program. This module helped us in plotting asteroids by using the data provided to us. It was very challenging finding a module that can take orbits from the orbital elements, plot planets, and multiple ones. This module helped us because we were able to understand the code as most of us were beginners and make some changes the way we wanted our project to turn out. We made small changes such as taking out the axis and creating a function that will create animations. I am happy we used Poliastro because it helped us directs us in a direction where we were satisfied of our final product.â—Nayeli Ju (2017)

âWe are a group of students at University of Illinois at Urbana-Champaign, United States. We are currently working on a student AIAA/AAS satellite competition to design a satellite perform some science missions on asteroid (469219) 2016 HO3. We are using your poliastro python package in designing and visualizing the trajectory from GEO into asteroidâs orbit. Thank you for your work on poliastro, especially the APIs that are very clear and informational, which helps us significantly.â—Yufeng Luo (University of Illinois at Urbana-Champaign, United States, 2017)

âWe, at the Institute of Space and Planetary Astrophysics (ISPA, University of Karachi), are using Poliastro as part of Space Flight Dynamics Text Book development program. The idea is to develop a book suitable for undergrad students which will not only cover theoretical background but will also focus on some computational tools. We chose Poliastro as one of the packages because it was very well written and provided results with good accuracy. It is especially useful in covering some key topics like the Lambertâs problem. We support the use of Poliastro and open source software because they are easily accessible to students (without any charges, unlike some other tools). A great plus point for Poliastro is that it is Python based and Python is now becoming a very important tool in areas related to Space Sciences and Technologies.â—Prof. Jawed iqbal, Syed Faisal ur Rahman (ISPA, University of Karachi, 2016)

## ContentsÂ¶

### About poliastroÂ¶

#### OverviewÂ¶

**poliastro** is an open source collection of Python subroutines for solving
problems in Astrodynamics and Orbital Mechanics.

poliastro combines cutting edge technologies like Python JIT compiling (using numba) with young, well developed astronomy packages (like astropy and jplephem) to provide a user friendly API for solving Astrodynamics problems. It is therefore a experiment to mix the best Python open source practices with my love for Orbital Mechanics.

Since I have only solved easy academic problems I cannot assess the suitability of the library for professional environments, though I am aware that at least a company that uses it.

#### HistoryÂ¶

I started poliastro as a wrapper of some MATLAB and Fortran algorithms that I needed for a University project: having good performance was a must, so pure Python was not an option. As a three language project, it was only known to work in my computer, and I had to fight against oct2py and f2py for long hours.

Later on, I enhanced poliastro plotting capabilities to serve me in further University tasks. I removed the MATLAB (Octave) code and kept only the Fortran algorithms. Finally, when numba was mature enough, I implemented everything in pure Python and poliastro 0.3 was born.

#### Future ideasÂ¶

These are some things that I would love to implement in poliastro to expand its capabilities:

- 3D plotting of orbits
- Continuous thrust maneuvers
- Tisserand graphs
- Porkchop plots

#### Note of the original authorÂ¶

I am Juan Luis Cano RodrĂguez (two names and two surnames, itâs the Spanish
way!), an Aerospace Engineer with a passion for Astrodynamics
and the Open Source world. Before poliastro started to be a truly community
project, I started it when I was an Erasmus student
at Politecnico di Milano, an important technical university in Italy which
deeply influenced my life and ambitions and gave name to the library itself.
It is and always will be my tiny tribute to a country that will always be in
my heart and to people that never ceased to inspire me. *Grazie mille!*

### Getting startedÂ¶

#### RequirementsÂ¶

poliastro requires the following Python packages:

- NumPy, for basic numerical routines
- Astropy, for physical units and time handling
- numba (optional), for accelerating the code
- jplephem, for the planetary ephemerides using SPICE kernels
- Plotly, for interactive orbit plotting
- matplotlib, for static orbit plotting
- SciPy, for root finding and numerical propagation
- pytest, for running the tests from the package

poliastro is usually tested on Linux, Windows and OS X on Python 3.5, 3.6 and 3.7 against latest NumPy.

#### InstallationÂ¶

The easiest and fastest way to get the package up and running is to install poliastro using conda:

```
$ conda install poliastro --channel conda-forge
```

Note

We encourage users to use conda and the conda-forge packages for convenience, especially when developing on Windows.

If the installation fails for any reason, please open an issue in the issue tracker.

##### Alternative installation methodsÂ¶

If you donât want to use conda you can install poliastro from PyPI using pip:

```
$ pip install numpy # Run this one first for pip 9 and older!
$ pip install poliastro[jupyter] pytest
```

Finally, you can also install the latest development version of poliastro directly from GitHub:

```
$ pip install https://github.com/poliastro/poliastro/archive/master.zip
```

This is useful if there is some feature that you want to try, but we did not release it yet as a stable version. Although you might find some unpolished details, these development installations should work without problems. If you find any, please open an issue in the issue tracker.

Warning

It is recommended that you **never ever use sudo** with distutils, pip,
setuptools and friends in Linux because you might seriously break your
system [1][2][3][4]. Options are per user directories, virtualenv
or local installations.

##### Using poliastro on JupyterLabÂ¶

After the release of Plotly 3.0, plotting orbits using poliastro is easier than ever.

You have to install three extensions of JupyterLab to make your experience smooth:

```
$ jupyter labextension install @jupyter-widgets/jupyterlab-manager@0.38
$ jupyter labextension install plotlywidget@0.7.0
$ jupyter labextension install @jupyterlab/plotly-extension@0.18.1
```

And as the documentation of JupyterLab Extensions states:

âIn order to install JupyterLab extensions, you need to have Node.js version 4 or later installed.â

If you face any further issues, you can refer to the installation guide by Plotly.

#### TestingÂ¶

If installed correctly, the tests can be run using pytest:

```
$ python -c "import poliastro.testing; poliastro.testing.test()"
===================================== test session starts =====================================
platform linux -- Python 3.7.1, pytest-4.2.0, py-1.7.0, pluggy-0.8.1
rootdir: /home/juanlu/.miniconda36/envs/_test37/lib/python3.7/site-packages/poliastro, inifile:
collected 747 items
[...]
========= 738 passed, 3 skipped, 5 xfailed, 1 xpassed, 13 warnings in 392.12 seconds ==========
$
```

If for some reason any test fails, please report it in the issue tracker.

### User guideÂ¶

#### Defining the orbit: `Orbit`

objectsÂ¶

The core of poliastro are the `Orbit`

objects
inside the `poliastro.twobody`

module. They store all the required
information to define an orbit:

- The body acting as the central body of the orbit, for example the Earth.
- The position and velocity vectors or the orbital elements.
- The time at which the orbit is defined.

First of all, we have to import the relevant modules and classes:

```
from astropy import units as u
from poliastro.bodies import Earth, Mars, Sun
from poliastro.twobody import Orbit
```

##### From position and velocityÂ¶

There are several methods available to create
`Orbit`

objects. For example, if we have the
position and velocity vectors we can use
`from_vectors()`

:

```
# Data from Curtis, example 4.3
r = [-6045, -3490, 2500] * u.km
v = [-3.457, 6.618, 2.533] * u.km / u.s
ss = Orbit.from_vectors(Earth, r, v)
```

And thatâs it! Notice a couple of things:

Defining vectorial physical quantities using Astropy units is very easy. The list is automatically converted to a

`astropy.units.Quantity`

, which is actually a subclass of NumPy arrays.If we display the orbit we just created, we get a string with the radius of pericenter, radius of apocenter, inclination, reference frame and attractor:

>>> ss 7283 x 10293 km x 153.2 deg (GCRS) orbit around Earth (â)

If no time is specified, then a default value is assigned:

>>> ss.epoch <Time object: scale='utc' format='jyear_str' value=J2000.000> >>> ss.epoch.iso '2000-01-01 12:00:00.000'

The reference frame of the orbit will be one pseudo-inertial frame around the attractor. You can retrieve it using the

`frame`

property:>>> ss.frame <GCRS Frame (obstime=J2000.000, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s)>

Note

At the moment, there is no explicit way to set the reference system of an orbit. This is the focus of our next releases, so we will likely introduce changes in the near future. Please subscribe to this issue to receive updates in your inbox.

##### Intermezzo: quick visualization of the orbitÂ¶

If weâre working on interactive mode (for example, using the wonderful Jupyter notebook) we can immediately plot the current orbit:

```
ss.plot()
```

This plot is made in the so called *perifocal frame*, which means:

- weâre visualizing the plane of the orbit itself,
- the \((x)\) axis points to the pericenter, and
- the \((y)\) axis is turned \(90 \mathrm{^\circ}\) in the direction of the orbit.

The dotted line represents the *osculating orbit*:
the instantaneous Keplerian orbit at that point. This is relevant in the
context of perturbations, when the object shall deviate from its Keplerian
orbit.

Note

This visualization uses Plotly under the hood and works best in a Jupyter notebook.
To use the old interface based on matplotlib,
which might be more useful for batch jobs and publication-quality plots,
check out the `poliastro.plotting.static.StaticOrbitPlotter`

.

##### From classical orbital elementsÂ¶

We can also define a `Orbit`

using a set of
six parameters called orbital elements. Although there are several of
these element sets, each one with its advantages and drawbacks, right now
poliastro supports the *classical orbital elements*:

- Semimajor axis \((a)\).
- Eccentricity \((e)\).
- Inclination \((i)\).
- Right ascension of the ascending node \((\Omega)\).
- Argument of pericenter \((\omega)\).
- True anomaly \((\nu)\).

In this case, weâd use the method
`from_classical()`

:

```
# Data for Mars at J2000 from JPL HORIZONS
a = 1.523679 * u.AU
ecc = 0.093315 * u.one
inc = 1.85 * u.deg
raan = 49.562 * u.deg
argp = 286.537 * u.deg
nu = 23.33 * u.deg
ss = Orbit.from_classical(Sun, a, ecc, inc, raan, argp, nu)
```

Notice that whether we create a `Orbit`

from \((r)\) and \((v)\) or from
elements we can access many mathematical properties of the orbit:

```
>>> ss.period.to(u.day)
<Quantity 686.9713888628166 d>
>>> ss.v
<Quantity [ 1.16420211, 26.29603612, 0.52229379] km / s>
```

To see a complete list of properties, check out the
`poliastro.twobody.orbit.Orbit`

class on the API reference.

Warning

Due to limitations in the internal design of poliastro, most orbital properties are not documented. Please subscribe to this GitHub issue to receive updates about it in your inbox.

#### Moving forward in time: propagationÂ¶

Now that we have defined an orbit, we might be interested in computing
how is it going to evolve in the future. In the context of orbital
mechanics, this process is known as **propagation**, and can be
performed with the `propagate`

method of
`Orbit`

objects:

```
>>> from poliastro.examples import iss
>>> iss
6772 x 6790 km x 51.6 deg (GCRS) orbit around Earth (â)
>>> iss.epoch
<Time object: scale='utc' format='iso' value=2013-03-18 12:00:00.000>
>>> iss.nu.to(u.deg)
<Quantity 46.595804677061956 deg>
>>> iss.n.to(u.deg / u.min)
<Quantity 3.887010576192155 deg / min>
```

Using the `propagate()`

method
we can now retrieve the position of the ISS after some time:

```
>>> iss_30m = iss.propagate(30 * u.min)
>>> iss_30m.epoch # Notice we advanced the epoch!
<Time object: scale='utc' format='iso' value=2013-03-18 12:30:00.000>
>>> iss_30m.nu.to(u.deg)
<Quantity 163.1409357544868 deg>
```

For more advanced propagation options, check out the
`poliastro.twobody.propagation`

module.

#### Studying non-keplerian orbits: perturbationsÂ¶

Apart from the Keplerian propagators, poliastro also allows the user to define custom perturbation accelerations to study non Keplerian orbits, thanks to Cowellâs method:

```
>>> from poliastro.twobody.propagation import cowell
>>> from numba import njit
>>> r0 = [-2384.46, 5729.01, 3050.46] * u.km
>>> v0 = [-7.36138, -2.98997, 1.64354] * u.km / u.s
>>> initial = Orbit.from_vectors(Earth, r0, v0)
>>> @njit
... def accel(t0, state, k):
... """Constant acceleration aligned with the velocity. """
... v_vec = state[3:]
... norm_v = (v_vec * v_vec).sum() ** .5
... return 1e-5 * v_vec / norm_v
...
>>> initial.propagate(3 * u.day, method=cowell, ad=accel)
18255 x 21848 km x 28.0 deg (GCRS) orbit around Earth (â)
```

Some natural perturbations are available in poliastro to be used directly in this way. For instance, let us examine the effect of J2 perturbation:

```
>>> from poliastro.core.perturbations import J2_perturbation
>>> tof = (48.0 * u.h).to(u.s)
>>> final = initial.propagate(tof, method=cowell, ad=J2_perturbation, J2=Earth.J2.value, R=Earth.R.to(u.km).value)
```

The J2 perturbation changes the orbit parameters (from Curtis example 12.2):

```
>>> ((final.raan - initial.raan) / tof).to(u.deg / u.h)
<Quantity -0.17232668 deg / h>
>>> ((final.argp - initial.argp) / tof).to(u.deg / u.h)
<Quantity 0.28220397 deg / h>
```

For more available perturbation options, see the
`poliastro.twobody.perturbations`

module.

#### Studying artificial perturbations: thrustÂ¶

In addition to natural perturbations, poliastro also has built-in artificial perturbations (thrusts) aimed at intentional change of some orbital elements. Let us simultaineously change eccentricy and inclination:

```
>>> from poliastro.twobody.thrust import change_inc_ecc
>>> from poliastro.twobody import Orbit
>>> from poliastro.bodies import Earth
>>> from poliastro.twobody.propagation import cowell
>>> from astropy import units as u
>>> from astropy.time import Time
>>> ecc_0, ecc_f = 0.4, 0.0
>>> a = 42164
>>> inc_0, inc_f = 0.0, (20.0 * u.deg).to(u.rad).value
>>> argp = 0.0
>>> f = 2.4e-7
>>> k = Earth.k.to(u.km**3 / u.s**2).value
>>> s0 = Orbit.from_classical(Earth, a * u.km, ecc_0 * u.one, inc_0 * u.deg, 0 * u.deg, argp * u.deg, 0 * u.deg, epoch=Time(0, format='jd', scale='tdb'))
>>> a_d, _, _, t_f = change_inc_ecc(s0, ecc_f, inc_f, f)
>>> sf = s0.propagate(t_f * u.s, method=cowell, ad=a_d, rtol=1e-8)
```

The thrust changes orbit parameters as desired (within errors):

```
>>> sf.inc, sf.ecc
(<Quantity 0.34719734 rad>, <Quantity 0.00894513>)
```

For more available perturbation options, see the
`poliastro.twobody.thrust`

module.

#### Changing the orbit: `Maneuver`

objectsÂ¶

poliastro helps us define several in-plane and general out-of-plane
maneuvers with the `Maneuver`

class inside the
`poliastro.maneuver`

module.

Each `Maneuver`

consists on a list of impulses \(\Delta v_i\)
(changes in velocity) each one applied at a certain instant \(t_i\). The
simplest maneuver is a single change of velocity without delay: you can
recreate it either using the `impulse()`

method or instantiating it directly.

```
from poliastro.maneuver import Maneuver
dv = [5, 0, 0] * u.m / u.s
man = Maneuver.impulse(dv)
man = Maneuver((0 * u.s, dv)) # Equivalent
```

There are other useful methods you can use to compute common in-plane
maneuvers, notably `hohmann()`

and
`bielliptic()`

for Hohmann and
bielliptic transfers respectively. Both return the corresponding
`Maneuver`

object, which in turn you can use to calculate the total cost
in terms of velocity change \(\sum |\Delta v_i|\) and the transfer
time:

```
>>> ss_i = Orbit.circular(Earth, alt=700 * u.km)
>>> ss_i
7078 x 7078 km x 0.0 deg (GCRS) orbit around Earth (â)
>>> hoh = Maneuver.hohmann(ss_i, 36000 * u.km)
>>> hoh.get_total_cost()
<Quantity 3.6173981270031357 km / s>
>>> hoh.get_total_time()
<Quantity 15729.741535747102 s>
```

You can also retrieve the individual vectorial impulses:

```
>>> hoh.impulses[0]
(<Quantity 0 s>, <Quantity [ 0. , 2.19739818, 0. ] km / s>)
>>> hoh[0] # Equivalent
(<Quantity 0 s>, <Quantity [ 0. , 2.19739818, 0. ] km / s>)
>>> tuple(val.decompose([u.km, u.s]) for val in hoh[1])
(<Quantity 15729.741535747102 s>, <Quantity [ 0. , 1.41999995, 0. ] km / s>)
```

To actually retrieve the resulting `Orbit`

after performing a maneuver, use
the method `apply_maneuver()`

:

```
>>> ss_f = ss_i.apply_maneuver(hoh)
>>> ss_f
36000 x 36000 km x 0.0 deg (GCRS) orbit around Earth (â)
```

#### More advanced plotting: `OrbitPlotter`

objectsÂ¶

We previously saw the `poliastro.plotting.plot()`

function to easily
plot orbits. Now weâd like to plot several orbits in one graph (for example,
the maneuver we computed in the previous section). For this purpose, we
have `OrbitPlotter`

objects in the
`plotting`

module.

These objects hold the perifocal plane of the first `Orbit`

we plot in
them, projecting any further trajectories on this plane. This allows to
easily visualize in two dimensions:

```
from poliastro.plotting import OrbitPlotter
op = OrbitPlotter2D()
ss_a, ss_f = ss_i.apply_maneuver(hoh, intermediate=True)
op.plot(ss_i, label="Initial orbit")
op.plot(ss_a, label="Transfer orbit")
op.plot(ss_f, label="Final orbit")
```

Which produces this beautiful plot:

#### Where are the planets? Computing ephemeridesÂ¶

New in version 0.3.0.

Thanks to Astropy and jplephem, poliastro can now read Satellite Planet Kernel (SPK) files, part of NASAâs SPICE toolkit. This means that we can query the position and velocity of the planets of the Solar System.

The method `get_body_ephem()`

will return
a planetary orbit using low precision ephemerides available in
Astropy and an `astropy.time.Time`

:

```
from astropy import time
epoch = time.Time("2015-05-09 10:43") # UTC by default
```

And finally, retrieve the planet orbit:

```
>>> from poliastro import ephem
>>> Orbit.from_body_ephem(Earth, epoch)
1 x 1 AU x 23.4 deg (ICRS) orbit around Sun (â)
```

This does not require any external download. If on the other hand we want to use higher precision ephemerides, we can tell Astropy to do so:

```
>>> from astropy.coordinates import solar_system_ephemeris
>>> solar_system_ephemeris.set("jpl")
Downloading http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de430.bsp
|==========>-------------------------------| 23M/119M (19.54%) ETA 59s22ss23
```

This in turn will download the ephemerides files from NASA and use them for future computations. For more information, check out Astropy documentation on ephemerides.

Note

The position and velocity vectors are given with respect to the
Solar System Barycenter in the **International Celestial Reference Frame**
(ICRF), which means approximately equatorial coordinates.

#### Traveling through space: solving the Lambert problemÂ¶

The determination of an orbit given two position vectors and the time of
flight is known in celestial mechanics as **Lambertâs problem**, also
known as two point boundary value problem. This contrasts with Keplerâs
problem or propagation, which is rather an initial value problem.

The package `poliastro.iod`

allows as to solve Lambertâs problem,
provided the main attractorâs gravitational constant, the two position
vectors and the time of flight. As you can imagine, being able to compute
the positions of the planets as we saw in the previous section is the
perfect complement to this feature!

For instance, this is a simplified version of the example Going to Mars with Python using poliastro, where the orbit of the Mars Science Laboratory mission (rover Curiosity) is determined:

```
date_launch = time.Time('2011-11-26 15:02', scale='utc')
date_arrival = time.Time('2012-08-06 05:17', scale='utc')
tof = date_arrival - date_launch
ss0 = Orbit.from_body_ephem(Earth, date_launch)
ssf = Orbit.from_body_ephem(Mars, date_arrival)
from poliastro import iod
(v0, v), = iod.lambert(Sun.k, ss0.r, ssf.r, tof)
```

And these are the results:

```
>>> v0
<Quantity [-29.29150998, 14.53326521, 5.41691336] km / s>
>>> v
<Quantity [ 17.6154992 ,-10.99830723, -4.20796062] km / s>
```

#### Working with NEOsÂ¶

NEOs (Near Earth Objects) are asteroids and comets whose orbits are near to earth (obvious, isnât it?).
More correctly, their perihelion (closest approach to the Sun) is less than 1.3 astronomical units (â 200 * 10^{6} km).
Currently, they are being an important subject of study for scientists around the world, due to their status as the relatively
unchanged remains from the solar system formation process.

Because of that, a new module related to NEOs has been added to `poliastro`

as part of SOCIS 2017 project.

For the moment, it is possible to search NEOs by name (also using wildcards),
and get their orbits straight from NASA APIs, using `orbit_from_name()`

.
For example, we can get Apophis asteroid (99942 Apophis) orbit with one command, and plot it:

```
from poliastro.neos import neows
apophis_orbit = neows.orbit_from_name('apophis') # Also '99942' or '99942 apophis' works
earth_orbit = Orbit.from_body_ephem(Earth)
op = OrbitPlotter()
op.plot(earth_orbit, label='Earth')
op.plot(apophis_orbit, label='Apophis')
```

*Per Python ad astra* ;)

### Jupyter notebooksÂ¶

#### Analyzing the Parker Solar Probe flybysÂ¶

##### 1. Modulus of the exit velocity, some features of Orbit #2Â¶

First, using the data available in the reports, we try to compute some of the properties of orbit #2. This is not enough to completely define the trajectory, but will give us information later on in the process.

```
[1]:
```

```
from astropy import units as u
```

```
[2]:
```

```
T_ref = 150 * u.day
T_ref
```

```
[2]:
```

```
[3]:
```

```
from poliastro.bodies import Earth, Sun, Venus
```

```
[4]:
```

```
k = Sun.k
k
```

```
[4]:
```

```
[5]:
```

```
import numpy as np
```

```
[6]:
```

```
a_ref = np.cbrt(k * T_ref**2 / (4 * np.pi**2)).to(u.km)
a_ref.to(u.au)
```

```
[6]:
```

```
[7]:
```

```
energy_ref = (-k / (2 * a_ref)).to(u.J / u.kg)
energy_ref
```

```
[7]:
```

```
[8]:
```

```
from poliastro.twobody import Orbit
from poliastro.util import norm
from astropy.time import Time
```

```
[9]:
```

```
flyby_1_time = Time("2018-09-28", scale="tdb")
flyby_1_time
```

```
[9]:
```

```
<Time object: scale='tdb' format='iso' value=2018-09-28 00:00:00.000>
```

```
[10]:
```

```
r_mag_ref = norm(Orbit.from_body_ephem(Venus, epoch=flyby_1_time).r)
r_mag_ref.to(u.au)
```

```
[10]:
```

```
[11]:
```

```
v_mag_ref = np.sqrt(2 * k / r_mag_ref - k / a_ref)
v_mag_ref.to(u.km / u.s)
```

```
[11]:
```

##### 2. Lambert arc between #0 and #1Â¶

To compute the arrival velocity to Venus at flyby #1, we have the necessary data to solve the boundary value problem.

```
[12]:
```

```
d_launch = Time("2018-08-11", scale="tdb")
d_launch
```

```
[12]:
```

```
<Time object: scale='tdb' format='iso' value=2018-08-11 00:00:00.000>
```

```
[13]:
```

```
ss0 = Orbit.from_body_ephem(Earth, d_launch)
ss1 = Orbit.from_body_ephem(Venus, epoch=flyby_1_time)
```

```
[14]:
```

```
tof = flyby_1_time - d_launch
```

```
[15]:
```

```
from poliastro import iod
```

```
[16]:
```

```
(v0, v1_pre), = iod.lambert(Sun.k, ss0.r, ss1.r, tof.to(u.s))
```

```
[17]:
```

```
v0
```

```
[17]:
```

```
[18]:
```

```
v1_pre
```

```
[18]:
```

```
[19]:
```

```
norm(v1_pre)
```

```
[19]:
```

##### 3. Flyby #1 around VenusÂ¶

We compute a flyby using poliastro with the default value of the entry angle, just to discover that the results do not match what we expected.

```
[20]:
```

```
from poliastro.threebody.flybys import compute_flyby
```

```
[21]:
```

```
V = Orbit.from_body_ephem(Venus, epoch=flyby_1_time).v
V
```

```
[21]:
```

```
[22]:
```

```
h = 2548 * u.km
```

```
[23]:
```

```
d_flyby_1 = Venus.R + h
d_flyby_1.to(u.km)
```

```
[23]:
```

```
[24]:
```

```
V_2_v_, delta_ = compute_flyby(v1_pre, V, Venus.k, d_flyby_1)
```

```
[25]:
```

```
norm(V_2_v_)
```

```
[25]:
```

##### 4. OptimizationÂ¶

Now we will try to find the value of \(\theta\) that satisfies our requirements.

```
[26]:
```

```
def func(theta):
V_2_v, _ = compute_flyby(v1_pre, V, Venus.k, d_flyby_1, theta * u.rad)
ss_1 = Orbit.from_vectors(Sun, ss1.r, V_2_v, epoch=flyby_1_time)
return (ss_1.period - T_ref).to(u.day).value
```

There are two solutions:

```
[27]:
```

```
import matplotlib.pyplot as plt
```

```
[28]:
```

```
theta_range = np.linspace(0, 2 * np.pi)
plt.plot(theta_range, [func(theta) for theta in theta_range])
plt.axhline(0, color="k", linestyle="dashed");
```

```
[29]:
```

```
func(0)
```

```
[29]:
```

```
-9.142672330001131
```

```
[30]:
```

```
func(1)
```

```
[30]:
```

```
7.09811543934556
```

```
[31]:
```

```
from scipy.optimize import brentq
```

```
[32]:
```

```
theta_opt_a = brentq(func, 0, 1) * u.rad
theta_opt_a.to(u.deg)
```

```
[32]:
```

```
[33]:
```

```
theta_opt_b = brentq(func, 4, 5) * u.rad
theta_opt_b.to(u.deg)
```

```
[33]:
```

```
[34]:
```

```
V_2_v_a, delta_a = compute_flyby(v1_pre, V, Venus.k, d_flyby_1, theta_opt_a)
V_2_v_b, delta_b = compute_flyby(v1_pre, V, Venus.k, d_flyby_1, theta_opt_b)
```

```
[35]:
```

```
norm(V_2_v_a)
```

```
[35]:
```

```
[36]:
```

```
norm(V_2_v_b)
```

```
[36]:
```

##### 5. Exit orbitÂ¶

And finally, we compute orbit #2 and check that the period is the expected one.

```
[37]:
```

```
ss01 = Orbit.from_vectors(Sun, ss1.r, v1_pre, epoch=flyby_1_time)
ss01
```

```
[37]:
```

```
0 x 1 AU x 18.8 deg (HCRS) orbit around Sun (â) at epoch 2018-09-28 00:00:00.000 (TDB)
```

The two solutions have different inclinations, so we still have to find out which is the good one. We can do this by computing the inclination over the ecliptic - however, as the original data was in the International Celestial Reference Frame (ICRF), whose fundamental plane is parallel to the Earth equator of a reference epoch, we have change the plane to the Earth **ecliptic**, which is what the original reports use.

```
[38]:
```

```
ss_1_a = Orbit.from_vectors(Sun, ss1.r, V_2_v_a, epoch=flyby_1_time)
ss_1_a
```

```
[38]:
```

```
0 x 1 AU x 25.0 deg (HCRS) orbit around Sun (â) at epoch 2018-09-28 00:00:00.000 (TDB)
```

```
[39]:
```

```
ss_1_b = Orbit.from_vectors(Sun, ss1.r, V_2_v_b, epoch=flyby_1_time)
ss_1_b
```

```
[39]:
```

```
0 x 1 AU x 13.1 deg (HCRS) orbit around Sun (â) at epoch 2018-09-28 00:00:00.000 (TDB)
```

Letâs define a function to do that quickly for us, using the ``get_frame`

<https://docs.poliastro.space/en/latest/safe.html#poliastro.frames.get_frame>`__ function from poliastro.frames:

```
[40]:
```

```
from astropy.coordinates import CartesianRepresentation
from poliastro.frames import Planes, get_frame
def change_plane(ss_orig, plane):
"""Changes the plane of the Orbit.
"""
ss_orig_rv = ss_orig.frame.realize_frame(
ss_orig.represent_as(CartesianRepresentation)
)
dest_frame = get_frame(ss_orig.attractor, plane, obstime=ss_orig.epoch)
ss_dest_rv = ss_orig_rv.transform_to(dest_frame)
ss_dest_rv.representation_type = CartesianRepresentation
ss_dest = Orbit.from_coords(ss_orig.attractor, ss_dest_rv, plane=plane)
return ss_dest
```

```
[41]:
```

```
change_plane(ss_1_a, Planes.EARTH_ECLIPTIC)
```

```
[41]:
```

```
0 x 1 AU x 3.5 deg (HeliocentricEclipticJ2000) orbit around Sun (â) at epoch 2018-09-28 00:00:00.000 (TDB)
```

```
[42]:
```

```
change_plane(ss_1_b, Planes.EARTH_ECLIPTIC)
```

```
[42]:
```

```
0 x 1 AU x 13.1 deg (HeliocentricEclipticJ2000) orbit around Sun (â) at epoch 2018-09-28 00:00:00.000 (TDB)
```

Therefore, **the correct option is the first one**.

```
[43]:
```

```
ss_1_a.period.to(u.day)
```

```
[43]:
```

```
[44]:
```

```
ss_1_a.a
```

```
[44]:
```

And, finally, we plot the solution:

```
[45]:
```

```
from poliastro.plotting import OrbitPlotter2D
frame = OrbitPlotter2D()
frame.plot(ss0, label=Earth)
frame.plot(ss1, label=Venus)
frame.plot(ss01, label="#0 to #1")
frame.plot(ss_1_a, label="#1 to #2")
```

```
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/propagation.py:230: UserWarning:
Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned
```

#### Catch that asteroid!Â¶

```
[1]:
```

```
from astropy import units as u
from astropy.time import Time
```

```
[2]:
```

```
from astropy.utils.data import conf
conf.dataurl
```

```
[2]:
```

```
'http://data.astropy.org/'
```

```
[3]:
```

```
conf.remote_timeout
```

```
[3]:
```

```
10.0
```

First, we need to increase the timeout time to allow the download of data occur properly

```
[4]:
```

```
conf.remote_timeout = 10000
```

Then, we do the rest of the imports and create our initial orbits.

```
[5]:
```

```
from astropy.coordinates import solar_system_ephemeris
solar_system_ephemeris.set("jpl")
from poliastro.bodies import *
from poliastro.twobody import Orbit
from poliastro.plotting import OrbitPlotter2D
from poliastro.plotting.misc import plot_solar_system
EPOCH = Time("2017-09-01 12:05:50", scale="tdb")
```

```
[6]:
```

```
earth = Orbit.from_body_ephem(Earth, EPOCH)
earth
```

```
[6]:
```

```
1 x 1 AU x 23.4 deg (ICRS) orbit around Sun (â) at epoch 2017-09-01 12:05:50.000 (TDB)
```

```
[7]:
```

```
earth.plot(label=Earth)
```

```
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/propagation.py:230: UserWarning:
Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned
```

```
[8]:
```

```
from poliastro.neos import neows
```

```
[9]:
```

```
florence = neows.orbit_from_name("Florence")
florence
```

```
[9]:
```

```
1 x 3 AU x 22.1 deg (HeliocentricEclipticJ2000) orbit around Sun (â) at epoch 2458600.5 (TDB)
```

Two problems: the epoch is not the one we desire, and the inclination is with respect to the ecliptic!

```
[10]:
```

```
florence.rv()
```

```
[10]:
```

```
(<Quantity [-2.76132873e+08, -1.71570015e+08, -1.09377634e+08] km>,
<Quantity [13.17478674, -9.82584125, -1.48126639] km / s>)
```

```
[11]:
```

```
florence.epoch
```

```
[11]:
```

```
<Time object: scale='tdb' format='jd' value=2458600.5>
```

```
[12]:
```

```
florence.epoch.iso
```

```
[12]:
```

```
'2019-04-27 00:00:00.000'
```

```
[13]:
```

```
florence.inc
```

```
[13]:
```

We first propagate:

```
[14]:
```

```
florence = florence.propagate(EPOCH)
florence.epoch.tdb.iso
```

```
[14]:
```

```
'2017-09-01 12:05:50.000'
```

And now we have to convert to the same frame that the planetary ephemerides are using to make consistent comparisons, which is ICRS:

```
[15]:
```

```
florence_icrs = florence.to_icrs()
florence_icrs.rv()
```

```
[15]:
```

```
(<Quantity [ 1.46404761e+08, -5.35736589e+07, -2.05640225e+07] km>,
<Quantity [ 7.33453685, 23.48471466, 24.12478169] km / s>)
```

Let us compute the distance between Florence and the Earth:

```
[16]:
```

```
from poliastro.util import norm
```

```
[17]:
```

```
norm(florence_icrs.r - earth.r) - Earth.R
```

```
[17]:
```

This value is consistent with what ESA says! \(7\,060\,160\) km

```
[18]:
```

```
abs(((norm(florence_icrs.r - earth.r) - Earth.R) - 7060160 * u.km) / (7060160 * u.km))
```

```
[18]:
```

```
[19]:
```

```
from IPython.display import HTML
HTML(
"""<blockquote class="twitter-tweet" data-lang="en"><p lang="es" dir="ltr">La <a href="https://twitter.com/esa_es">@esa_es</a> ha preparado un resumen del asteroide <a href="https://twitter.com/hashtag/Florence?src=hash">#Florence</a> đ <a href="https://t.co/Sk1lb7Kz0j">pic.twitter.com/Sk1lb7Kz0j</a></p>— AeroPython (@AeroPython) <a href="https://twitter.com/AeroPython/status/903197147914543105">August 31, 2017</a></blockquote>
<script src="//platform.twitter.com/widgets.js" charset="utf-8"></script>"""
)
```

```
[19]:
```

La @esa_es ha preparado un resumen del asteroide #Florence đ pic.twitter.com/Sk1lb7Kz0j

— AeroPython (@AeroPython) August 31, 2017

And now we can plot!

```
[20]:
```

```
frame = plot_solar_system(outer=False, epoch=EPOCH)
frame.plot(florence_icrs, label="Florence")
```

```
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/propagation.py:230: UserWarning:
Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned
```

The difference between doing it well and doing it wrong is clearly visible:

```
[21]:
```

```
frame = OrbitPlotter2D()
frame.plot(earth, label="Earth")
frame.plot(florence, label="Florence (Ecliptic)")
frame.plot(florence_icrs, label="Florence (ICRS)")
```

And now letâs do something more complicated: express our orbit with respect to the Earth! For that, we will use GCRS, with care of setting the correct observation time:

```
[22]:
```

```
from astropy.coordinates import GCRS, CartesianRepresentation
```

```
[23]:
```

```
florence_heclip = florence.frame.realize_frame(
florence.represent_as(CartesianRepresentation)
)
```

```
[24]:
```

```
florence_gcrs_trans_cart = florence_heclip.transform_to(
GCRS(obstime=EPOCH)
).represent_as(CartesianRepresentation)
florence_gcrs_trans_cart
```

```
[24]:
```

```
<CartesianRepresentation (x, y, z) in km
(5099783.89297201, -4712095.67900917, 640778.03222918)
(has differentials w.r.t.: 's')>
```

```
[25]:
```

```
florence_hyper = Orbit.from_vectors(
Earth,
r=florence_gcrs_trans_cart.xyz,
v=florence_gcrs_trans_cart.differentials['s'].d_xyz,
epoch=EPOCH
)
florence_hyper
```

```
[25]:
```

```
6969288 x -6973635 km x 104.2 deg (GCRS) orbit around Earth (â) at epoch 2017-09-01 12:05 (TDB)
```

We now retrieve the ephemerides of the Moon, which are given directly in GCRS:

```
[26]:
```

```
moon = Orbit.from_body_ephem(Moon, EPOCH)
moon
```

```
[26]:
```

```
367937 x 405209 km x 19.4 deg (GCRS) orbit around Earth (â) at epoch 2017-09-01 12:05 (TDB)
```

```
[27]:
```

```
moon.plot(label=Moon)
```

And now for the final plot:

```
[28]:
```

```
import matplotlib.pyplot as plt
plt.ion()
from poliastro.plotting.static import StaticOrbitPlotter
```

```
[29]:
```

```
frame = StaticOrbitPlotter()
# This first plot sets the frame
frame.plot(florence_hyper, label="Florence")
# And then we add the Moon
frame.plot(moon, label=Moon)
plt.xlim(-1000000, 8000000)
plt.ylim(-5000000, 5000000)
plt.gcf().autofmt_xdate()
```

Per Python ad astra!

#### Comparing Hohmann and bielliptic transfersÂ¶

```
[1]:
```

```
import numpy as np
import matplotlib.pyplot as plt
plt.ion()
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
from astropy import units as u
from poliastro.bodies import Earth
from poliastro.twobody import Orbit
from poliastro.maneuver import Maneuver
```

```
[2]:
```

```
ZOOM = True
R = np.linspace(2, 75, num=100)
Rstar = [15.58, 40, 60, 100, 200, np.inf]
hohmann_data = np.zeros_like(R)
bielliptic_data = np.zeros((len(R), len(Rstar)))
ss_i = Orbit.circular(Earth, 1.8 * u.km)
r_i = ss_i.a
v_i = np.sqrt(ss_i.v.dot(ss_i.v))
for ii, r in enumerate(R):
r_f = r * r_i
man = Maneuver.hohmann(ss_i, r_f)
hohmann_data[ii] = (man.get_total_cost() / v_i).decompose().value
for jj, rstar in enumerate(Rstar):
r_b = rstar * r_i
man = Maneuver.bielliptic(ss_i, r_b, r_f)
bielliptic_data[ii, jj] = (man.get_total_cost() / v_i).decompose().value
idx_max = np.argmax(hohmann_data)
ylims = (0.35, 0.6)
```

```
/home/juanlu/.miniconda36/envs/poliastro37/lib/python3.7/site-packages/astropy/units/quantity.py:461: RuntimeWarning:
invalid value encountered in true_divide
```

```
[3]:
```

```
fig, ax = plt.subplots(figsize=(8, 6))
l, = ax.plot(R, hohmann_data, lw=2)
for jj in range(len(Rstar)):
ax.plot(R, bielliptic_data[:, jj], color=l.get_color())
ax.vlines([11.94, R[idx_max]], *ylims, color='0.6')
if ZOOM:
ax_zoom = zoomed_inset_axes(ax, 4, loc=4, axes_kwargs={'facecolor': '0.97'})
ax_zoom.plot(R, hohmann_data, lw=2)
for jj in range(len(Rstar)):
ax_zoom.plot(R, bielliptic_data[:, jj], color=l.get_color())
ax_zoom.vlines([11.94, R[idx_max]], *ylims, color='0.6')
ax_zoom.set_xlim(11.0, 16.0)
ax_zoom.set_ylim(0.52, 0.545)
ax_zoom.set_xticks([])
ax_zoom.set_yticks([])
ax_zoom.grid(False)
ax_zoom.set_title("4x zoom")
mark_inset(ax, ax_zoom, loc1=1, loc2=3, fc="none", ec='0.3')
ax.set_xlabel("R")
ax.set_ylabel("Relative change in velocity")
ax.set_ylim(*ylims)
ax.set_xlim(2, 75)
ax.set_title("Hohmann vs bielliptic transfers")
fig.savefig("hohmann-bielliptic-transfers.png")
```

#### Customising static orbit plotsÂ¶

The default styling for plots works pretty well however sometimes you may need to change things. The following will show you how to change the style of your plots and have different types of lines and dots

This is the default plot we will start with:

```
[1]:
```

```
from astropy.time import Time
import matplotlib.pyplot as plt
plt.ion()
from poliastro.plotting.static import StaticOrbitPlotter
from poliastro.bodies import Earth, Mars, Jupiter, Sun
from poliastro.twobody import Orbit
```

```
[2]:
```

```
epoch = Time("2018-08-17 12:05:50", scale="tdb")
plotter = StaticOrbitPlotter()
plotter.plot(Orbit.from_body_ephem(Earth, epoch), label="Earth")
plotter.plot(Orbit.from_body_ephem(Mars, epoch), label="Mars")
plotter.plot(Orbit.from_body_ephem(Jupiter, epoch), label="Jupiter");
```

```
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/propagation.py:230: UserWarning:
Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned
```

```
[3]:
```

```
epoch = Time("2018-08-17 12:05:50", scale="tdb")
plotter = StaticOrbitPlotter()
earth_plots = plotter.plot(Orbit.from_body_ephem(Earth, epoch), label=Earth)
earth_plots[0].set_linestyle("-") # solid line
earth_plots[0].set_linewidth(0.5)
earth_plots[1].set_marker("H") # Hexagon
earth_plots[1].set_markersize(15)
mars_plots = plotter.plot(Orbit.from_body_ephem(Mars, epoch), label=Mars)
jupiter_plots = plotter.plot(Orbit.from_body_ephem(Jupiter, epoch), label=Jupiter)
```

Here we get hold of the lines list from the `OrbitPlotter.plot`

method this is a list of lines. The first is the orbit line. The second is the current position marker. With the matplotlib lines objects we can start changing the style. First we make the line solid but thin line. Then we change the current position marker to a large hexagon.

More details of the style options for the markers can be found here: https://matplotlib.org/2.0.2/api/markers_api.html#module-matplotlib.markers More details of the style options on lines can be found here: https://matplotlib.org/2.0.2/api/lines_api.html However make sure that you use the set methods rather than just changing the attributes as the methods will force a re-draw of the plot.

Next we will make some changes to the other two orbits.

```
[4]:
```

```
epoch = Time("2018-08-17 12:05:50", scale="tdb")
plotter = StaticOrbitPlotter()
earth_plots = plotter.plot(Orbit.from_body_ephem(Earth, epoch), label=Earth)
earth_plots[0].set_linestyle("-") # solid line
earth_plots[0].set_linewidth(0.5)
earth_plots[1].set_marker("H") # Hexagon
earth_plots[1].set_markersize(15)
mars_plots = plotter.plot(Orbit.from_body_ephem(Mars, epoch), label=Mars)
mars_plots[0].set_dashes([0, 1, 0, 1, 1, 0])
mars_plots[0].set_linewidth(2)
mars_plots[1].set_marker("D") # Diamond
mars_plots[1].set_markersize(15)
mars_plots[1].set_fillstyle("none")
# make sure this is set if you use fillstyle 'none'
mars_plots[1].set_markeredgewidth(1)
jupiter_plots = plotter.plot(Orbit.from_body_ephem(Jupiter, epoch), label=Jupiter)
jupiter_plots[0].set_linestyle("") # No line
jupiter_plots[1].set_marker("*") # star
jupiter_plots[1].set_markersize(15)
```

You can also change the style of the plot using the matplotlib axis which can be aquired from the OrbitPlotter()

See the folling example that creates a grid, adds a title, and makes the background transparent. To make the changes clearer it goes back to the inital example.

```
[5]:
```

```
epoch = Time("2018-08-17 12:05:50", scale="tdb")
fig, ax = plt.subplots()
ax.grid(True)
ax.set_title("Earth, Mars, and Jupiter")
ax.set_facecolor("None")
plotter = StaticOrbitPlotter(ax)
plotter.plot(Orbit.from_body_ephem(Earth, epoch), label=Earth)
plotter.plot(Orbit.from_body_ephem(Mars, epoch), label=Mars)
plotter.plot(Orbit.from_body_ephem(Jupiter, epoch), label=Jupiter)
```

#### New Horizons launch and trajectoryÂ¶

Main data source: Guo & Farquhar âNew Horizons Mission Designâ http://www.boulder.swri.edu/pkb/ssr/ssr-mission-design.pdf

```
[1]:
```

```
from astropy import time
from astropy import units as u
from poliastro.bodies import Sun, Earth, Jupiter
from poliastro.twobody import Orbit
from poliastro.plotting import OrbitPlotter2D
from poliastro import iod
from poliastro.util import norm
```

##### Parking orbitÂ¶

Quoting from âNew Horizons Mission Designâ:

It was first inserted into an elliptical Earth parking orbit ofperigee altitude 165 kmandapogee altitude 215 km. [Emphasis mine]

```
[2]:
```

```
r_p = Earth.R + 165 * u.km
r_a = Earth.R + 215 * u.km
a_parking = (r_p + r_a) / 2
ecc_parking = 1 - r_p / a_parking
parking = Orbit.from_classical(
Earth,
a_parking,
ecc_parking,
0 * u.deg,
0 * u.deg,
0 * u.deg,
0 * u.deg, # We don't mind
time.Time("2006-01-19", scale="utc"),
)
print(parking.v)
parking.plot()
```

```
[0. 7.81989358 0. ] km / s
```

##### Hyperbolic exitÂ¶

Hyperbolic excess velocity:

Relation between orbital velocity \(v\), local escape velocity \(v_e\) and hyperbolic excess velocity \(v_{\infty}\):

###### Option a): Insert \(C_3\) from report, check \(v_e\) at parking perigeeÂ¶

```
[3]:
```

```
C_3_A = 157.6561 * u.km ** 2 / u.s ** 2 # Designed
a_exit = -(Earth.k / C_3_A).to(u.km)
ecc_exit = 1 - r_p / a_exit
exit = Orbit.from_classical(
Earth,
a_exit,
ecc_exit,
0 * u.deg,
0 * u.deg,
0 * u.deg,
0 * u.deg, # We don't mind
time.Time("2006-01-19", scale="utc"),
)
norm(exit.v).to(u.km / u.s)
```

```
[3]:
```

Quoting âNew Horizons Mission Designâ:

After a short coast in the parking orbit, the spacecraft was then injected into the desired heliocentric orbit by the Centaur second stage and Star 48B third stage. At the Star 48B burnout, the New Horizons spacecraft reached the highest Earth departure speed,estimated at 16.2 km/s, becoming the fastest spacecraft ever launched from Earth. [Emphasis mine]

```
[4]:
```

```
v_estimated = 16.2 * u.km / u.s
print(
"Relative error of {:.2f} %".format(
(norm(exit.v) - v_estimated) / v_estimated * 100
)
)
```

```
Relative error of 3.20 %
```

So it stays within the same order of magnitude. Which is reasonable, because real life burns are not instantaneous.

```
[5]:
```

```
from plotly.graph_objs import FigureWidget
```

```
[6]:
```

```
fig = FigureWidget()
op = OrbitPlotter2D(figure=fig)
op.plot(parking)
op.plot(exit)
fig.layout.xaxis.range = -8000, 8000
fig.layout.yaxis.range = -20000, 20000
fig
```

###### Option b): Compute \(v_{\infty}\) using the Jupyter flybyÂ¶

According to Wikipedia, the closest approach occurred at 05:43:40 UTC. We can use this data to compute the solution of the Lambert problem between the Earth and Jupiter.

```
[7]:
```

```
nh_date = time.Time("2006-01-19 19:00", scale="utc")
nh_flyby_date = time.Time("2007-02-28 05:43:40", scale="utc")
nh_tof = nh_flyby_date - nh_date
nh_earth = Orbit.from_body_ephem(Earth, nh_date)
nh_r_0, v_earth = nh_earth.rv()
nh_jup = Orbit.from_body_ephem(Jupiter, nh_flyby_date)
nh_r_f, v_jup = nh_jup.rv()
(nh_v_0, nh_v_f), = iod.lambert(Sun.k, nh_r_0, nh_r_f, nh_tof)
```

```
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:375: TimeScaleWarning:
Input time was converted to scale='tdb' with value 2006-01-19 19:01:05.184. Use Time(..., scale='tdb') instead.
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:375: TimeScaleWarning:
Input time was converted to scale='tdb' with value 2007-02-28 05:44:45.185. Use Time(..., scale='tdb') instead.
```

The hyperbolic excess velocity is measured with respect to the Earth:

```
[8]:
```

```
C_3_lambert = (norm(nh_v_0 - v_earth)).to(u.km / u.s) ** 2
C_3_lambert
```

```
[8]:
```

```
[9]:
```

```
print("Relative error of {:.2f} %".format((C_3_lambert - C_3_A) / C_3_A * 100))
```

```
Relative error of 0.51 %
```

Which again, stays within the same order of magnitude of the figure given to the Guo & Farquhar report.

##### From Earth to JupiterÂ¶

```
[10]:
```

```
nh = Orbit.from_vectors(Sun, nh_r_0.to(u.km), nh_v_0.to(u.km / u.s), nh_date)
op = OrbitPlotter2D()
op.plot(nh_jup, label=Jupiter)
op.plot(nh_earth, label=Earth)
op.plot(nh, label="New Horizons")
```

```
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/propagation.py:230: UserWarning:
Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned
```

#### 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**, initial orbit determination using the solution of the **Lambertâs problem** and **orbit plotting**.

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).

First of all, we import the necessary modules. Apart from poliastro we will make use of astropy to deal with physical units and time definitions and jplephem to compute the positions and velocities of the planets.

```
[1]:
```

```
import numpy as np
import astropy.units as u
from astropy import time
from poliastro import iod
from poliastro.bodies import Sun
from poliastro.twobody import Orbit
from poliastro.twobody.propagation import propagate
from poliastro.util import time_range
```

We need a binary file from NASA called *SPICE kernel* to compute the position and velocities of the planets. Astropy downloads it for us:

```
[2]:
```

```
from astropy.coordinates import solar_system_ephemeris, get_body_barycentric_posvel
solar_system_ephemeris.set("jpl")
```

```
[2]:
```

```
<ScienceState solar_system_ephemeris: 'jpl'>
```

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.

```
[3]:
```

```
# Initial data
N = 50
date_launch = time.Time("2011-11-26 15:02", scale="utc")
date_arrival = time.Time("2012-08-06 05:17", scale="utc")
tof = date_arrival - date_launch
tof.to(u.h)
```

```
[3]:
```

Once we have the vector of times we can use `get_body_barycentric_posvel`

to compute the array of positions and velocities of the Earth and Mars.

```
[4]:
```

```
times_vector = time_range(date_launch, end=date_arrival, periods=N)
times_vector[:5]
```

```
[4]:
```

```
<Time object: scale='utc' format='iso' value=['2011-11-26 15:02:00.000' '2011-12-01 19:14:33.082'
'2011-12-06 23:27:06.163' '2011-12-12 03:39:39.245'
'2011-12-17 07:52:12.327']>
```

```
[5]:
```

```
rr_earth, vv_earth = get_body_barycentric_posvel("earth", times_vector)
```

```
[6]:
```

```
rr_earth[:3]
```

```
[6]:
```

```
<CartesianRepresentation (x, y, z) in km
[(64600643.37167563, 1.21424866e+08, 52640047.33041222),
(52175250.21264037, 1.26254284e+08, 54733247.42732787),
(39319701.40598051, 1.30036609e+08, 56373071.6065251 )]>
```

```
[7]:
```

```
vv_earth[:3]
```

```
[7]:
```

```
<CartesianRepresentation (x, y, z) in km / d
[(-2352414.27027126, 1032013.3380897 , 447276.92493007),
(-2445842.68494247, 833043.95148986, 361105.49364196),
(-2518740.18681062, 627712.9008316 , 272197.06320273)]>
```

```
[8]:
```

```
rr_mars, vv_mars = get_body_barycentric_posvel("mars", times_vector)
```

```
[9]:
```

```
rr_mars[:3]
```

```
[9]:
```

```
<CartesianRepresentation (x, y, z) in km
[(-1.23149631e+08, 1.90752511e+08, 90809654.2669948 ),
(-1.31992428e+08, 1.86383187e+08, 89044491.25204735),
(-1.40598005e+08, 1.81677346e+08, 87118570.32883616)]>
```

```
[10]:
```

```
vv_mars[:3]
```

```
[10]:
```

```
<CartesianRepresentation (x, y, z) in km / d
[(-1730626.66251077, -811069.96095538, -325255.37513281),
(-1686163.26853493, -877100.53950512, -356742.77622963),
(-1638971.32577256, -941103.98130274, -387374.07466148)]>
```

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!

```
[11]:
```

```
# Compute the transfer orbit!
r0 = rr_earth[0].xyz
rf = rr_mars[-1].xyz
(va, vb), = iod.lambert(Sun.k, r0, rf, tof)
ss0_trans = Orbit.from_vectors(Sun, r0, va, date_launch)
ssf_trans = Orbit.from_vectors(Sun, rf, vb, date_arrival)
```

```
[12]:
```

```
from poliastro.plotting import OrbitPlotter3D
from poliastro.bodies import Earth, Mars
from plotly.graph_objs import FigureWidget
```

```
[13]:
```

```
# I like color
color_earth0 = "#3d4cd5"
color_earthf = "#525fd5"
color_mars0 = "#ec3941"
color_marsf = "#ec1f28"
color_sun = "#ffcc00"
color_orbit = "#888888"
color_trans = "#444444"
fig = FigureWidget()
frame = OrbitPlotter3D(figure=fig)
frame.set_attractor(Sun)
frame.plot_trajectory(rr_earth, label=Earth, color=color_earth0)
frame.plot_trajectory(rr_mars, label=Mars, color=color_marsf)
frame.plot_trajectory(
propagate(ss0_trans, time.TimeDelta(times_vector - ss0_trans.epoch)),
label="MSL trajectory",
color=color_trans,
)
frame.set_view(30 * u.deg, 260 * u.deg, distance=3 * u.km)
fig.layout.title = "MSL Mission: from Earth to Mars"
fig
```

This line opens a new browser tab and saves the resulting image:

```
[14]:
```

```
#frame.savefig("msl3d.png", title="MSL Mission: from Earth to Mars")
```

Not bad! Letâs celebrate with some music!

```
[15]:
```

```
from IPython.display import YouTubeVideo
YouTubeVideo('zSgiXGELjbc')
```

```
[15]:
```

#### Going to Jupiter with Python using Jupyter and poliastroÂ¶

```
[1]:
```

```
import numpy as np
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 OrbitPlotter2D
from poliastro.util import norm
solar_system_ephemeris.set("jpl")
```

```
[1]:
```

```
<ScienceState solar_system_ephemeris: 'jpl'>
```

```
[2]:
```

```
## 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")
```

```
[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:379: TimeScaleWarning:
Input time was converted to scale='tdb' with value 2011-08-05 16:26:06.183. Use Time(..., scale='tdb') instead.
```

```
[4]:
```

```
r_e0
```

```
[4]:
```

```
[5]:
```

```
v_e0
```

```
[5]:
```

```
[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:379: TimeScaleWarning:
Input time was converted to scale='tdb' with value 2013-10-09 19:22:07.182. Use Time(..., scale='tdb') instead.
```

```
[7]:
```

```
# Assume that the insertion velocity is tangential to that of the Earth
dv = C_3 ** 0.5 * v_e0 / norm(v_e0)
man = Maneuver.impulse(dv)
```

```
[8]:
```

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

```
[8]:
```

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

```
[9]:
```

```
ic1.period.to(u.year)
```

```
[9]:
```

```
[10]:
```

```
op = OrbitPlotter2D()
op.plot(ss_e0)
op.plot(ic1)
```

```
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/propagation.py:232: UserWarning:
Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned
```

```
[11]:
```

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

```
[11]:
```

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

```
[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
```

```
[12]:
```

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

```
[13]:
```

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

```
[14]:
```

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

```
[14]:
```

```
[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)
```

```
[16]:
```

```
op = OrbitPlotter2D()
op.plot(ss_e0, label="Earth")
op.plot_trajectory(ic1.sample(50, max_anomaly=180 * u.deg), label="Inner Cruise 1")
op.plot(ss_aph_post, label="Back to Earth")
```

```
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/propagation.py:232: UserWarning:
Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned
```

```
[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:379: TimeScaleWarning:
Input time was converted to scale='tdb' with value 2016-07-05 03:19:08.184. Use Time(..., scale='tdb') instead.
```

```
[18]:
```

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

```
[19]:
```

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

```
[20]:
```

```
op = OrbitPlotter2D()
# Plotting approximation, suggestions welcome
op.plot(ss_e0, label="Earth")
op.plot_trajectory(ic1.sample(50, max_anomaly=180 * u.deg), label="Inner Cruise 1")
op.plot_trajectory(
ss_aph_post.sample(50, min_anomaly=180 * u.deg, max_anomaly=400 * u.deg),
label="Back to Earth",
)
op.plot_trajectory(
ss_oip.sample(50, min_anomaly=10 * u.deg, max_anomaly=180 * u.deg),
label="Jupiter Orbit Insertion Phase",
)
op.plot(ss_j, label="Jupiter")
```

#### Plotting in 3DÂ¶

```
[1]:
```

```
import numpy as np
from poliastro.examples import *
from poliastro.plotting import *
```

```
[2]:
```

```
churi.plot(use_3d=True)
```

```
WARNING: ErfaWarning: ERFA function "taiutc" yielded 3 of "dubious year (Note 4)" [astropy._erfa.core]
```

```
[3]:
```

```
frame = OrbitPlotter3D()
frame.plot(churi)
frame.plot(Orbit.from_body_ephem(Earth))
```

```
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/propagation.py:230: UserWarning:
Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned
```

```
[4]:
```

```
frame = OrbitPlotter3D()
frame.plot(molniya)
frame.plot(Orbit.from_body_ephem(Earth))
```

```
---------------------------------------------------------------------------
NotImplementedError Traceback (most recent call last)
<ipython-input-4-c9d0aa645ff7> in <module>
2
3 frame.plot(molniya)
----> 4 frame.plot(Orbit.from_body_ephem(Earth))
~/Development/poliastro/poliastro-library/src/poliastro/plotting/_base.py in plot(self, orbit, label, color)
134 color = next(self._color_cycle)
135
--> 136 self._set_attractor(orbit.attractor)
137
138 label = generate_label(orbit, label)
~/Development/poliastro/poliastro-library/src/poliastro/plotting/_base.py in _set_attractor(self, attractor)
44 elif attractor is not self._attractor:
45 raise NotImplementedError(
---> 46 "Attractor has already been set to {}.".format(self._attractor.name)
47 )
48
NotImplementedError: Attractor has already been set to Earth.
```

```
[5]:
```

```
frame = OrbitPlotter3D()
frame.plot(molniya)
frame.plot(iss)
```

```
[6]:
```

```
from poliastro.neos import neows
from poliastro.examples import iss
eros = neows.orbit_from_name("eros")
```

```
[7]:
```

```
frame = OrbitPlotter3D()
frame.plot(Orbit.from_body_ephem(Earth), label=Earth)
frame.plot(eros, label="eros")
```

```
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/propagation.py:230: UserWarning:
Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned
```

```
[8]:
```

```
from astropy.coordinates import get_body_barycentric_posvel
from poliastro.util import time_range
```

```
[9]:
```

```
date_launch = time.Time("2011-11-26 15:02", scale="utc")
date_arrival = time.Time("2012-08-06 05:17", scale="utc")
rr_earth, _ = get_body_barycentric_posvel(
"earth", time_range(date_launch, end=date_arrival, periods=50)
)
```

```
[10]:
```

```
frame = OrbitPlotter3D()
frame.set_attractor(Sun)
frame.plot(Orbit.from_body_ephem(Earth), label=Earth)
frame.plot_trajectory(rr_earth, label=Earth)
```

```
[11]:
```

```
frame = OrbitPlotter3D()
frame.plot(eros, label="eros")
frame.plot_trajectory(rr_earth, label=Earth)
```