Precise barycentric radial velocity corrections. by StuartLittlefair · Pull Request #6861 · astropy/astropy

StuartLittlefair · 2017-11-16T10:29:31Z

This PR builds on the improvements to EarthLocation.get_gcrs_posvel in #6699, to offer radial velocity corrections compatible with Wright & Eastman (2014) to 10 mm/s (for sources at infinite distance).

Prior to this PR, radial velocity corrections between astropy and WE14 differed by several metres/second. This is due to ignoring the gravitational redshift at the observatory and a slight bug in SkyCoord.radial_velocity_correction. This PR fixes both of these issues to give much more accurate radial velocity corrections. The plot before shows the situation before (top) and after (below). The yellow star marks the position of the Sun.

As discussed in 5752 improvements to the precision of these corrections should probably be thought of as bug fixes.

The only remaining physics to include is that caused by the proper motion and parallax of the source, but I propose only including those terms once velocity support is added to SkyCoord.

This PR fixes #6697 and also fixes #6680, subject to the precision required by that user.

astropy-bot · 2017-11-16T10:29:33Z

Hi there @StuartLittlefair 👋 - thanks for the pull request! I'm just a friendly 🤖 that checks for issues related to the changelog and making sure that this pull request is milestoned and labeled correctly. This is mainly intended for the maintainers, so if you are not a maintainer you can ignore this, and a maintainer will let you know if any action is required on your part 😃.

Everything looks good from my point of view! 👍

If there are any issues with this message, please report them here.

mhvk · 2017-11-18T19:44:36Z

-        obsgeovel = CartesianRepresentation(vel_x, vel_y, vel_z)
+
+        # CIRS velocity
+        rotation_vector = CartesianRepresentation(0/u.s, 0/u.s, OMEGA_EARTH)


Now that we have velocity support in frames at least, it should be possible to not do anything at all, except to add zero velocity to the ITRS position... PR soon...

mhvk · 2017-11-18T19:51:09Z

This is great! But looking at the code, I realized that you are just reimplementing part of the frame transfer graph, and it seems to make even more sense to just use it completely. It turns out all that is needed is add a velocity of 0 to the itrs frame, and then transform to gcrs. I did this as a new get_gcrs() function, as I think this is a more logical way to present it than get_gcrs_posvel() (so the latter now just uses it internally).

I also couldn't resist some reorganizing of the redshift function, probably for no net benefit...

See https://github.com/StuartLittlefair/astropy/pull/4

mhvk · 2017-11-18T19:52:24Z

With new methods, I don't know that we can backport this...

Use frame transformations / simplify redshift

StuartLittlefair · 2017-11-20T10:05:06Z

Thanks @mhvk! Lovely tidying in that PR.

mhvk

Mostly small things, but it would seem that this might not be quite correct if the coordinates are not in ICRS.

mhvk · 2017-11-20T12:50:42Z

 from ..utils.compat.numpycompat import NUMPY_LT_1_12
 from .angles import Longitude, Latitude
-from .representation import CartesianRepresentation
+from .matrix_utilities import matrix_transpose


Should have seen this when I made the PR to your PR, but we don't need this any more.

yep, good point

mhvk · 2017-11-20T12:55:10Z

+        else:
+            # skycoord has distances so apply parallax
+            obs_icrs_cart = get_body_barycentric('earth', obstime) + gcrs_p
+            icrs_cart = self.cartesian


Don't we have to transform ourselves to ICRS for this to hold?

damn - we do! This was a bug in the original version too. Also needs conversion in the if part of the if...else... clause.

mhvk · 2017-11-20T12:58:11Z

+            targcart = self.represent_as(UnitSphericalRepresentation).to_cartesian()
+        else:
+            # skycoord has distances so apply parallax
+            obs_icrs_cart = get_body_barycentric('earth', obstime) + gcrs_p


It seems a bit wasteful to re-calculate Earth here, when it was already done for the velocity above.
Probably best to just have pos_earth, vel_earth = get_body_barycentric('earth', obstime) early on.

mhvk · 2017-11-20T13:00:44Z

+            # skycoord has distances so apply parallax
+            obs_icrs_cart = get_body_barycentric('earth', obstime) + gcrs_p
+            icrs_cart = self.cartesian
+            targcart = (icrs_cart - obs_icrs_cart).represent_as(UnitSphericalRepresentation).to_cartesian()


This is a bit slow, since you have to do all the needless trig. I'd do:

targcart = icrs_cart - obs_icrs_cart targcart /= targcart.norm()

Much better.

shbhuk · 2017-11-21T15:31:35Z

Hi,
We are working on the same code. As discussed earlier with @StuartLittlefair , the main reason for us for developing it outside of Astropy was the leap second management. Plus factors like the Earth Orientation Parameters and Grav. Time Dilation. I see that the Time dilation bit has been added now.

This is our repo
Just need to complete some final testing before we release it.

Shubham

mhvk

One typo, but perhaps more importantly, could you update the docstring of radial_velocity_correction to note that (1) barycentric is more precise, and (2) that it should be applied multiplicatively for full precision (as you describe in the full documentation).

mhvk · 2017-11-21T17:03:56Z

+            gamma_obs = 1 / np.sqrt(1 - beta_obs.norm()**2)
+            gr = location.gravitational_redshift(obstime)
+            # barycentric redshift according to eq 28 in Wright & Eastmann (2014),
+            # neglectic Shapiro delay and effects of the star's own motion


typo neglectic -> neglecting

mhvk · 2017-11-21T17:08:04Z

Approval was on the implementation (very nice!); doc-string would still be good to do (can be a [ci skip])

mhvk · 2017-11-22T13:48:38Z

OK, this all looks good from my perspective. The main question is about the milestone: @eteq, @astrofrog: this PR introduces two new methods, which are used in the radial-velocity correction. Is this OK for 2.0.3? It does correct a bug in the initial implementation.

shbhuk · 2017-11-22T20:57:35Z

What is the status of issue discussed in #5783 for keeping Astropy updated with the latest leap seconds?
If the leap second is not updated, then the correction will have a discrepancy upto 4 cm/s for each new leap second not considered. This will be due to the offset between UTC and TDB, which will affect the ephemeris and observatory posvel. This seems to be a fix which is needed in the Time conversion part of the code (get_jd12).

Until that is fixed, we should add a warning in the documentation about this potential loss of precision.

@eteq

astrofrog · 2017-11-22T23:05:31Z

I would personally vote for 2.0.3 - @eteq?

StuartLittlefair · 2017-11-23T15:56:10Z

@shbhuk we won't let the leap seconds in astropy get out of date. At the moment this means releasing a bug fix to astropy, but that will get done when necessary

indiajoe · 2017-11-23T17:45:19Z

@shbhuk Thanks for pointing this out.
@StuartLittlefair We have an observatory level pipeline which needs bary centric velocity calculation at cm/sec. The idea of updating astropy everyday automatically is a bit scary. We do update our software stacks, but prefer to do that manually to make sure nothing in pipeline gets broken.
I do understand there is a plan in future to use the leap seconds from IERS table instead. That would be much easier to update automatically.
Until then raising a warning during execution if astropy is outdated, or mentioning this issue in the documentation would be a good idea. It would at least help in debugging what happened before a wrong result is published! :)

mhvk · 2017-11-23T18:57:03Z

Agreed we really need to do this automatically, ideally from the IERS tables (which are already downloaded when needed anyway). Indeed, even to check those we do need to download the IERS tables - we cannot know when a leap second will occur when we release a new version of astropy. Help with this much appreciated! (see #5783 and links there)

eteq

@StuartLittlefair - thanks for all this - a huge improvement, and the core algorithms all LGTM. I also agree this is ok as 2.0.3 (with my caveat below).

The main concern I have is the introduction of two new (and useful) methods in a bugfix. As you'll see inline I have an idea for how to address that. If you're cool with what I suggest there, I'm happy to do the legwork on that myself so as not to delay this any further.

eteq · 2017-12-06T15:44:45Z

                                     for the location of this object at the
                                     default ``obstime``.""")

+    def get_gcrs(self, obstime):


This is nominally a new feature and should be documented as such, but then we run into the awkwardness of "new feature" despite this being a bugfix. I wonder if the way to go is to have it be private (e.g. _get_gcrs) for this PR but then add a follow-on (trivial) PR for 3.0 that just de-privatizes it and has a changelog entry mentioning it?

eteq · 2017-12-06T15:46:21Z

+        obsgeovel = gcrs_data.differentials['s'].to_cartesian()
+        return obsgeopos, obsgeovel
+
+    def gravitational_redshift(self, obstime):


similar here with the above re: privitizing this for this PR. One could definitely see people wanting to use this on its own, but we should then have it in the changelog.

eteq · 2017-12-06T15:48:11Z

+        """
+        # needs to be here to avoid circular imports
+        from .solar_system import get_body_barycentric
+        names = ('sun', 'jupiter', 'moon', 'earth')


Perhaps this should be a keyword argument with this as the default instead of hard-coded? That way if say, you thought you really cared about Saturn for whatever reason, you could include it.

(If we go with my private->public scheme this could be done in the follow-on PR. )

I think this is a great idea but will leave it to a follow-on PR to avoid holding this up, since one wants to be careful about documenting which names are OK, and matching masses with names...

eteq · 2017-12-06T15:54:08Z

+
+           sc = SkyCoord(1*u.deg, 2*u.deg)
+           vcorr = sc.rv_correction(kind='barycentric', obstime=t, location=loc)
+           rv = rv + vcorr + rv * vcorr / consts.c


Now that it's actually useful, we should add something to do this automatically. That's definitely a follow-on PR since it would be a new function/method, but now it's definitely justified.

agreed! For 3.0 I'd like to do a PR that includes the effects of the proper motion and distance if they're present on the SkyCoord so I could add a function/method to that PR

eteq · 2017-12-06T15:56:44Z

@shbhuk @indiajoe - can you create a new issue about the idea of raising a warning re: leap seconds? I personally don't agree that there should be a warning, but it's open for discussion, and I would rather not block this PR while we hash that out. Bear in mind that this problem already existed, just was masked by the missing pieces that this PR fixes, so it's not a "new" problem that needs addressing now. And more generally we definitely do not want warnings to suddenly appear due to a bugfix release.

(Or alternatively, put something about this into #5783 as stronger motivation for making #5783 happen... As detailed there we definitely want that to happen but it's non-trivial work)

eteq · 2017-12-06T16:44:00Z

Great, looks good now once the tests pass. I'll make the PR to do the publicizing (already have it ready at https://github.com/eteq/astropy/tree/publicize_vcorrs_bits but will wait on this so that no extra rebase is needed).

shbhuk · 2017-12-06T19:32:58Z

@eteq Shall post it on the #5783 thread soon.

Precise barycentric radial velocity corrections.

StuartLittlefair added 15 commits November 16, 2017 10:09

correction of get_gcrs_posvel to use CIRS velocity first

3d1cfb5

add regression test

4269fc3

changelog entry

62acd18

correct issue number

7570900

typo in changelog

5d7e98e

changelog moved to 2.0.3

207c245

use PR number, not issue number in changelog

56e7f27

remove unused import of matmul

b5eabb1

add gravitational redshift calculation to EarthLocation

851e6c8

improve precision of radial_velocity_correction

be93808

fix GR calculation for non-scalar times

18c9058

only do accurate correction calcs for barycentric corrs

1c370c7

give the comparison to barycorr teeth!

6b8f913

add docs on precision

1c187e9

minor doc tweaks

b58a254

StuartLittlefair added the coordinates label Nov 16, 2017

StuartLittlefair requested review from eteq and mhvk November 16, 2017 10:34

StuartLittlefair added 2 commits November 16, 2017 10:35

changelog entry

0a633d6

update documentation values

dc60cb6

StuartLittlefair added this to the v2.0.3 milestone Nov 16, 2017

mhvk reviewed Nov 18, 2017

View reviewed changes

mhvk added 2 commits November 18, 2017 14:47

Use frame tranformations in EarthLocation.get_gcrs_posvel.

723a1dd

Simplify gravitational_redshift a little.

beba8ea

Merge pull request #4 from mhvk/pr-6861

075ad72

Use frame transformations / simplify redshift

mhvk requested changes Nov 20, 2017

View reviewed changes

improvements as a result of mhvk review

b8ea6f2

mhvk reviewed Nov 21, 2017

View reviewed changes

mhvk approved these changes Nov 21, 2017

View reviewed changes

docstring changes and typo fix

02a024e

StuartLittlefair mentioned this pull request Nov 23, 2017

Fix gcrs posvel #6699

Closed

eteq requested changes Dec 6, 2017

View reviewed changes

privatise new methods

c781b75

eteq approved these changes Dec 6, 2017

View reviewed changes

eteq merged commit 9bfcc92 into astropy:master Dec 6, 2017

eteq mentioned this pull request Dec 6, 2017

Make barycentric vcorr-related methods public #6935

Merged

bsipocz pushed a commit that referenced this pull request Dec 8, 2017

Merge pull request #6861 from StuartLittlefair/precise_vcorrs

bd8d214

Precise barycentric radial velocity corrections.

mhvk mentioned this pull request Jan 10, 2018

TST: Solar System doctest failure on v3 #7063

Closed

This was referenced Jan 11, 2018

Implement UTC_BJD converter as additional feature. shbhuk/barycorrpy#5

Closed

Consider adding a programatic way to add leap seconds liberfa/erfa#43

Closed

eteq mentioned this pull request Feb 13, 2018

Add vmeasured option for radial velocity correction #7191

Closed

mhvk mentioned this pull request Sep 20, 2019

Potential trap with the accuracy of velocity correction #9269

Closed

mhvk mentioned this pull request Nov 20, 2020

Speed up GCRS<->(TETE, CIRS) by smarter location conversions #11069

Merged

pllim mentioned this pull request Apr 1, 2026

DOC: Remove dead link for barycorr reference #19515

Merged

1 task

Conversation

StuartLittlefair commented Nov 16, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

astropy-bot Bot commented Nov 16, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mhvk commented Nov 18, 2017

Uh oh!

mhvk commented Nov 18, 2017

Uh oh!

StuartLittlefair commented Nov 20, 2017

Uh oh!

mhvk left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

shbhuk commented Nov 21, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mhvk left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mhvk commented Nov 21, 2017

Uh oh!

mhvk commented Nov 22, 2017

Uh oh!

shbhuk commented Nov 22, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

astrofrog commented Nov 22, 2017

Uh oh!

StuartLittlefair commented Nov 23, 2017

Uh oh!

indiajoe commented Nov 23, 2017

Uh oh!

mhvk commented Nov 23, 2017

Uh oh!

eteq left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

eteq commented Dec 6, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

StuartLittlefair commented Nov 16, 2017 •

edited

Loading

astropy-bot Bot commented Nov 16, 2017 •

edited

Loading

shbhuk commented Nov 21, 2017 •

edited

Loading

shbhuk commented Nov 22, 2017 •

edited

Loading

eteq left a comment •

edited

Loading

eteq commented Dec 6, 2017 •

edited

Loading