diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md
index 0d9114e57..b5a26d4a0 100644
--- a/CHANGELOG-unreleased.md
+++ b/CHANGELOG-unreleased.md
@@ -14,7 +14,8 @@ the released changes.
 - Made `Residuals` independent of `GLSFitter` (GLS chi2 is now computed using the new function `Residuals._calc_gls_chi2()`).
 ### Added
 - Added `WaveX` model as a `DelayComponent` with Fourier amplitudes as fitted parameters
-- Added DMWaveX model (Fourier representation of DM noise)
+- Added radial velocity methods for binary models
+- Added `DMWaveX` model (Fourier representation of DM noise)
 ### Fixed
 - Wave model `validate()` can correctly use PEPOCH to assign WAVEEPOCH parameter
 - Fixed RTD by specifying theme explicitly.
diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py
index 1f02b8c75..6cc845412 100644
--- a/src/pint/models/timing_model.py
+++ b/src/pint/models/timing_model.py
@@ -799,6 +799,88 @@ def orbital_phase(self, barytimes, anom="mean", radians=True):
         # return with radian units or return as unitless cycles from 0-1
         return anoms * u.rad if radians else anoms / (2 * np.pi)
 
+    def pulsar_radial_velocity(self, barytimes):
+        """Return line-of-sight velocity of the pulsar relative to the system barycenter at barycentric MJD times.
+
+        Parameters
+        ----------
+        barytimes: Time, TOAs, array-like, or float
+            MJD barycentric time(s). The times to compute the
+            orbital phases.  Needs to be a barycentric time in TDB.
+            If a TOAs instance is passed, the barycentering will happen
+            automatically.  If an astropy Time object is passed, it must
+            be in scale='tdb'.  If an array-like object is passed or
+            a simple float, the time must be in MJD format.
+
+        Raises
+        ------
+        ValueError
+            If an astropy Time object is passed with scale!="tdb".
+
+        Returns
+        -------
+        array
+            The line-of-sight velocity
+
+        Notes
+        -----
+        This is the radial velocity of the pulsar.
+
+        See [1]_
+
+        .. [1] Lorimer & Kramer, 2008, "The Handbook of Pulsar Astronomy", Eqn. 8.24
+        """
+        # this should also update the binary instance
+        nu = self.orbital_phase(barytimes, anom="true")
+        b = self.components[
+            [x for x in self.components.keys() if x.startswith("Binary")][0]
+        ]
+        bbi = b.binary_instance  # shorthand
+        psi = nu + bbi.omega()
+        return (
+            2
+            * np.pi
+            * bbi.a1()
+            / (bbi.pb() * np.sqrt(1 - bbi.ecc() ** 2))
+            * (np.cos(psi) + bbi.ecc() * np.cos(bbi.omega()))
+        ).cgs
+
+    def companion_radial_velocity(self, barytimes, massratio):
+        """Return line-of-sight velocity of the companion relative to the system barycenter at barycentric MJD times.
+
+        Parameters
+        ----------
+        barytimes: Time, TOAs, array-like, or float
+            MJD barycentric time(s). The times to compute the
+            orbital phases.  Needs to be a barycentric time in TDB.
+            If a TOAs instance is passed, the barycentering will happen
+            automatically.  If an astropy Time object is passed, it must
+            be in scale='tdb'.  If an array-like object is passed or
+            a simple float, the time must be in MJD format.
+        massratio : float
+            Ratio of pulsar mass to companion mass
+
+
+        Raises
+        ------
+        ValueError
+            If an astropy Time object is passed with scale!="tdb".
+
+        Returns
+        -------
+        array
+            The line-of-sight velocity
+
+        Notes
+        -----
+        This is the radial velocity of the companion.
+
+        See [1]_
+
+        .. [1] Lorimer & Kramer, 2008, "The Handbook of Pulsar Astronomy", Eqn. 8.24
+        """
+        return -self.pulsar_radial_velocity(barytimes) * massratio
+
     def conjunction(self, baryMJD):
         """Return the time(s) of the first superior conjunction(s) after baryMJD.
 
diff --git a/tests/test_rv.py b/tests/test_rv.py
new file mode 100644
index 000000000..6ea9d4293
--- /dev/null
+++ b/tests/test_rv.py
@@ -0,0 +1,43 @@
+import pytest
+import numpy as np
+from astropy import units as u, constants as c
+
+from pint.models import get_model
+import pint.simulation
+import pint.binaryconvert
+from pinttestdata import datadir
+import os
+
+
+class TestRV:
+    def setup_method(self):
+        self.m = get_model(
+            os.path.join(datadir, "B1855+09_NANOGrav_dfg+12_modified_DD.par")
+        )
+        # set the eccentricity to nonzero
+        # but not too high so that the models with various approximations won't fail
+        self.m.ECC.value = 0.01
+        self.ts = self.m.T0.value + np.linspace(0, 2 * self.m.PB.value, 101)
+        self.ts = pint.simulation.make_fake_toas_fromMJDs(self.ts, self.m)
+
+        self.v = self.m.pulsar_radial_velocity(self.ts)
+
+    def test_rv_basemodel(self):
+        nu = self.m.orbital_phase(self.ts, anom="true", radians=True)
+        # HBOPA 8.24
+        v = (
+            (2 * np.pi / self.m.PB.quantity)
+            * (self.m.A1.quantity / np.sqrt(1 - self.m.ECC.quantity**2))
+            * (
+                np.cos(nu + self.m.OM.quantity)
+                + self.m.ECC.quantity * np.cos(self.m.OM.quantity)
+            )
+        )
+        assert np.allclose(self.v, v)
+
+    @pytest.mark.parametrize("othermodel", ["ELL1", "ELL1H", "DDS", "BT"])
+    def test_rv_othermodels(self, othermodel):
+        mc = pint.binaryconvert.convert_binary(self.m, othermodel)
+        vc = mc.pulsar_radial_velocity(self.ts)
+        # have a generous tolerance here since some of the models don't work well for high ECC
+        assert np.allclose(vc, self.v, atol=20 * u.km / u.s, rtol=1e-2)