Merge pull request astropy#9980 from StuartLittlefair/use_motion_in_rvs

Use motion in rvs
lpsinger · Mar 11, 2020 · 034fc83 · 034fc83
2 parents fef6440 + d27b743
commit 034fc83
Show file tree

Hide file tree

Showing 5 changed files with 133 additions and 18 deletions.
diff --git a/CHANGES.rst b/CHANGES.rst
@@ -337,6 +337,9 @@ astropy.coordinates
 - Read-only longitudes can now be passed in to ``EarthLocation`` even if
   they include angles outside of the range of -180 to 180 degrees. [#9900]
 
+- ```SkyCoord.radial_velocity_correction``` no longer raises an Exception
+  when space motion information is present on the SkyCoord. [#9980]
+
 astropy.cosmology
 ^^^^^^^^^^^^^^^^^
 

diff --git a/astropy/coordinates/sky_coordinate.py b/astropy/coordinates/sky_coordinate.py
@@ -1,5 +1,6 @@
 import re
 import copy
+import warnings
 
 import numpy as np
 
@@ -11,6 +12,7 @@
 from astropy.utils.data_info import MixinInfo
 from astropy.utils import ShapedLikeNDArray
 from astropy.time import Time
+from astropy.utils.exceptions import AstropyUserWarning
 
 from .distances import Distance
 from .angles import Angle
@@ -1482,7 +1484,8 @@ def radial_velocity_correction(self, kind='barycentric', obstime=None,
         Use barycentric corrections if m/s precision is required.
 
         The algorithm here is sufficient to perform corrections at the mm/s level, but
-        care is needed in application. Strictly speaking, the barycentric correction is
+        care is needed in application. The barycentric correction returned uses the optical
+        approximation v = z * c. Strictly speaking, the barycentric correction is
         multiplicative and should be applied as::
 
           >>> from astropy.time import Time
@@ -1494,9 +1497,6 @@ def radial_velocity_correction(self, kind='barycentric', obstime=None,
           >>> vcorr = sc.radial_velocity_correction(kind='barycentric', obstime=t, location=loc)  # doctest: +REMOTE_DATA
           >>> rv = rv + vcorr + rv * vcorr / c  # doctest: +SKIP
 
-        If your target is nearby and/or has finite proper motion you may need to account
-        for terms arising from this. See Wright & Eastman (2014) for details.
-
         Also note that this method returns the correction velocity in the so-called
         *optical convention*::
 
@@ -1593,13 +1593,14 @@ def radial_velocity_correction(self, kind='barycentric', obstime=None,
         gcrs_p, gcrs_v = location.get_gcrs_posvel(obstime)
         # transforming to GCRS is not the correct thing to do here, since we don't want to
         # include aberration (or light deflection)? Instead, only apply parallax if necessary
+        icrs_cart = self.icrs.cartesian
+        icrs_cart_novel = icrs_cart.without_differentials()
         if self.data.__class__ is UnitSphericalRepresentation:
-            targcart = self.icrs.cartesian
+            targcart = icrs_cart_novel
         else:
             # skycoord has distances so apply parallax
             obs_icrs_cart = pos_earth + gcrs_p
-            icrs_cart = self.icrs.cartesian
-            targcart = icrs_cart - obs_icrs_cart
+            targcart = icrs_cart_novel - obs_icrs_cart
             targcart /= targcart.norm()
 
         if kind == 'barycentric':
@@ -1608,7 +1609,24 @@ def radial_velocity_correction(self, kind='barycentric', obstime=None,
             gr = location.gravitational_redshift(obstime)
             # barycentric redshift according to eq 28 in Wright & Eastmann (2014),
             # neglecting Shapiro delay and effects of the star's own motion
-            zb = gamma_obs * (1 + targcart.dot(beta_obs)) / (1 + gr/speed_of_light) - 1
+            zb = gamma_obs * (1 + beta_obs.dot(targcart)) / (1 + gr/speed_of_light)
+            # try and get terms corresponding to stellar motion.
+            if icrs_cart.differentials:
+                try:
+                    ro = self.icrs.cartesian
+                    beta_star = ro.differentials['s'].to_cartesian() / speed_of_light
+                    # ICRS unit vector at coordinate epoch
+                    ro = ro.without_differentials()
+                    ro /= ro.norm()
+                    zb *= (1 + beta_star.dot(ro)) / (1 + beta_star.dot(targcart))
+                except u.UnitConversionError:
+                    warnings.warn("SkyCoord contains some velocity information, but not enough to "
+                                  "calculate the full space motion of the source, and so this has "
+                                  "been ignored for the purposes of calculating the radial velocity "
+                                  "correction. This can lead to errors on the order of metres/second.",
+                                  AstropyUserWarning)
+
+            zb = zb - 1
             return zb * speed_of_light
         else:
             # do a simpler correction ignoring time dilation and gravitational redshift

diff --git a/astropy/coordinates/tests/test_velocity_corrs.py b/astropy/coordinates/tests/test_velocity_corrs.py
@@ -301,3 +301,95 @@ def test_invalid_argument_combos():
 
     with pytest.raises(ValueError):
         scwattrs.radial_velocity_correction(timel)
+
+
+@pytest.mark.remote_data
+def test_regression_9645():
+    sc = SkyCoord(10*u.deg, 20*u.deg, distance=5*u.pc,
+                  pm_ra_cosdec=0*u.mas/u.yr, pm_dec=0*u.mas/u.yr, radial_velocity=0*u.km/u.s)
+    sc_novel = SkyCoord(10*u.deg, 20*u.deg, distance=5*u.pc)
+    corr = sc.radial_velocity_correction(obstime=test_input_time, location=test_input_loc)
+    corr_novel = sc_novel.radial_velocity_correction(obstime=test_input_time, location=test_input_loc)
+    assert_quantity_allclose(corr, corr_novel)
+
+
+@pytest.mark.remote_data
+def test_barycorr_withvels():
+    # this is the result of calling _get_barycorr_bvcs_withvels
+    barycorr_bvcs = u.Quantity(
+        [-10335.94901398, -14198.49074135, -2237.58603184,
+         -14198.49030979, -17425.47907834, -17131.72454389,
+         2424.38453928, 2130.62856716, -17425.47852064,
+         -19872.51346598, -24442.3888212, -11017.09452282,
+         6978.07440406, 11547.94841586, -1877.34538839,
+         -19872.51256396, -21430.09440631, -27669.15820551,
+         -16917.09266796, 2729.57728906, 16476.50782874,
+         13971.98025823, -2898.04441526, -21430.09317332,
+         -22028.52360593, -29301.93584671, -21481.13829465,
+         -3147.44812726, 14959.5072031, 22232.91879544,
+         14412.12149724, -3921.56784245, -22028.52196585,
+         -21641.02245985, -29373.06002768, -24205.91188603,
+         -8557.34397026, 10250.50467491, 23417.23250086,
+         24781.98159696, 13706.17139042, -4627.70476083,
+         -21641.02032556, -20284.9303734, -28193.92149786,
+         -22908.5203847, -6901.82469681, 12336.45463077,
+         25804.51284952, 27200.49627844, 15871.20965329,
+         -2882.25075298, -20284.92789334, -18020.92912708,
+         -25752.96566033, -20585.82158138, -4937.26072294,
+         13870.58130079, 27037.30623272, 28402.05764536,
+         17326.25428387, -1007.62298192, -18020.92629577,
+         -14950.3271854, -22223.73822407, -14402.95127754,
+         3930.72291792, 22037.66667233, 29311.07829935,
+         21490.29301495, 3156.62395395, -14950.32432903,
+         -11210.52709262, -17449.59090902, -6697.54579712,
+         12949.09890328, 26696.01926168, 24191.50484946,
+         7321.507315, -11210.5243125, -6968.87637548,
+         -11538.75489728, 1886.5050118, 19881.64323346,
+         24451.52233565, 11026.26444463, -6968.87377902,
+         -2415.17913524, -2121.44617927, 17434.60437916,
+         17140.87243824, -2415.17725544, 2246.79681159,
+         14207.61324454, 2246.79782408, 6808.43882788], u.m/u.s)
+    coos = _get_test_input_radecvels()
+    bvcs_astropy = coos.radial_velocity_correction(obstime=test_input_time, location=test_input_loc)
+    assert_quantity_allclose(bvcs_astropy, barycorr_bvcs, atol=10*u.mm/u.s)
+    return bvcs_astropy, barycorr_bvcs  # for interactively examination
+
+
+def _get_test_input_radecvels():
+    coos = _get_test_input_radecs()
+    ras = coos.ra
+    decs = coos.dec
+    pmra = np.linspace(-1000, 1000, coos.size)*u.mas/u.yr
+    pmdec = np.linspace(0, 1000, coos.size)*u.mas/u.yr
+    rvs = np.linspace(0, 100, coos.size)*u.km/u.s
+    distance = np.linspace(10, 1000, coos.size)*u.pc
+    return SkyCoord(ras, decs, pm_ra_cosdec=pmra, pm_dec=pmdec, radial_velocity=rvs, distance=distance)
+
+
+def _get_barycorr_bvcs_withvels(coos, loc, injupyter=False):
+    """
+    Gets the barycentric correction of the test data from the
+    http://astroutils.astronomy.ohio-state.edu/exofast/barycorr.html web site.
+    Requires the https://github.com/tronsgaard/barycorr python interface to that
+    site.
+
+    Provided to reproduce the test data above, but not required to actually run
+    the tests.
+    """
+    import barycorr
+    from astropy.utils.console import ProgressBar
+
+    bvcs = []
+    for coo in ProgressBar(coos, ipython_widget=injupyter):
+        res = barycorr.bvc(test_input_time.utc.jd,
+                           coo.ra.deg, coo.dec.deg,
+                           lat=loc.geodetic[1].deg,
+                           lon=loc.geodetic[0].deg,
+                           pmra=coo.pm_ra_cosdec.to_value(u.mas/u.yr),
+                           pmdec=coo.pm_dec.to_value(u.mas/u.yr),
+                           parallax=coo.distance.to_value(u.mas, equivalencies=u.parallax()),
+                           rv=coo.radial_velocity.to_value(u.m/u.s),
+                           epoch=test_input_time.utc.jd,
+                           elevation=loc.geodetic[2].to(u.m).value)
+        bvcs.append(res)
+    return bvcs*u.m/u.s
diff --git a/docs/coordinates/velocities.rst b/docs/coordinates/velocities.rst
@@ -518,9 +518,11 @@ Precision of `~astropy.coordinates.SkyCoord.radial_velocity_correction`
 ------------------------------------------------------------------------
 
 The correction computed by `~astropy.coordinates.SkyCoord.radial_velocity_correction`
-can be added to any observed radial velocity to provide a correction that is
-accurate to a level of approximately 3 m/s. If you need more precise
-corrections, there are a number of subtleties of which you must be aware.
+uses the optical approximation :math:`v = zc` (see :ref:`astropy-units-dopper-equivalencies`
+for details). The corretion can be added to any observed radial velocity
+to provide a correction that is accurate to a level of approximately 3 m/s.
+If you need more precise corrections, there are a number of subtleties of
+which you must be aware.
 
 The first is that you should always use a barycentric correction, as the
 barycenter is a fixed point where gravity is constant. Since the heliocenter
@@ -534,7 +536,7 @@ surface. For these reasons, the barycentric correction in
 be used for high precision work.
 
 Other considerations necessary for radial velocity corrections at the cm/s
-level are outlined in `Wright & Eastmann (2014) <http://adsabs.harvard.edu/abs/2014PASP..126..838W>`_.
+level are outlined in `Wright & Eastman (2014) <http://adsabs.harvard.edu/abs/2014PASP..126..838W>`_.
 Most important is that the barycentric correction is, strictly speaking,
 *multiplicative*, so that you should apply it as:
 
@@ -550,9 +552,7 @@ the barycentric correction in this way leads to errors of order 3 m/s.
 The barycentric correction in `~astropy.coordinates.SkyCoord.radial_velocity_correction` is consistent
 with the `IDL implementation <http://astroutils.astronomy.ohio-state.edu/exofast/barycorr.html>`_ of
 the Wright & Eastmann (2014) paper to a level of 10 mm/s for a source at
-infinite distance. We do not include the Shapiro delay, nor any effect related
-to the finite distance or proper motion of the source. The Shapiro delay is
-unlikely to be important unless you seek mm/s precision, but the effects of the
-source's parallax and proper motion can be important at the cm/s level. These
-effects are likely to be added to future versions of Astropy along with
-velocity support for |skycoord|, but in the meantime see `Wright & Eastmann (2014) <http://adsabs.harvard.edu/abs/2014PASP..126..838W>`_.
+infinite distance. We do not include the Shapiro delay nor the light
+travel time correction from equation 28 of that paper. The neglected terms
+are not important unless you require accuracies of better than 1 cm/s.
+If you do require that precision, see `Wright & Eastmann (2014) <http://adsabs.harvard.edu/abs/2014PASP..126..838W>`_.
diff --git a/docs/units/equivalencies.rst b/docs/units/equivalencies.rst
@@ -143,6 +143,8 @@ These equivalencies even work with non-base units::
   >>> imperial.inch.to(imperial.Cal, equivalencies=u.spectral())  # doctest: +FLOAT_CMP
   1.869180759162485e-27
 
+.. _astropy-units-dopper-equivalencies:
+
 Spectral (Doppler) equivalencies
 --------------------------------