Enable calculation of barycentric velocities in get_body_barycentric

[mhvk]

This is a first step towards being able to compute barycentric corrections (#3544). It allows one to calculate barycentric velocities of solar system bodies, including the Earth, either based on builtin ephemerides or JPL ones.

Probably need some discussion on how exactly to implement getting a barycentric correction (in particular, should it be a method of SkyCoord or perhaps of Time, like light_travel_time).

[mhvk]

Just to give a sense of what would be needed to calculate barycentric velocity corrections:

from astropy.time import Time
from astropy.coordinates import SkyCoord, solar_system
from astropy import units as u
t = Time('2016-03-20T12:30:00')  # near Spring equinox
ep, ev = solar_system.get_body_barycentric('earth', t, get_velocity=True)
sc = SkyCoord(['12h00m00s', '18h00m00s'], ['+0d00m00s', '-23d00m00s'])
sc.cartesian.xyz.T.dot(ev.xyz).to(u.km/u.s)
# <Quantity [  0.42148212, 29.87764459] km / s>

Here, for the coordinate, a real routine would need to ensure the skycoord is in ICRS and that it has unit length (as in Time.light_travel_time).

[mhvk]

@adrn, @eteq, @StuartLittlefair - would you be able to review this? While it is only a small step towards having barycentric velocities, it is an essential one, and I think fairly independent.

One question I wavered about was whether to given get_body_barycentric the option to return velocities as well, or whether to create a new routine that always returns two items (but with the velocity not necessarily filled) and then let get_body_barycentric use that. I in fact tried both and ended up with the first option, but can be convinced otherwise (arguments against: the number of results a function returns should not depend on arguments; in favour: otherwise get_body_barycentric is just a shell).

[StuartLittlefair]

Is there a reason for using a different array here? Can't you just use body_pv_bary throughout?

[StuartLittlefair]

@mhvk This looks good! I don't know what you want to do about the bug I just found. My instinct is just to fix it here, but I could also raise a separate issue for that.

I don't have a strong opinion on functions return different numbers of items depending on arguments. It's pretty explicit that specifying get_velocity=True means you also get a velocity so I'm happy with the way this has been done.

[mhvk]

@StuartLittlefair - indeed, multi-dimensional times broke with jplephem. I agree one might as well fix it here, so did that. Also followed your suggestion of not introducing new names unnecessarily.

[StuartLittlefair]

I still think this should be posvel = spk.compute(jd1, jd2), otherwise you're ignoring the pairs in kernel_spec

[mhvk]

Duh, that's the problem with untested code paths. But I don't really know any ephemeris with data_type=3...

[StuartLittlefair]

Hi @mhvk - Just one minor thing outstanding and then I'm perfectly happy with this.

[mhvk]

OK, I rebased to get rid of the astropy-helpers version change that somehow snuck in there. Otherwise, I think this is done.

@eteq - I think @StuartLittlefair did a thorough review, but the one question we have not quite answered is whether the present approach of adding a get_velocity argument that causes the output to change from a single position to a tuple of pos & vel is the best one (see #5231 (comment)). What do you think?

[eteq]

This "note" part should perhaps go in the "Notes" section instead?

[eteq]

@mhvk - sorry for the delay on this. I'm very much against a function that returns substantially different things based on a keyword argument. It tends to make for less readable code because ideally you can look at the name of a function and have that tell you what it's doing. So I would change this to get_body_barycentric_posvel or similar, and do just as you suggested by adding back in a get_body_barycentric which uses what's written here. The one caveat is how much of performance impact this has. If it's substantial another option is to change this function to _get_body_barycentric_all and keep the get_velocity argument, and then have get_body_barycentric_posvel and get_body_barycentric by small "shell" functions that use that as the implementation.

I reviewed the rest of it, though, and I'm happy aside from the small docstring suggestion.

[mhvk]

@eteq - OK, I made it two routines instead, which both use the function as I had it (which now is private). I agree it is cleaner, even though I'm not all that happy with a very long name like get_body_barycentric_posvel. But hopefully, once SkyCoord can handle velocities, this may become a classmethod or so.

p.s. I also made the docstring change, though it is in the private routine now.

[mhvk]

@eteq - I think I addressed your comments. OK to merge?

[eteq]

@mhvk - yep, I'm with you that get_body_barycentric_posvel is awkwardly long, but I think it's better than the alternative.

One oddity, though: it appears that only appveyor ran... weird, but I think I can get the tests to run (immediately) by just pushing it to my own fork. Will report back when it finishes (and if it passes I'll go ahead and merge)

[mhvk]

@eteq - did the travis run work on your branch? Otherwise, I think reopening and closing this PR would do the trick, or I can just do a rebase and force-push.

[eteq]

@mhvk - yes: https://travis-ci.org/eteq/astropy/builds/158024088

The coverage error is ignorable, but there's a sphinx warning. Since I'm already looking at it, I'll see if I can fix it.

[eteq]

Ah, I think I figured it out, @mhvk - see mhvk#19 . The test for that is running in https://travis-ci.org/eteq/astropy/builds/158215451 so if that passes you can merge that and then we can merge this.

[mhvk]

Cancelled the builds here, since this is already being tested at https://travis-ci.org/eteq/astropy/builds/158215451

[mhvk]

Those travis builds now passed; since appveyor was already OK, I'm merging...