Replace ephemerides functions with Astropy equivalents

[astrojuanlu]

This replacement cannot be done until Astropy supports proper motions and velocities, see this bug:

astropy/astropy#4344

For instance, get_body_barycentric does not return the velocity of the planets, and therefore it's useless for our purposes.

[astrojuanlu]

First steps: barycentric velocities in Astropy astropy/astropy#5231

[astrojuanlu]

Astropy 1.3 is out 🎉

http://www.astropy.org/announcements/release-1.3.html

In particular:

• Vector arithmetic for coordinates, including taking the norm http://docs.astropy.org/en/stable/coordinates/representations.html#astropy-coordinates-representations-arithmetic
• Barycentric position and velocity of a Solar System body http://docs.astropy.org/en/stable/api/astropy.coordinates.get_body_barycentric_posvel.html

[astrojuanlu]

Using the class CartesianRepresentation is needed for vector arithmetic, but it doesn't work with normal Quantity objects for obvious reasons.

[astrojuanlu]

For some reason I'm not getting the same results with poliastro and the jpl ephemeris of Astropy. More debugging required.

In [1]: from poliastro.ephem import planet_ephem

In [2]: from astropy.time import Time

In [5]: from poliastro.ephem import EARTH

In [6]: t = Time("2014-09-22 23:22")

In [7]: planet_ephem(EARTH, t)
Out[7]: 
(<Quantity [  1.50492843e+08, -9.68589576e+05, -4.41002632e+05] km>,
 <Quantity [  -26318.45001712, 2353112.69264584, 1020089.67574404] km / d>)

In [8]: from astropy.coordinates import solar_system_ephemeris

In [9]: solar_system_ephemeris.set("jpl")
Downloading http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de430.bsp
|=====================================================================================| 119M/119M (100.00%)      1m21s
Out[9]: <ScienceState solar_system_ephemeris: 'jpl'>

In [10]: print(solar_system_ephemeris.bodies)
('sun', 'mercury', 'venus', 'earth-moon-barycenter', 'earth', 'moon', 'mars', 'jupiter', 'saturn', 'uranus', 'neptune', 'pluto')

In [11]: from astropy.coordinates import get_body_barycentric_posvel

In [12]: get_body_barycentric_posvel("earth", t)
Out[12]: 
(<CartesianRepresentation (x, y, z) in km
     (  1.50497562e+08, -967989.50119331, -440470.76265497)>,
 <CartesianRepresentation (x, y, z) in km / d
     (-26112.51286015,  2354057.67909512,  1020412.96649401)>)

In [13]: get_body_barycentric_posvel("earth-moon-barycenter", t)
Out[13]: 
(<CartesianRepresentation (x, y, z) in km
     (  1.50492823e+08, -966759.85640148, -440209.43710451)>,
 <CartesianRepresentation (x, y, z) in km / d
     (-26352.62592805,  2353112.87308378,  1020089.75407521)>)

[astrojuanlu]

Reason: it's using TDB scale in the background, see https://github.com/astropy/astropy/blob/v1.3/astropy/coordinates/solar_system.py#L218

In [28]: list(solar_system_ephemeris.kernel[0, 3].generate(t.jd1, t.jd2))
Out[28]: 
[array([  1.50492843e+08,  -9.68589576e+05,  -4.41002632e+05]),
 array([  -26318.45001712,  2353112.69264584,  1020089.67574404])]

In [29]: t.scale
Out[29]: 'utc'

In [30]: t.utc
Out[30]: <Time object: scale='utc' format='iso' value=2014-09-22 23:22:00.000>

In [31]: t.utc.jd2
Out[31]: 0.9736111111111111

In [32]: t.tdb.jd2
Out[32]: 0.9743886852002873

In [33]: list(solar_system_ephemeris.kernel[0, 3].generate(t.tdb.jd1, t.tdb.jd2))
Out[33]: 
[array([  1.50492823e+08,  -9.66759856e+05,  -4.40209437e+05]),
 array([  -26352.62592805,  2353112.87308378,  1020089.75407521])]

[astrojuanlu]

Quoting the original report:

The coordinate time scale used for DE430 and DE431 is Barycentric Dynamical Time (TDB)
as defined in terms of Barycentric Coordinate Time (TCB).

http://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf

Therefore, poliastro has been almost correct, but wrong, since #42 was merged. I'm tagging this as a bug.

[astrojuanlu]

After this change, we might want to add regression tests with strong tolerance checks. Idea: parametrize the tolerance locally so we can try a range of possible values and easily see the lowest one that fulfills the requirements.