Add get_moon function to astropy.coordinates using the erfa.plan94 - #4736

@eteq , here is another implementation for the same topic as #4724. FYI.

This needs test and change log.

@eteq I am trying to run the following test. The function is working fine for an array of times. But, since we are using the data for the EMB, the error gets magnified by a factor of (M_earth / M_moon + 1) = 82.25. Due to this, the position given by this get_moon function is quite inaccurate.

from astropy.coordinates import get_moon as gm
from astropy.time import Time

t = Time(["2016-03-27 11:30:00", "2016-03-27 17:00:00"])
moon_pos = gm(t)
print moon_pos

Output:

<SkyCoord (GCRS: obstime=['2016-03-27 11:30:00.000' '2016-03-27 17:00:00.000'],
obsgeoloc=[ 0.  0.  0.] m, 
obsgeovel=[ 0.  0.  0.] m / s): (ra, dec, distance) in (deg, deg, AU)
[(248.16744252, -20.41678509, 0.00381568),
 (250.1811754, -20.61313744, 0.00385118)]>

Thus, I am trying to find a better alternative in the erfa functions. Do you have any suggestions?

@sharadmirani - I think for the moon there is a better approach using erfa.faf03 being pursued in #4724. As I wrote on the mailing list, I think for this PR I would encourage you to change it to a get_planet function, as I think that would be very useful to have too.

@mhvk @sharadmirani Actually, I think this is currently a better approach than #4724 beause #4724 is not correct, at least as it is currently implemented.

That said, the way this particular function is implemented also seems problematic for the moon itself: this relies on the difference between fits to the Earth-Moon barycenter and the location of the earth. I'm pretty sure this is going to yield substantial errors due to machine-level errors in the subtraction, and I think that's exactly what you're running into, @sharadmirani. So I agree with @mhvk that this implementation is better for something like get_planet, which is what plan94 was actually designed for.

A better way to do get_moon would be something like what #4724 is trying to do: use ERFA functions to extract enough information about the orbit of the moon to locate it relative to the Earth, and then let the GCRS->whatever transformation do the rest of the job. The problem is that I'm pretty sure (although not certain) that the necessary orbital elements just aren't in ERFA. ERFA doesn't have a "where is the moon function", the moon stuff is there only to enable computing of libration/precession/nutation. So I think the path forward might be to extract the relevant polynomials from the IERS conventions document for the moon's orbital elements and implement them directly in python (instead of using ERFA). That will ensure the get_moon and ERFA functions are consistent, because ERFA uses the IERS standards.

cc @bmorris3, as he has a partial (or maybe complete?) version of get_moon in astroplan - not the same implementation, but still useful for comparison.

@eteq The get_moon implementation in astroplan at present uses pyephem, and therefore won't be useful to this code review.

However, I'm soon going to make a PR that uses jplephem. I need to refine the API a bit, but it's a working means of querying the DE430 JPL ephemerides. If astropy chooses to support the dependency on jplephem, then my PR will implement the moon and planets in a self-consistent way (but will it be consistent with the Sun/Earth reference frames as implemented in ERFA, I don't know).

This can be closed since #4890 is merged.

Indeed, will close this. But thanks @sharadmirani for the trial; we'll also still get a version using plan94.