Solar system plan94 - #4920

p.s. The tests were simplified somewhat but two classes still contain the same function twice, so it can obviously be improved. Also needs some tests for sun and `earth'.

@eteq, @StuartLittlefair - I added some tests for the Sun. While doing that, I found a strange discrepancy with get_sun. Might there be something wrong with that routine, or is there a difference in GCRS that I do not understand? (for reference, the new tests check consistency of get_body('sun') with horizons, so it seems more likely a problem with get_sun...)

In [3]: get_body('sun', t)
Out[3]: 
<SkyCoord (GCRS: obstime=2016-09-25 00:00:00.000, obsgeoloc=[ 0.  0.  0.] m, obsgeovel=[ 0.  0.  0.] m / s): (ra, dec, distance) in (deg, deg, AU)
    (181.94539939, -0.84295538, 1.00285985)>

In [4]: get_body('sun', t, ephemeris='de430')
Out[4]: 
<SkyCoord (GCRS: obstime=2016-09-25 00:00:00.000, obsgeoloc=[ 0.  0.  0.] m, obsgeovel=[ 0.  0.  0.] m / s): (ra, dec, distance) in (deg, deg, km)
    (181.94540454, -0.84296274, 150025696.65830812)>

In [5]: get_sun(t)
Out[5]: 
<SkyCoord (GCRS: obstime=2016-09-25 00:00:00.000, obsgeoloc=[ 0.  0.  0.] m, obsgeovel=[ 0.  0.  0.] m / s): (ra, dec, distance) in (deg, deg, AU)
    (182.85249255, -1.23568342, 1.00285986)>

EDIT: the following does again give a consistent result. Maybe get_sun should use this?

SkyCoord(HeliocentricTrueEcliptic(CartesianRepresentation([0., 0., 0.]*u.km), equinox=t - 1.00285985 * u.au / c.c)).transform_to(GCRS(obstime=t))Out[20]: 
<SkyCoord (GCRS: obstime=2016-09-25 00:00:00.000, obsgeoloc=[ 0.  0.  0.] m, obsgeovel=[ 0.  0.  0.] m / s): (ra, dec, distance) in (deg, deg, km)
    (181.94539939, -0.84295538, 150025698.42827976)>

@StuartLittlefair, @eteq - the comment about inconsistencies in SOLAR_SYSTEM_BODIES reminds me that I should have been explicit that this was put up perhaps in a bit too much of a hurry, but I wanted to be sure there was something to look at, as I really hope to also have support for plan94 in 1.2.0, and be sure that users explicitly ask for downloading a JPL ephemeris.

@mhvk, @eteq I agree that something supporting plan94 is a good idea for 1.2.0. Apart from consistency, is the primary motivation avoiding the large download? I do wonder if we should make de432 the default. OK, it's a restricted date range but on the other hand it's only a 10MB download and lasts until 2050...

@mhvk -sorry this is a bit out of order, I can't see how to reply to a comment directly.

The discrepancy between the different methods for calculating the Sun's position is worrying. Given that the tests against Horizons pass I'm thinking that either the aberration calculation is wrong (I didn't look too closely). Alternatively, the transform between ICRS and GCRS also includes deflection by large bodies in the Solar system. I don't suppose Jupiter is getting in the way at the time you've chosen is it?

@StuartLittlefair - on get_sun: I just tried getting the Sun using HeliocentricTrueEcliptic(0.) (see edited comment above) and this does give a result consistent with Horizons. I think this would be a pretty logical way to implement get_sun. It might also be the more logical way to implement plan94 planet positions...

@mhvk works well for get_sun, but for the planets you'd have to convert the ICRS planet positions into ecliptic coordinates which might be more work than it's worth.

@StuartLittlefair - actually, plan94 (EDIT: annoying automatic return) gives heliocentric positions, though I am still unsure exactly which frame it is. From erfa/plan94.c:

**  4) On successful return, the array pv contains the following:
**
**        pv[0][0]   x      }
**        pv[0][1]   y      } heliocentric position, AU
**        pv[0][2]   z      }
**
**        pv[1][0]   xdot   }
**        pv[1][1]   ydot   } heliocentric velocity, AU/d
**        pv[1][2]   zdot   }
**
**     The reference frame is equatorial and is with respect to the
**     mean equator and equinox of epoch J2000.0.

p.s. And epv00 gives also the heliocentric position of the Earth, so in principle one could go from there.

They're in the FK5(equinox=j2000), so almost ICRS (to about 80mas).

So although you've got heliocentric coordinates it's not easy to get the coordinates to make a HeliocentricTrueEcliptic frame from them.

You could use the coordinates from plan94 directly in a HCRS frame...

Actually, the transformation to equatorial is done at the very end of plan94 and easy to undo. But I think it is probably easiest to just keep it as I have it, i.e., tranform to barycentric using the epv00 output.

Sorry: I did not read right: yes, I had thought about HCRS yesterday, but for some reason the conversion away from HCRS does not seem possible. Should have raised an issue about this earlier. See #4922.

@mhvk - Hmm, that get_sun discrepency is odd - I'm seeing the same. Out of curiosity, did you try removing the abberation trick to see if that changes things?

Another possibility to consider: does HORIZONS include the aberration correction? It might be that their output doesn't include that for whatever reason, which would mean we need to add it back in to alll the ephemerides calculations.

And there's also the light deflection question: I can't remember off hand if there's an ICRS->GCRS transform here, but if so, it might be that we're getting the sun's "deflection"? If so, that might mean we need to add some logic into the ICRS->GCRS transform to decide when to use deflection...

Oh, and I like the idea of using ScienceState for this, but in that case it should be easier to specify jpl ephemerides. Maybe something like the kernel_url, but with more convenient strings? I.e. emphemerides_type, which would support:

• None: the default, which raises an exception that helpfully explains how to choose a kernel
• erfa/`plan94': this would be the "no-download" option
• jpl: this would be an alias for de430, which I think is the best default JPL because it's pretty generically usable without being an insane download (431 is ~4GB).
• de430/de432: pretty obvious
• de431: this is a bit more awkward to support because it comes in two parts from JPL, which would need to be downloaded and joined on-the-fly. But I imagine we can special-case it without too much trouble.
• Anything else would be treated as a URL

OK, back in. I'll go ahead and rebase. My general feeling is to keep ephemeris as a keyword, but make the default be whatever science state one has defined. If it is set, it is taken instead of the science state, but does not override it.

Also, I agree the state should be something explicit, like a string, so I'll implement that.

Finally, to me it makes most sense to actually download the ephemeris when setting the state, as this is the point the user might expect a delay, and not wait until an attempt is made to use the ephemeris. This has the advantage that validation can be done at the setting stage as well.

I got some way towards a science state code that I thought was rather elegant, but students happened so I didn't finish it. Also, I have social commitments this evening, so I won't be able to contribute much. However, I'll post what I had as a comment here. Maybe it will be useful.

The idea is that the code would be used as:

import astropy.coordinates as coord
from astropy.time import Time

with coord.solar_system_ephemeris.set('plan94'):
    print('Allowed bodies ', coord.solar_system_ephemeris.bodies)
    pos = coord.get_body_barycentric('venus', Time.now())

Anyway, here it is....

__all__ = ["solar_system_ephemeris","get_body_barycentric"]
"""
each value in the BODIES dictionary a list of kernel pairs needed
to find the barycentric position of that object from a JPL kernel.

Every JPL 'planet' ephemeris kernel can be assumed to have these 
kernel pairs contained within, so these pairings should work with
all SPK files from JPL.
"""
BODY_NAME_TO_KERNEL_SPEC = OrderedDict(
                                      (('sun', [(0, 10)]),
                                       ('mercury', [(0, 1), (1, 199)]),
                                       ('venus', [(0, 2), (2, 299)]),
                                       ('earth-moon-barycenter', [(0, 3)]),
                                       ('earth',  [(0, 3), (3, 399)]),
                                       ('moon', [(0, 3), (3, 301)]),
                                       ('mars', [(0, 4)]),
                                       ('jupiter', [(0, 5)]),
                                       ('saturn', [(0, 6)]),
                                       ('uranus', [(0, 7)]),
                                       ('neptune', [(0, 8)]),
                                       ('pluto', [(0, 9)]))
                                      )

"""
The ERFA `plan94` routines have a limited subset of objects.
"""                                  
ERFA_BODY_NAME_TO_PLANET_INDEX = OrderedDict(
    (('mercury', 1),
     ('venus', 2),
     ('earth-moon-barycenter', 3),
     ('mars', 4),
     ('jupiter', 5),
     ('saturn', 6),
     ('uranus', 7),
     ('neptune', 8)))

SOLAR_SYSTEM_BODIES = None

class classproperty(object):
    def __init__(self, getter):
        self.getter = getter
    def __get__(self, instance, owner):
        return self.getter(owner)    

class solar_system_ephemeris(ScienceState):
    """
    The ephemeris to use for calculating solar system body positions.

    The default is `None`, in which case any attempt to calculate solar system
    body positions will raise an Exception.

    Valid choices (other than `None`) are:
        plan94 : Uses the ERFA routines to find approximate positions of a limited
                 subset of bodies. Accuracy w.r.t JPL ephemeris is about 1 arcminute.
                 No moon, or earth positions available.
        jpl : Uses the NASA JPL Satellite Planet Kernel (SPK) file DE430. This is ~120 MB
              and will be downloaded the first time a position is calculated and 
              cached from then on. It covers the date range ~1550-2650 CE.
        de430 : An alias for jpl.
        de432 : Uses the DE432 SPK file instead. This in only 10MB, but covers a limited
                date range of ~1950-2050 CE.
        url : You can also supply the URL to any valid SPK file.
    """

    _value = None

    @classmethod
    def validate(cls, value):
        global KERNEL
        if value is None:
            KERNEL = None
            return value
        if value.lower() not in ['jpl', 'de430', 'de432', 'plan94']:
            # CUSTOM URL
            try:
                six.moves.urllib.parse.urlparse(value)
            except:
                # this is not a valid url or one of the approved items                
                raise ValueError('{} is not a valid URL, or one of the pre-defined '
                                 'ephemerides'.format(value))
        KERNEL = None
        return value

    @classproperty    
    def url(cls):
        if cls._value is None or cls._value.lower() == 'plan94':
            url = None
        elif cls._value in ['jpl', 'de430', 'de432']:
            if cls._value.lower() == 'jpl':
                de_string = 'de430'
            else:
                de_string = cls._value.lower()
            de_number = int(de_string[2:])
            url = ('http://naif.jpl.nasa.gov/pub/naif/'
                   'generic_kernels/spk/planets/de{0:03d}.bsp'.format(de_number))
        else:
            url = cls._value
        return url     

    @classproperty
    def bodies(cls):
        if cls._value is None:
            return None
        if cls._value.lower() == 'plan94':
            return tuple(ERFA_BODY_NAME_TO_PLANET_INDEX.keys())
        else:
            return tuple(BODY_NAME_TO_KERNEL_SPEC.keys())     


def _get_kernel():
    """
    Try importing jplephem, download/retrieve from cache the Satellite Planet
    Kernel corresponding to the given ephemeris.
    """
    global KERNEL
    if KERNEL is None and solar_system_ephemeris.url is not None:
        try:
            from jplephem.spk import SPK
    except ImportError:
        raise ImportError("Solar system JPL ephemeris calculations require "
                          "the jplephem package "
                          "(https://pypi.python.org/pypi/jplephem)")

    KERNEL = SPK.open(_download_spk_file(solar_system_ephemeris.url))
    return KERNEL


def get_body_barycentric(body, time):
    """
    Calculate the barycentric position of the solar system body ``body``
    in cartesian coordinates.
    Parameters
    ----------
    body : str
        The solar system body to calculate. The allowed values for ``body``
        depend on the ephemeris, but standard ones can be found in
        ``astropy.coordinates.SOLAR_SYSTEM_BODIES``.
    time : `~astropy.time.Time`
        Time at which the position should be calculated.
    Returns
    -------
    cartesian_position : `~astropy.coordinates.CartesianRepresentation`
        Barycentric (ICRS) position of the body in cartesian coordinates
    """
    if solar_system_ephemeris.get() is None:
        raise ValueError("""You must set the solar_system_ephemeris first. For example:
                                >>> from astropy import coordinates as coord
                                >>> coord.solar_system_ephemeris.set('plan94')""")

    jd1, jd2 = get_jd12(time, 'tdb')
    body = body.lower()

    if solar_system_ephemeris.url is None:
        earth_pv_helio, earth_pv_bary = erfa.epv00(jd1, jd2)
        if body == 'earth':
            cartesian_position_body = earth_pv_bary[..., 0, :]
        else:
            sun_bary = earth_pv_bary[..., 0, :] - earth_pv_helio[..., 0, :]
            if body == 'sun':
                cartesian_position_body = sun_bary
            else:
                body_index = ERFA_BODY_NAME_TO_PLANET_INDEX[body]
                body_pv_helio = erfa.plan94(jd1, jd2, body_index)
                cartesian_position_body = body_pv_helio[..., 0, :] + sun_bary

        barycen_to_body_vector = u.Quantity(
            np.rollaxis(cartesian_position_body, -1, 0), u.au)

    else:
        # Lookup chain for JPL ephemeris.
        chain = BODY_NAME_TO_KERNEL_SPEC[body.lower()]

        kernel = _get_kernel()

        cartesian_position_body = sum([kernel[pair].compute(jd1, jd2)
                                       for pair in chain])

        barycen_to_body_vector = u.Quantity(cartesian_position_body, unit=u.km)

    return CartesianRepresentation(barycen_to_body_vector)

OK, finally a complete rework. Please have a look... The remote data test may currently fail since the times used in the tests are sufficiently nearby that IERS A is downloaded even for the erfa ephemeris test (see #4933; I may have to change the "future Kitt Peak" test too, since once that is no longer future and better UT1-UTC become available, the positions will change).

@mhvk Is the motivation behind allowing the ephemeris as an argument to get_body etc to save typing for one-off calculations? If so, I wonder if we are better off making plan94 the default ephemeris instead?

My reasoning is that I would expect setting a science state to affect all calls to get_body, get_body_barycentric etc, even if I didn't personally write the code that calls them. Suppose @bmorris3 hard-codes an ephemeris into some calls in astroplan. If I set the science state at the start of a notebook, say, those calls in astroplan would ignore my intention to use a specific ephemeris.

Other than the comment above about the ephemeris argument I think this all looks very good!

If you want to change the Kitt Peak tests to 2014-09-25T00:00, then the appropriate numbers for Horizons (as generated on 2016-05-17) are:

Jupiter:
09:09:55.55, +16:51:57.8 (1 way light travel time = 49.244937 mins)
Mercury:
13 38 58.50 -13 34 42.6 (1 way light travel time = 7.699020 mins)
Moon:
12 33 12.85 -05 17 54.4 (1 way light travel time = 0.022054 mins)

@eteq - I'm done implementing your other changes. The one I did not do was change 'approximate' to 'erfa'. My main reason is that #4963 will hopefully soon give us approximate positions for the moon, and I think ideally we just use those by default. But since #4963 is not based on erfa, I felt another name was needed...

p.s. I found the problem with FLOAT_CMP, so there are now two changes to the relevant regexps: e7d271b

Note that the docs build failed because I tried to link to solar_system_ephemeris.set, but by default methods inherited from a parent class are skipped in the documentation. I think we should change this default -- see issue #4978 -- but thought I better not do this in this PR.

@mhvk - OK, gotcha re: approximate. I wonder if there's some better name than "approximate", though... "approximate" implies the others are somehow not approximate, when really they're just less approximate. What about 'builtin'? As in "using what's built-in Astropy, so you don't need to download anything"?

Other than that, 👍 from me on this now!

Yes, builtin is much better. Will change...

Tests passed, so am merging this! Thanks @StuartLittlefair and @eteq for the great collaboration in getting this in!