Orbit.from_body_ephem returns wrong orbit for the Moon

[Juanlu001]

🐞 Problem

The method Orbit.from_body_ephem returns the wrong orbit for bodies not orbiting the Sun, in particular for the Moon:

In [4]: Orbit.from_body_ephem(Moon)
Out[4]: 150821918 x -150822832 km x 23.5 deg orbit around Earth (♁)

This is because from_body_ephem uses get_body_barycentric, which returns ICRS coordinates, but then sets the .parent to the Earth:

poliastro/src/poliastro/twobody/orbit.py

Lines 166 to 167 in b659a94

	r, v = get_body_barycentric_posvel(body.name, epoch)
	return cls.from_vectors(body.parent, r.xyz.to(u.km), v.xyz.to(u.km / u.day), epoch)

🎯 Goal

At least for the bodies included in poliastro.bodies, Orbit.from_body_ephem should always return a correct result.

💡 Possible solutions

Setting the .parent to the Sun in the case of the Moon would be absurd, and give incorrect results for the propagation. Therefore, we have to look for other alternatives:

• Just do not allow bodies that do not orbit the Sun. get_body_barycentric assumes bodies from the Solar System, so this would not even work for other bodies, unless a proper ephemerides is provided.
• Make a special case for the Moon, and convert from ICRS to GCRS before returning the Orbit. This would increase the number of possible frames returned by poliastro to three: ICRS for most bodies, false ecliptic for NEOS, and GCRS for the Moon. Still, it's the only correct thing to do.
• ?

In the long term, it's clear that we need strong coupling between bodies, Orbit.from_body_ephem, and reference frames (see #257, #109).

📋 Steps to solve the problem

• Comment below about what you've started working on.
• Add, commit, push your changes
• Submit a pull request and add this in comments - Addresses #<put issue number here>
• Ask for a review in comments section of pull request
• Celebrate your contribution to this project 🎉

[Juanlu001]

In the long term...

I will provide a short term fix for this next week, after the beta release.

[shreyasbapat]

http://docs.astropy.org/en/stable/api/astropy.coordinates.get_body_barycentric_posvel.html

The velocity cannot be calculated for the Moon. To just get the position, use get_body_barycentric().

Here, it is writter that we can not calculate the velocity of Moon using get_body_barycentric_posvel
Now, in that case, we can convert the ICRS to GCRS, but the problem remains the same i.e. What about the velocity.
I am researching this further.

[Juanlu001]

The documentation of Astropy is incorrect at that point: it cannot be computed with the builtin ephemeris:

In [1]: from astropy.time import Time

In [2]: from astropy.coordinates import get_body_barycentric_posvel

In [4]: get_body_barycentric_posvel("moon", Time.now())
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-4-6ace003c92bf> in <module>()
----> 1 get_body_barycentric_posvel("moon", Time.now())

~/.miniconda36/envs/py36/lib/python3.6/site-packages/astropy/coordinates/solar_system.py in get_body_barycentric_posvel(body, time, ephemeris)
    328 
    329     """
--> 330     return _get_body_barycentric_posvel(body, time, ephemeris)
    331 
    332 

~/.miniconda36/envs/py36/lib/python3.6/site-packages/astropy/coordinates/solar_system.py in _get_body_barycentric_posvel(body, time, ephemeris, get_velocity)
    223             if get_velocity:
    224                 raise KeyError("the Moon's velocity cannot be calculated with "
--> 225                                "the '{0}' ephemeris.".format(ephemeris))
    226             return calc_moon(time).cartesian
    227 

KeyError: "the Moon's velocity cannot be calculated with the 'builtin' ephemeris."

In [5]: from astropy.coordinates import solar_system_ephemeris

In [7]: solar_system_ephemeris.set("jpl")
Out[7]: <ScienceState solar_system_ephemeris: 'jpl'>

In [8]: get_body_barycentric_posvel("moon", Time.now())
Out[8]: 
(<CartesianRepresentation (x, y, z) in km
     (61039211.98626414, -1.26609874e+08, -54878422.25654943)>,
 <CartesianRepresentation (x, y, z) in km / d
     (2280846.19963135, 862485.61026938, 381518.78266656)>)

[shreyasbapat]

So, okay, only changing the ICRS to GCRS will work? As the velocity we get it also in the frame of Sun. Converting it to the frame of Earth is probably not a trivial task? Or it is? Just subtracting two vectors?

[Juanlu001]

Yes, changing to ICRS to GCRS is all that's needed. You have examples of how to do that in two notebooks:

• http://docs.poliastro.space/en/latest/examples/Catch%20that%20asteroid!.html
• http://docs.poliastro.space/en/latest/examples/Visualizing%20the%20SpaceX%20Tesla%20Roadster%20trip%20to%20Mars.html

[shreyasbapat]

Oh. This is amazing. I will do this!

[Juanlu001]

Fixed in #410.