From 3e940af1a44bc8290dd1803eac4992f068a9af30 Mon Sep 17 00:00:00 2001
From: "busko@stsci.edu" <New1trilha>
Date: Wed, 24 Jul 2019 11:16:08 -0400
Subject: [PATCH 01/12] Add support for redshift in conversions to velocity.

---
 specutils/spectra/spectrum_mixin.py |  7 +++--
 specutils/tests/test_spectrum1d.py  | 42 +++++++++++++++++++++++++++++
 2 files changed, 47 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index 0b3da9253..8d1216507 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -2,8 +2,8 @@
 from copy import deepcopy
 
 import numpy as np
-from astropy.wcs import WCSSUB_SPECTRAL
-from astropy.units import Unit
+from astropy.units import Quantity
+from astropy.constants import c as C
 from astropy import units as u
 import astropy.units.equivalencies as eq
 from astropy.utils.decorators import lazyproperty
@@ -185,6 +185,9 @@ def velocity(self):
 
         new_data = self.spectral_axis.to(u.km/u.s, equivalencies=equiv)
 
+        # Redshift-correct velocities using approximate low-z correction.
+        new_data += Quantity((self._redshift * C), new_data.unit)
+
         return new_data
 
     def with_spectral_unit(self, unit, velocity_convention=None,
diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py
index 2deacb6d8..77edfde16 100644
--- a/specutils/tests/test_spectrum1d.py
+++ b/specutils/tests/test_spectrum1d.py
@@ -121,6 +121,48 @@ def test_spectral_axis_conversions():
     new_spec = spec.with_spectral_unit(u.GHz)
 
 
+def test_redshift():
+    spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy,
+                      spectral_axis=np.array([4000, 8000]) * u.AA,
+                      velocity_convention='optical',
+                      rest_value=6000 * u.AA)
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(-99930.8, 0.5)
+    assert spec.velocity[1].value == pytest.approx(99930.8, 0.5)
+
+    spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy,
+                      spectral_axis=np.array([4000, 8000]) * u.AA,
+                      velocity_convention='optical',
+                      rest_value=6000 * u.AA,
+                      redshift = 0.1)
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(-69930.8, 0.5)
+    assert spec.velocity[1].value == pytest.approx(129930.8, 0.5)
+
+    #-------------------------
+
+    spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy,
+                      spectral_axis=np.array([10.5, 11.5]) * u.GHz,
+                      velocity_convention='radio',
+                      rest_value=11. * u.GHz)
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(13626., 1.)
+    assert spec.velocity[1].value == pytest.approx(-13626., 1.)
+
+    spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy,
+                      spectral_axis=np.array([10.5, 11.5]) * u.GHz,
+                      velocity_convention='radio',
+                      rest_value=11. * u.GHz,
+                      redshift = 0.1)
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(43606., 1.)
+    assert spec.velocity[1].value == pytest.approx(16352., 1.)
+
+
 def test_flux_unit_conversion():
     # By default the flux units should be set to Jy
     s = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy,

From 0fadb66a921fd1661dcade8ca223a77a12b67601 Mon Sep 17 00:00:00 2001
From: "busko@stsci.edu" <New1trilha>
Date: Wed, 24 Jul 2019 11:17:02 -0400
Subject: [PATCH 02/12] Add support for redshift in conversions to velocity.

---
 setup.cfg                       |  2 +-
 specutils/spectra/spectrum1d.py | 10 +++++++++-
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/setup.cfg b/setup.cfg
index dc5c2bffe..02765bf05 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -21,7 +21,7 @@ remote_data_strict = True
 asdf_schema_root = specutils/io/asdf/schemas
 # The remote data tests will run by default. Passing --remote-data=none on the
 # command line will override this setting.
-addopts = --remote-data=any --doctest-rst
+#addopts = --remote-data=any --doctest-rst
 
 [ah_bootstrap]
 auto_use = True
diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index 3832f9684..291272144 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -33,6 +33,9 @@ class Spectrum1D(OneDSpectrumMixin, NDDataRef):
         Any quantity supported by the standard spectral equivalencies
         (wavelength, energy, frequency, wave number). Describes the rest value
         of the spectral axis for use with velocity conversions.
+    redshift : float
+        The redshift correction to be applied to the spectral axis for use
+        with velocity conversions.
     uncertainty : `~astropy.nddata.NDUncertainty`
         Contains uncertainty information along with propagation rules for
         spectrum arithmetic. Can take a unit, but if none is given, will use
@@ -42,7 +45,8 @@ class Spectrum1D(OneDSpectrumMixin, NDDataRef):
         around with the spectrum container object.
     """
     def __init__(self, flux=None, spectral_axis=None, wcs=None,
-                 velocity_convention=None, rest_value=None, **kwargs):
+                 velocity_convention=None, rest_value=None, redshift=None,
+                 **kwargs):
         # Check for pre-defined entries in the kwargs dictionary.
         unknown_kwargs = set(kwargs).difference(
             {'data', 'unit', 'uncertainty', 'meta', 'mask', 'copy'})
@@ -123,6 +127,10 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None,
             elif not self._rest_value.unit.is_equivalent(u.AA) and not self._rest_value.unit.is_equivalent(u.Hz):
                 raise u.UnitsError("Rest value must be energy/wavelength/frequency equivalent.")
 
+        self._redshift = redshift
+        if redshift is None:
+            self._redshift = 0.
+
         super(Spectrum1D, self).__init__(
             data=flux.value if isinstance(flux, u.Quantity) else flux,
             wcs=wcs, **kwargs)

From 1bd075cd3dc4a15616af2850f7c3394eae6bc902 Mon Sep 17 00:00:00 2001
From: "busko@stsci.edu" <New1trilha>
Date: Wed, 24 Jul 2019 11:19:09 -0400
Subject: [PATCH 03/12] Revert change in config file.

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index 02765bf05..dc5c2bffe 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -21,7 +21,7 @@ remote_data_strict = True
 asdf_schema_root = specutils/io/asdf/schemas
 # The remote data tests will run by default. Passing --remote-data=none on the
 # command line will override this setting.
-#addopts = --remote-data=any --doctest-rst
+addopts = --remote-data=any --doctest-rst
 
 [ah_bootstrap]
 auto_use = True

From a049a29aa1e037354724cbae28ac78b477ca814c Mon Sep 17 00:00:00 2001
From: "busko@stsci.edu" <New1trilha>
Date: Wed, 24 Jul 2019 11:58:08 -0400
Subject: [PATCH 04/12] Fine tune tests.

---
 specutils/tests/test_spectrum1d.py | 40 ++++++++++++++++--------------
 1 file changed, 22 insertions(+), 18 deletions(-)

diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py
index 77edfde16..42125235f 100644
--- a/specutils/tests/test_spectrum1d.py
+++ b/specutils/tests/test_spectrum1d.py
@@ -122,45 +122,49 @@ def test_spectral_axis_conversions():
 
 
 def test_redshift():
-    spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy,
-                      spectral_axis=np.array([4000, 8000]) * u.AA,
+    spec = Spectrum1D(flux=np.array([26.0, 30., 44.5]) * u.Jy,
+                      spectral_axis=np.array([4000, 6000, 8000]) * u.AA,
                       velocity_convention='optical',
                       rest_value=6000 * u.AA)
 
     assert spec.velocity.unit == u.Unit('km/s')
-    assert spec.velocity[0].value == pytest.approx(-99930.8, 0.5)
-    assert spec.velocity[1].value == pytest.approx(99930.8, 0.5)
+    assert spec.velocity[0].value == pytest.approx(-99930.8, abs=0.5)
+    assert spec.velocity[1].value == pytest.approx(0.0, abs=0.5)
+    assert spec.velocity[2].value == pytest.approx(99930.8, abs=0.5)
 
-    spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy,
-                      spectral_axis=np.array([4000, 8000]) * u.AA,
+    spec = Spectrum1D(flux=np.array([26.0, 30., 44.5]) * u.Jy,
+                      spectral_axis=np.array([4000, 6000, 8000]) * u.AA,
                       velocity_convention='optical',
                       rest_value=6000 * u.AA,
                       redshift = 0.1)
 
     assert spec.velocity.unit == u.Unit('km/s')
-    assert spec.velocity[0].value == pytest.approx(-69930.8, 0.5)
-    assert spec.velocity[1].value == pytest.approx(129930.8, 0.5)
+    assert spec.velocity[0].value == pytest.approx(-69951.3, abs=0.5)
+    assert spec.velocity[1].value == pytest.approx(29979.2, abs=0.5)
+    assert spec.velocity[2].value == pytest.approx(129910.1, abs=0.5)
 
     #-------------------------
 
-    spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy,
-                      spectral_axis=np.array([10.5, 11.5]) * u.GHz,
+    spec = Spectrum1D(flux=np.array([26.0, 30.0, 44.5]) * u.Jy,
+                      spectral_axis=np.array([10.5, 11.0, 11.5]) * u.GHz,
                       velocity_convention='radio',
-                      rest_value=11. * u.GHz)
+                      rest_value=11.0 * u.GHz)
 
     assert spec.velocity.unit == u.Unit('km/s')
-    assert spec.velocity[0].value == pytest.approx(13626., 1.)
-    assert spec.velocity[1].value == pytest.approx(-13626., 1.)
+    assert spec.velocity[0].value == pytest.approx(13626., abs=1.)
+    assert spec.velocity[1].value == pytest.approx(0., abs=1.)
+    assert spec.velocity[2].value == pytest.approx(-13626., abs=1.)
 
-    spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy,
-                      spectral_axis=np.array([10.5, 11.5]) * u.GHz,
+    spec = Spectrum1D(flux=np.array([26.0, 30.0, 44.5]) * u.Jy,
+                      spectral_axis=np.array([10.5, 11.0, 11.5]) * u.GHz,
                       velocity_convention='radio',
-                      rest_value=11. * u.GHz,
+                      rest_value=11.0 * u.GHz,
                       redshift = 0.1)
 
     assert spec.velocity.unit == u.Unit('km/s')
-    assert spec.velocity[0].value == pytest.approx(43606., 1.)
-    assert spec.velocity[1].value == pytest.approx(16352., 1.)
+    assert spec.velocity[0].value == pytest.approx(43606., abs=1.)
+    assert spec.velocity[1].value == pytest.approx(29979., abs=1.)
+    assert spec.velocity[2].value == pytest.approx(16352., abs=1.)
 
 
 def test_flux_unit_conversion():

From 09f66747f119cf5248c726dcefbedf66f85e89f3 Mon Sep 17 00:00:00 2001
From: busko <busko@stsci.edu>
Date: Wed, 7 Aug 2019 11:36:45 -0400
Subject: [PATCH 05/12] Add support for radial velocity.

---
 specutils/spectra/spectrum1d.py     | 16 +++++++++-------
 specutils/spectra/spectrum_mixin.py |  5 ++++-
 specutils/tests/test_spectrum1d.py  | 28 ++++++++++++++++++++++++++--
 3 files changed, 39 insertions(+), 10 deletions(-)

diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index 291272144..69fa9b00f 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -33,9 +33,11 @@ class Spectrum1D(OneDSpectrumMixin, NDDataRef):
         Any quantity supported by the standard spectral equivalencies
         (wavelength, energy, frequency, wave number). Describes the rest value
         of the spectral axis for use with velocity conversions.
-    redshift : float
-        The redshift correction to be applied to the spectral axis for use
-        with velocity conversions.
+    redshift : float or `~astropy.units.Quantity`
+        The redshift/RV correction to be applied to the spectral axis for use
+        with velocity conversions. If float, the value is interpreted as the
+        redshift z. If an instance of `~astropy.units.Quantity` is passed,
+        it is interpreted as a radial velocity with the appropriate velocity units.
     uncertainty : `~astropy.nddata.NDUncertainty`
         Contains uncertainty information along with propagation rules for
         spectrum arithmetic. Can take a unit, but if none is given, will use
@@ -45,7 +47,7 @@ class Spectrum1D(OneDSpectrumMixin, NDDataRef):
         around with the spectrum container object.
     """
     def __init__(self, flux=None, spectral_axis=None, wcs=None,
-                 velocity_convention=None, rest_value=None, redshift=None,
+                 velocity_convention=None, rest_value=None, redshift_rv=None,
                  **kwargs):
         # Check for pre-defined entries in the kwargs dictionary.
         unknown_kwargs = set(kwargs).difference(
@@ -127,9 +129,9 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None,
             elif not self._rest_value.unit.is_equivalent(u.AA) and not self._rest_value.unit.is_equivalent(u.Hz):
                 raise u.UnitsError("Rest value must be energy/wavelength/frequency equivalent.")
 
-        self._redshift = redshift
-        if redshift is None:
-            self._redshift = 0.
+        self._redshift_rv = redshift_rv
+        if redshift_rv is None:
+            self._redshift_rv = 0.
 
         super(Spectrum1D, self).__init__(
             data=flux.value if isinstance(flux, u.Quantity) else flux,
diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index 8d1216507..ea46fc020 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -186,7 +186,10 @@ def velocity(self):
         new_data = self.spectral_axis.to(u.km/u.s, equivalencies=equiv)
 
         # Redshift-correct velocities using approximate low-z correction.
-        new_data += Quantity((self._redshift * C), new_data.unit)
+        if isinstance(self._redshift_rv, Quantity):
+            new_data += self._redshift_rv
+        else:
+            new_data += Quantity((self._redshift_rv * C), new_data.unit)
 
         return new_data
 
diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py
index 42125235f..78ac731df 100644
--- a/specutils/tests/test_spectrum1d.py
+++ b/specutils/tests/test_spectrum1d.py
@@ -3,6 +3,7 @@
 import gwcs
 import numpy as np
 import pytest
+from astropy.units import Quantity
 from astropy.nddata import StdDevUncertainty
 
 from .conftest import remote_access
@@ -136,7 +137,7 @@ def test_redshift():
                       spectral_axis=np.array([4000, 6000, 8000]) * u.AA,
                       velocity_convention='optical',
                       rest_value=6000 * u.AA,
-                      redshift = 0.1)
+                      redshift_rv= 0.1)
 
     assert spec.velocity.unit == u.Unit('km/s')
     assert spec.velocity[0].value == pytest.approx(-69951.3, abs=0.5)
@@ -159,13 +160,36 @@ def test_redshift():
                       spectral_axis=np.array([10.5, 11.0, 11.5]) * u.GHz,
                       velocity_convention='radio',
                       rest_value=11.0 * u.GHz,
-                      redshift = 0.1)
+                      redshift_rv= 0.1)
 
     assert spec.velocity.unit == u.Unit('km/s')
     assert spec.velocity[0].value == pytest.approx(43606., abs=1.)
     assert spec.velocity[1].value == pytest.approx(29979., abs=1.)
     assert spec.velocity[2].value == pytest.approx(16352., abs=1.)
 
+    #------------------------- radial velocity mode
+
+    spec = Spectrum1D(flux=np.array([26.0, 30., 44.5]) * u.Jy,
+                      spectral_axis=np.array([4000, 6000, 8000]) * u.AA,
+                      velocity_convention='optical',
+                      rest_value=6000 * u.AA)
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(-99930.8, abs=0.5)
+    assert spec.velocity[1].value == pytest.approx(0.0, abs=0.5)
+    assert spec.velocity[2].value == pytest.approx(99930.8, abs=0.5)
+
+    spec = Spectrum1D(flux=np.array([26.0, 30., 44.5]) * u.Jy,
+                      spectral_axis=np.array([4000, 6000, 8000]) * u.AA,
+                      velocity_convention='optical',
+                      rest_value=6000 * u.AA,
+                      redshift_rv=Quantity(1000., 'km/s'))
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(-98930.8, abs=0.5)
+    assert spec.velocity[1].value == pytest.approx(1000.0, abs=0.5)
+    assert spec.velocity[2].value == pytest.approx(100930.8, abs=0.5)
+
 
 def test_flux_unit_conversion():
     # By default the flux units should be set to Jy

From abd12572ae7e091d0268c949b94af76a22f682d1 Mon Sep 17 00:00:00 2001
From: Ivo Busko <busko@stsci.edu>
Date: Fri, 9 Aug 2019 11:03:25 -0400
Subject: [PATCH 06/12] Fix comment and doc string

---
 specutils/spectra/spectrum1d.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index 69fa9b00f..86404851b 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -33,7 +33,7 @@ class Spectrum1D(OneDSpectrumMixin, NDDataRef):
         Any quantity supported by the standard spectral equivalencies
         (wavelength, energy, frequency, wave number). Describes the rest value
         of the spectral axis for use with velocity conversions.
-    redshift : float or `~astropy.units.Quantity`
+    redshift_rv : float or `~astropy.units.Quantity`
         The redshift/RV correction to be applied to the spectral axis for use
         with velocity conversions. If float, the value is interpreted as the
         redshift z. If an instance of `~astropy.units.Quantity` is passed,
@@ -129,6 +129,8 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None,
             elif not self._rest_value.unit.is_equivalent(u.AA) and not self._rest_value.unit.is_equivalent(u.Hz):
                 raise u.UnitsError("Rest value must be energy/wavelength/frequency equivalent.")
 
+        # redshift may be moved to somewhere else in the future,
+        # perhaps to the object that handles coordinates. 
         self._redshift_rv = redshift_rv
         if redshift_rv is None:
             self._redshift_rv = 0.

From b6736adb574314130afc04025510e9ce3ba3465c Mon Sep 17 00:00:00 2001
From: Ivo Busko <busko@stsci.edu>
Date: Fri, 9 Aug 2019 11:38:05 -0400
Subject: [PATCH 07/12] Attempt to fix whitespace

---
 specutils/spectra/spectrum1d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index 86404851b..12ba46bc5 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -134,7 +134,7 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None,
         self._redshift_rv = redshift_rv
         if redshift_rv is None:
             self._redshift_rv = 0.
-
+            
         super(Spectrum1D, self).__init__(
             data=flux.value if isinstance(flux, u.Quantity) else flux,
             wcs=wcs, **kwargs)

From 780333825b8d9a3b9ec829066350791dedc9415f Mon Sep 17 00:00:00 2001
From: Ivo Busko <busko@stsci.edu>
Date: Fri, 9 Aug 2019 11:57:20 -0400
Subject: [PATCH 08/12] Fixing whitespace.

---
 specutils/spectra/spectrum1d.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index 12ba46bc5..2dd40d3b8 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -130,11 +130,11 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None,
                 raise u.UnitsError("Rest value must be energy/wavelength/frequency equivalent.")
 
         # redshift may be moved to somewhere else in the future,
-        # perhaps to the object that handles coordinates. 
+        # perhaps to the object that handles coordinates.
         self._redshift_rv = redshift_rv
         if redshift_rv is None:
             self._redshift_rv = 0.
-            
+
         super(Spectrum1D, self).__init__(
             data=flux.value if isinstance(flux, u.Quantity) else flux,
             wcs=wcs, **kwargs)

From bcc2597fdf90374c8baa1cde1268a208ecadc3a9 Mon Sep 17 00:00:00 2001
From: Erik Tollerud <erik.tollerud@gmail.com>
Date: Fri, 13 Sep 2019 13:40:11 -0400
Subject: [PATCH 09/12] switch to redshift and rv properties

---
 specutils/spectra/spectrum1d.py     | 75 ++++++++++++++++++++++++-----
 specutils/spectra/spectrum_mixin.py |  8 ---
 2 files changed, 62 insertions(+), 21 deletions(-)

diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index 2dd40d3b8..e046281cd 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -3,7 +3,8 @@
 
 import numpy as np
 from astropy import units as u
-from astropy.nddata import NDDataRef, NDUncertainty
+from astropy import constants as cnst
+from astropy.nddata import NDDataRef
 from astropy.utils.decorators import lazyproperty
 
 from ..wcs import WCSAdapter, WCSWrapper
@@ -33,11 +34,10 @@ class Spectrum1D(OneDSpectrumMixin, NDDataRef):
         Any quantity supported by the standard spectral equivalencies
         (wavelength, energy, frequency, wave number). Describes the rest value
         of the spectral axis for use with velocity conversions.
-    redshift_rv : float or `~astropy.units.Quantity`
-        The redshift/RV correction to be applied to the spectral axis for use
-        with velocity conversions. If float, the value is interpreted as the
-        redshift z. If an instance of `~astropy.units.Quantity` is passed,
-        it is interpreted as a radial velocity with the appropriate velocity units.
+    redshift
+        See `redshift` for more information.
+    radial_velocity
+        See `radial_velocity` for more information.
     uncertainty : `~astropy.nddata.NDUncertainty`
         Contains uncertainty information along with propagation rules for
         spectrum arithmetic. Can take a unit, but if none is given, will use
@@ -47,8 +47,8 @@ class Spectrum1D(OneDSpectrumMixin, NDDataRef):
         around with the spectrum container object.
     """
     def __init__(self, flux=None, spectral_axis=None, wcs=None,
-                 velocity_convention=None, rest_value=None, redshift_rv=None,
-                 **kwargs):
+                 velocity_convention=None, rest_value=None, redshift=None,
+                 radial_velocity=None, **kwargs):
         # Check for pre-defined entries in the kwargs dictionary.
         unknown_kwargs = set(kwargs).difference(
             {'data', 'unit', 'uncertainty', 'meta', 'mask', 'copy'})
@@ -129,11 +129,13 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None,
             elif not self._rest_value.unit.is_equivalent(u.AA) and not self._rest_value.unit.is_equivalent(u.Hz):
                 raise u.UnitsError("Rest value must be energy/wavelength/frequency equivalent.")
 
-        # redshift may be moved to somewhere else in the future,
-        # perhaps to the object that handles coordinates.
-        self._redshift_rv = redshift_rv
-        if redshift_rv is None:
-            self._redshift_rv = 0.
+        if redshift is None:
+            self.radial_velocity = radial_velocity
+        elif radial_velocity is None:
+            self.redshift = redshift
+        else:
+            raise ValueError('cannot set both radial_velocity and redshift at '
+                             'the same time.')
 
         super(Spectrum1D, self).__init__(
             data=flux.value if isinstance(flux, u.Quantity) else flux,
@@ -227,6 +229,53 @@ def bin_edges(self):
     def shape(self):
         return self.flux.shape
 
+    @property
+    def redshift(self):
+        """
+        The redshift(s) of the objects represented by this spectrum.  May be
+        scalar (if this spectrum's `flux` is 1D) or vector.  Note that
+        the concept of "redshift of a spectrum" can be ambiguous, so the
+        interpretation is set to some extent by either the user, or operations
+        (like template fitting) that set this attribute when they are run on
+        a spectrum.
+        """
+        return self._radial_velocity/cnst.c
+    @redshift.setter
+    def redshift(self, val):
+        if val is None:
+            self._radial_velocity = None
+        else:
+            self._radial_velocity = val * cnst.c
+
+    @property
+    def radial_velocity(self):
+        """
+        The radial velocity(s) of the objects represented by this spectrum.  May
+        be scalar (if this spectrum's `flux` is 1D) or vector.  Note that
+        the concept of "RV of a spectrum" can be ambiguous, so the
+        interpretation is set to some extent by either the user, or operations
+        (like template fitting) that set this attribute when they are run on
+        a spectrum.
+        """
+        return self._radial_velocity
+    @radial_velocity.setter
+    def radial_velocity(self, val):
+        if val is None:
+            self._radial_velocity = None
+        else:
+            if not val.unit.is_equivalent(u.km/u.s):
+                raise u.UnitsError('radial_velocity must be a velocity')
+
+            # the trick below checks if the two shapes given are broadcastable onto
+            # each other. See https://stackoverflow.com/questions/47243451/checking-if-two-arrays-are-broadcastable-in-python
+            input_shape = val.shape
+            flux_shape = self.flux.shape[:-1]
+            if all((m == n) or (m == 1) or (n == 1)
+                   for m, n in zip(input_shape[::-1], flux_shape)):
+                raise ValueError("radial_velocity or redshift must have shape that "
+                                 "is compatible with this spectrum's flux array")
+            self._radial_velocity = val
+
     @staticmethod
     def _compare_wcs(this_operand, other_operand):
         """
diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index ea46fc020..fd7290cb5 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -2,8 +2,6 @@
 from copy import deepcopy
 
 import numpy as np
-from astropy.units import Quantity
-from astropy.constants import c as C
 from astropy import units as u
 import astropy.units.equivalencies as eq
 from astropy.utils.decorators import lazyproperty
@@ -185,12 +183,6 @@ def velocity(self):
 
         new_data = self.spectral_axis.to(u.km/u.s, equivalencies=equiv)
 
-        # Redshift-correct velocities using approximate low-z correction.
-        if isinstance(self._redshift_rv, Quantity):
-            new_data += self._redshift_rv
-        else:
-            new_data += Quantity((self._redshift_rv * C), new_data.unit)
-
         return new_data
 
     def with_spectral_unit(self, unit, velocity_convention=None,

From 549e4e41dcb02f99a52917ed41cfb14feba78768 Mon Sep 17 00:00:00 2001
From: Erik Tollerud <erik.tollerud@gmail.com>
Date: Fri, 13 Sep 2019 13:47:17 -0400
Subject: [PATCH 10/12] fix incorrect docstring invokation

---
 specutils/spectra/spectrum1d.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index e046281cd..d2c3ff710 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -71,7 +71,7 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None,
             raise ValueError("Flux must be a `Quantity` object.")
 
         # In cases of slicing, new objects will be initialized with `data`
-        # instead of `flux`. Ensure we grab the `data` argument.
+        # instead of ``flux``. Ensure we grab the `data` argument.
         if flux is None and 'data' in kwargs:
             flux = kwargs.pop('data')
 
@@ -233,7 +233,7 @@ def shape(self):
     def redshift(self):
         """
         The redshift(s) of the objects represented by this spectrum.  May be
-        scalar (if this spectrum's `flux` is 1D) or vector.  Note that
+        scalar (if this spectrum's ``flux`` is 1D) or vector.  Note that
         the concept of "redshift of a spectrum" can be ambiguous, so the
         interpretation is set to some extent by either the user, or operations
         (like template fitting) that set this attribute when they are run on
@@ -251,7 +251,7 @@ def redshift(self, val):
     def radial_velocity(self):
         """
         The radial velocity(s) of the objects represented by this spectrum.  May
-        be scalar (if this spectrum's `flux` is 1D) or vector.  Note that
+        be scalar (if this spectrum's ``flux`` is 1D) or vector.  Note that
         the concept of "RV of a spectrum" can be ambiguous, so the
         interpretation is set to some extent by either the user, or operations
         (like template fitting) that set this attribute when they are run on

From 994118ef01a7ce3233ec86fcff64eda7258bb3f4 Mon Sep 17 00:00:00 2001
From: Erik Tollerud <erik.tollerud@gmail.com>
Date: Fri, 13 Sep 2019 14:15:59 -0400
Subject: [PATCH 11/12] update tests and redshift-setting

---
 specutils/spectra/spectrum_mixin.py |  7 +++++
 specutils/tests/test_spectrum1d.py  | 43 ++++++++++-------------------
 2 files changed, 22 insertions(+), 28 deletions(-)

diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index fd7290cb5..77868bcd7 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -165,6 +165,9 @@ def velocity(self):
             (wavelength, energy, frequency, wave number).
         type : {"doppler_relativistic", "doppler_optical", "doppler_radio"}
             The type of doppler spectral equivalency.
+        redshift or radial_velocity
+            If present, this shift is applied to the final output velocity to
+            get into the rest frame of the object.
 
         Returns
         -------
@@ -183,6 +186,10 @@ def velocity(self):
 
         new_data = self.spectral_axis.to(u.km/u.s, equivalencies=equiv)
 
+        # if redshift/rv is present, apply it:
+        if self._radial_velocity is not None:
+            new_data += self.radial_velocity
+
         return new_data
 
     def with_spectral_unit(self, unit, velocity_convention=None,
diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py
index 78ac731df..5c1e3cc58 100644
--- a/specutils/tests/test_spectrum1d.py
+++ b/specutils/tests/test_spectrum1d.py
@@ -3,7 +3,6 @@
 import gwcs
 import numpy as np
 import pytest
-from astropy.units import Quantity
 from astropy.nddata import StdDevUncertainty
 
 from .conftest import remote_access
@@ -128,21 +127,17 @@ def test_redshift():
                       velocity_convention='optical',
                       rest_value=6000 * u.AA)
 
-    assert spec.velocity.unit == u.Unit('km/s')
-    assert spec.velocity[0].value == pytest.approx(-99930.8, abs=0.5)
-    assert spec.velocity[1].value == pytest.approx(0.0, abs=0.5)
-    assert spec.velocity[2].value == pytest.approx(99930.8, abs=0.5)
+    assert u.allclose(spec.velocity, [-99930.8, 0, 99930.8]*u.km/u.s,
+                     atol=0.5*u.km/u.s)
 
     spec = Spectrum1D(flux=np.array([26.0, 30., 44.5]) * u.Jy,
                       spectral_axis=np.array([4000, 6000, 8000]) * u.AA,
                       velocity_convention='optical',
                       rest_value=6000 * u.AA,
-                      redshift_rv= 0.1)
+                      redshift= 0.1)
 
-    assert spec.velocity.unit == u.Unit('km/s')
-    assert spec.velocity[0].value == pytest.approx(-69951.3, abs=0.5)
-    assert spec.velocity[1].value == pytest.approx(29979.2, abs=0.5)
-    assert spec.velocity[2].value == pytest.approx(129910.1, abs=0.5)
+    assert u.allclose(spec.velocity, [-69951.3, 29979.2, 129910.1]*u.km/u.s,
+                     atol=0.5*u.km/u.s)
 
     #-------------------------
 
@@ -151,21 +146,17 @@ def test_redshift():
                       velocity_convention='radio',
                       rest_value=11.0 * u.GHz)
 
-    assert spec.velocity.unit == u.Unit('km/s')
-    assert spec.velocity[0].value == pytest.approx(13626., abs=1.)
-    assert spec.velocity[1].value == pytest.approx(0., abs=1.)
-    assert spec.velocity[2].value == pytest.approx(-13626., abs=1.)
+    assert u.allclose(spec.velocity, [13626., 0, -13626]*u.km/u.s,
+                     atol=1*u.km/u.s)
 
     spec = Spectrum1D(flux=np.array([26.0, 30.0, 44.5]) * u.Jy,
                       spectral_axis=np.array([10.5, 11.0, 11.5]) * u.GHz,
                       velocity_convention='radio',
                       rest_value=11.0 * u.GHz,
-                      redshift_rv= 0.1)
+                      redshift= 0.1)
 
-    assert spec.velocity.unit == u.Unit('km/s')
-    assert spec.velocity[0].value == pytest.approx(43606., abs=1.)
-    assert spec.velocity[1].value == pytest.approx(29979., abs=1.)
-    assert spec.velocity[2].value == pytest.approx(16352., abs=1.)
+    assert u.allclose(spec.velocity, [43606., 29979., 16352.]*u.km/u.s,
+                      atol=1*u.km/u.s)
 
     #------------------------- radial velocity mode
 
@@ -174,21 +165,17 @@ def test_redshift():
                       velocity_convention='optical',
                       rest_value=6000 * u.AA)
 
-    assert spec.velocity.unit == u.Unit('km/s')
-    assert spec.velocity[0].value == pytest.approx(-99930.8, abs=0.5)
-    assert spec.velocity[1].value == pytest.approx(0.0, abs=0.5)
-    assert spec.velocity[2].value == pytest.approx(99930.8, abs=0.5)
+    assert u.allclose(spec.velocity, [-99930.8, 0.0, 99930.8]*u.km/u.s,
+                      atol=0.5*u.km/u.s)
 
     spec = Spectrum1D(flux=np.array([26.0, 30., 44.5]) * u.Jy,
                       spectral_axis=np.array([4000, 6000, 8000]) * u.AA,
                       velocity_convention='optical',
                       rest_value=6000 * u.AA,
-                      redshift_rv=Quantity(1000., 'km/s'))
+                      radial_velocity=1000.*u.km/u.s)
 
-    assert spec.velocity.unit == u.Unit('km/s')
-    assert spec.velocity[0].value == pytest.approx(-98930.8, abs=0.5)
-    assert spec.velocity[1].value == pytest.approx(1000.0, abs=0.5)
-    assert spec.velocity[2].value == pytest.approx(100930.8, abs=0.5)
+    assert u.allclose(spec.velocity, [-98930.8, 1000.0, 100930.8]*u.km/u.s,
+                      atol=0.5*u.km/u.s)
 
 
 def test_flux_unit_conversion():

From 3c6f70abacde878441f499c769dd6fd8900e7c59 Mon Sep 17 00:00:00 2001
From: Erik Tollerud <erik.tollerud@gmail.com>
Date: Fri, 13 Sep 2019 14:24:32 -0400
Subject: [PATCH 12/12] fix initializer based on tests

---
 specutils/spectra/spectrum1d.py | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index d2c3ff710..ce0c3a46b 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -129,6 +129,13 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None,
             elif not self._rest_value.unit.is_equivalent(u.AA) and not self._rest_value.unit.is_equivalent(u.Hz):
                 raise u.UnitsError("Rest value must be energy/wavelength/frequency equivalent.")
 
+        super(Spectrum1D, self).__init__(
+            data=flux.value if isinstance(flux, u.Quantity) else flux,
+            wcs=wcs, **kwargs)
+
+        # set redshift after super() - necessary because the shape-checking
+        # requires that the flux be initialized
+
         if redshift is None:
             self.radial_velocity = radial_velocity
         elif radial_velocity is None:
@@ -137,10 +144,6 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None,
             raise ValueError('cannot set both radial_velocity and redshift at '
                              'the same time.')
 
-        super(Spectrum1D, self).__init__(
-            data=flux.value if isinstance(flux, u.Quantity) else flux,
-            wcs=wcs, **kwargs)
-
         if hasattr(self, 'uncertainty') and self.uncertainty is not None:
             if not flux.shape == self.uncertainty.array.shape:
                 raise ValueError('Flux axis ({}) and uncertainty ({}) shapes must be the same.'.format(
@@ -245,7 +248,7 @@ def redshift(self, val):
         if val is None:
             self._radial_velocity = None
         else:
-            self._radial_velocity = val * cnst.c
+            self.radial_velocity = val * cnst.c
 
     @property
     def radial_velocity(self):
@@ -270,7 +273,7 @@ def radial_velocity(self, val):
             # each other. See https://stackoverflow.com/questions/47243451/checking-if-two-arrays-are-broadcastable-in-python
             input_shape = val.shape
             flux_shape = self.flux.shape[:-1]
-            if all((m == n) or (m == 1) or (n == 1)
+            if not all((m == n) or (m == 1) or (n == 1)
                    for m, n in zip(input_shape[::-1], flux_shape)):
                 raise ValueError("radial_velocity or redshift must have shape that "
                                  "is compatible with this spectrum's flux array")