# Catch that asteroid!¶

[1]:

import matplotlib.pyplot as plt
plt.ion()

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 OrbitPlotter, plot

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 (☉)

[7]:

plot(earth, label=Earth);

/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:475: 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 (☉)


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

[10]:

florence.epoch

[10]:

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

[11]:

florence.epoch.iso

[11]:

'2018-03-23 00:00:00.000'

[12]:

florence.inc

[12]:

$22.144811 \; \mathrm{{}^{\circ}}$

We first propagate:

[13]:

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

[13]:

'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:

[14]:

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

[14]:

(<Quantity [ 1.46265478e+08, -5.38817374e+07, -2.08986005e+07] km>,
<Quantity [ 7.3998822 , 23.46299461, 24.12028277] km / s>)


Let us compute the distance between Florence and the Earth:

[15]:

from poliastro.util import norm

[16]:

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

[16]:

$7057830.8 \; \mathrm{km}$

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

[17]:

from IPython.display import HTML

HTML(
)

[17]:


And now we can plot!

[18]:

frame = OrbitPlotter()

frame.plot(earth, label="Earth")

frame.plot(Orbit.from_body_ephem(Mars, EPOCH))
frame.plot(Orbit.from_body_ephem(Venus, EPOCH))
frame.plot(Orbit.from_body_ephem(Mercury, EPOCH))

frame.plot(florence_icrs, label="Florence");

/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:475: 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:

[19]:

frame = OrbitPlotter()

frame.plot(earth, label="Earth")

frame.plot(florence, label="Florence (Ecliptic)")
frame.plot(florence_icrs, label="Florence (ICRS)");

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

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



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:

[20]:

from astropy.coordinates import GCRS, CartesianRepresentation

[21]:

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

[22]:

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

[22]:

<CartesianRepresentation (x, y, z) in km
(4960528.40227817, -5020204.24301458, 306195.40673516)
(has differentials w.r.t.: 's')>

[23]:

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

[23]:

7064205 x -7068561 km x 104.3 deg (GCRS) orbit around Earth (♁)


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

[24]:

moon = Orbit.from_body_ephem(Moon, EPOCH)
moon

[24]:

367937 x 405209 km x 19.4 deg (GCRS) orbit around Earth (♁)

[25]:

plot(moon, label=Moon)
plt.gcf().autofmt_xdate()


And now for the final plot:

[26]:

frame = OrbitPlotter()

# 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()