# Catch that asteroid!¶

In [1]:

%matplotlib inline
import matplotlib.pyplot as plt

from astropy import units as u
from astropy.time import Time
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")

In [2]:

earth = Orbit.from_body_ephem(Earth, EPOCH)
earth

Out[2]:

1 x 1 AU x 23.4 deg orbit around Sun (☉)

In [3]:

plot(earth, label=Earth)

Out[3]:

[<matplotlib.lines.Line2D at 0x7f3e35e60198>,
<matplotlib.lines.Line2D at 0x7f3e35e603c8>]

In [4]:

from poliastro.neos import neows

In [5]:

florence = neows.orbit_from_name("Florence")
florence

Out[5]:

1 x 3 AU x 22.2 deg orbit around Sun (☉)


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

In [6]:

florence.epoch

Out[6]:

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

In [7]:

florence.epoch.iso

Out[7]:

'2017-09-04 00:00:00.000'

In [8]:

florence.inc

Out[8]:

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

We first propagate:

In [9]:

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

Out[9]:

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


And now we have to convert to another reference frame, using http://docs.astropy.org/en/stable/coordinates/.

In [10]:

from astropy.coordinates import (
ICRS, GCRS,
CartesianRepresentation, CartesianDifferential
)
from poliastro.frames import HeliocentricEclipticJ2000


The NASA servers give the orbital elements of the asteroids in an Heliocentric Ecliptic frame. Fortunately, it is already defined in Astropy:

In [11]:

florence_heclip = HeliocentricEclipticJ2000(
x=florence.r[0], y=florence.r[1], z=florence.r[2],
d_x=florence.v[0], d_y=florence.v[1], d_z=florence.v[2],
representation=CartesianRepresentation,
differential_cls=CartesianDifferential,
obstime=EPOCH
)
florence_heclip

Out[11]:

<HeliocentricEclipticJ2000 Coordinate (obstime=2017-09-01 12:05): (x, y, z) in km
(  1.45904366e+08, -58569290.31320047,  2270778.95771309)
(d_x, d_y, d_z) in km / s
( 7.40819577,  31.11060241,  12.80050223)>


Now we just have to convert to ICRS, which is the “standard” reference in which poliastro works:

In [12]:

florence_icrs_trans = florence_heclip.transform_to(ICRS)
florence_icrs_trans.representation = CartesianRepresentation
florence_icrs_trans

Out[12]:

<ICRS Coordinate: (x, y, z) in km
(  1.46271269e+08, -53880006.88710973, -20906928.0521954)
(v_x, v_y, v_z) in km / s
( 7.39978737,  23.46064313,  24.1234135)>

In [13]:

florence_icrs = Orbit.from_vectors(
Sun,
r=[florence_icrs_trans.x, florence_icrs_trans.y, florence_icrs_trans.z] * u.km,
v=[florence_icrs_trans.v_x, florence_icrs_trans.v_y, florence_icrs_trans.v_z] * (u.km / u.s),
epoch=florence.epoch
)
florence_icrs

Out[13]:

1 x 3 AU x 44.6 deg orbit around Sun (☉)

In [14]:

florence_icrs.rv()

Out[14]:

(<Quantity [  1.46271269e+08, -5.38800069e+07, -2.09069281e+07] km>,
<Quantity [  7.39978737, 23.46064313, 24.1234135 ] km / s>)


Let us compute the distance between Florence and the Earth:

In [15]:

from poliastro.util import norm

In [16]:

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

Out[16]:

$7060313.3 \; \mathrm{km}$

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

In [17]:

from IPython.display import HTML

HTML(
)

Out[17]:


And now we can plot!

In [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")

Out[18]:

[<matplotlib.lines.Line2D at 0x7f3e2bf66320>,
<matplotlib.lines.Line2D at 0x7f3e2bf66518>]


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

In [19]:

frame = OrbitPlotter()

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

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

Out[19]:

[<matplotlib.lines.Line2D at 0x7f3e2be71048>,
<matplotlib.lines.Line2D at 0x7f3e2be71240>]


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:

In [20]:

florence_gcrs_trans = florence_heclip.transform_to(GCRS(obstime=EPOCH))
florence_gcrs_trans.representation = CartesianRepresentation
florence_gcrs_trans

Out[20]:

<GCRS Coordinate (obstime=2017-09-01 12:05, obsgeoloc=( 0.,  0.,  0.) m, obsgeovel=( 0.,  0.,  0.) m / s): (x, y, z) in km
( 4966319.35958239, -5018473.35356456,  297867.61376881)
(v_x, v_y, v_z) in km / s
(-2.76873111, -1.96008601,  13.10279932)>

In [21]:

florence_hyper = Orbit.from_vectors(
Earth,
r=[florence_gcrs_trans.x, florence_gcrs_trans.y, florence_gcrs_trans.z] * u.km,
v=[florence_gcrs_trans.v_x, florence_gcrs_trans.v_y, florence_gcrs_trans.v_z] * (u.km / u.s),
epoch=EPOCH
)
florence_hyper

Out[21]:

7066691 x -7071046 km x 104.3 deg orbit around Earth (♁)


Notice that the ephemerides of the Moon is also given in ICRS, and therefore yields a weird hyperbolic orbit!

In [22]:

moon = Orbit.from_body_ephem(Moon, EPOCH)
moon

Out[22]:

151218466 x -151219347 km x 23.3 deg orbit around Earth (♁)

In [23]:

moon.a

Out[23]:

$-440.42131 \; \mathrm{km}$
In [24]:

moon.ecc

Out[24]:

$343350.57 \; \mathrm{}$

So we have to convert again.

In [25]:

moon_icrs = ICRS(
x=moon.r[0], y=moon.r[1], z=moon.r[2],
v_x=moon.v[0], v_y=moon.v[1], v_z=moon.v[2],
representation=CartesianRepresentation,
differential_cls=CartesianDifferential
)
moon_icrs

Out[25]:

<ICRS Coordinate: (x, y, z) in km
(  1.41399531e+08, -49228391.42507221, -21337616.62766309)
(v_x, v_y, v_z) in km / s
( 11.10890252,  25.6785744,  11.0567569)>

In [26]:

moon_gcrs = moon_icrs.transform_to(GCRS(obstime=EPOCH))
moon_gcrs.representation = CartesianRepresentation
moon_gcrs

Out[26]:

<GCRS Coordinate (obstime=2017-09-01 12:05, obsgeoloc=( 0.,  0.,  0.) m, obsgeovel=( 0.,  0.,  0.) m / s): (x, y, z) in km
( 94189.90120828, -367278.24304992, -133087.21297573)
(v_x, v_y, v_z) in km / s
( 0.94073662,  0.25786326,  0.03569047)>

In [27]:

moon = Orbit.from_vectors(
Earth,
[moon_gcrs.x, moon_gcrs.y, moon_gcrs.z] * u.km,
[moon_gcrs.v_x, moon_gcrs.v_y, moon_gcrs.v_z] * (u.km / u.s),
epoch=EPOCH
)
moon

Out[27]:

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


And finally, we plot the Moon:

In [28]:

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


And now for the final plot:

In [29]:

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