Catch that asteroid!

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

[1]:
from astropy.utils.data import conf
conf.dataurl
[1]:
'http://data.astropy.org/'
[2]:
conf.remote_timeout
[2]:
10.0
[3]:
conf.remote_timeout = 10000

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

[4]:
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 StaticOrbitPlotter
from poliastro.plotting.misc import plot_solar_system

EPOCH = Time("2017-09-01 12:05:50", scale="tdb")
[5]:
earth = Orbit.from_body_ephem(Earth, EPOCH)
earth
[5]:
1 x 1 AU x 23.4 deg (ICRS) orbit around Sun (☉) at epoch 2017-09-01 12:05:50.000 (TDB)
[6]:
earth.plot(label=Earth);
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:1090: UserWarning:

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

../_images/examples_Catch_that_asteroid!_8_1.png
[7]:
florence = Orbit.from_sbdb("Florence")
florence
[7]:
1 x 3 AU x 22.1 deg (HeliocentricEclipticJ2000) orbit around Sun (☉) at epoch 2458600.5008007586 (TDB)

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

[8]:
florence.rv()
[8]:
(<Quantity [-2.76132873e+08, -1.71570015e+08, -1.09377634e+08] km>,
 <Quantity [13.17478676, -9.82584122, -1.48126638] km / s>)
[9]:
florence.epoch
[9]:
<Time object: scale='tdb' format='jd' value=2458600.5008007586>
[10]:
florence.epoch.iso
[10]:
'2019-04-27 00:01:09.186'
[11]:
florence.inc
[11]:
$22.142394 \; \mathrm{{}^{\circ}}$

We first propagate:

[12]:
florence = florence.propagate(EPOCH)
florence.epoch.tdb.iso
[12]:
'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:

[13]:
florence_icrs = florence.to_icrs()
florence_icrs.rv()
[13]:
(<Quantity [ 1.46404253e+08, -5.35752830e+07, -2.05656912e+07] km>,
 <Quantity [ 7.3348819 , 23.48458627, 24.12473237] km / s>)

Let us compute the distance between Florence and the Earth:

[14]:
from poliastro.util import norm
[15]:
norm(florence_icrs.r - earth.r) - Earth.R
[15]:
$6967159.9 \; \mathrm{km}$

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

[16]:
abs(((norm(florence_icrs.r - earth.r) - Earth.R) - 7060160 * u.km) / (7060160 * u.km))
[16]:
$0.013172523 \; \mathrm{}$
[17]:
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>&mdash; 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>"""
)
[17]:

And now we can plot!

[18]:
frame = plot_solar_system(outer=False, epoch=EPOCH)
frame.plot(florence_icrs, label="Florence");
/home/juanlu/Development/poliastro/poliastro-library/src/poliastro/twobody/orbit.py:1090: UserWarning:

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

../_images/examples_Catch_that_asteroid!_26_1.png

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

[19]:
frame = StaticOrbitPlotter()

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

frame.plot(florence, label="Florence (Ecliptic)")
frame.plot(florence_icrs, label="Florence (ICRS)");
../_images/examples_Catch_that_asteroid!_28_0.png

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
    (5099276.00977803, -4713719.97367401, 639109.24663496)
 (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]:
6969897 x -6974244 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:

[24]:
moon = Orbit.from_body_ephem(Moon, EPOCH)
moon
[24]:
367937 x 405209 km x 19.4 deg (GCRS) orbit around Earth (♁) at epoch 2017-09-01 12:05 (TDB)
[25]:
moon.plot(label=Moon);
../_images/examples_Catch_that_asteroid!_36_0.png

And now for the final plot:

[26]:
import matplotlib.pyplot as plt

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

anomaly outside range, clipping

../_images/examples_Catch_that_asteroid!_38_1.png

Per Python ad astra!