diff --git a/src/exoplanet/orbits/keplerian.py b/src/exoplanet/orbits/keplerian.py
index 1f3ef98b..25626a13 100644
--- a/src/exoplanet/orbits/keplerian.py
+++ b/src/exoplanet/orbits/keplerian.py
@@ -567,6 +567,30 @@ def get_relative_angles(self, t, parallax=None, light_delay=False):
 
         return (rho, theta)
 
+    def get_star_relative_angles(self, t, parallax=None, light_delay=False):
+        """The stars' relative position to the star in the sky plane, in
+        separation, position angle coordinates.
+        .. note:: This treats each planet independently and does not take the
+            other planets into account when computing the position of the
+            star. This is fine as long as the planet masses are small.
+        Args:
+            t: The times where the position should be evaluated.
+            light_delay: account for the light travel time delay? Default is
+                False.
+        Returns:
+            The separation (arcseconds) and position angle (radians,
+            measured east of north) of the star relative to the barycentric frame.
+        """
+        X, Y, Z = self._get_position(
+            -self.a_star, t, parallax, light_delay=light_delay
+        )
+
+        # calculate rho and theta
+        rho = tt.squeeze(tt.sqrt(X**2 + Y**2))  # arcsec
+        theta = tt.squeeze(tt.arctan2(Y, X))  # radians between [-pi, pi]
+
+        return (rho, theta)
+
     def _get_velocity(self, m, t):
         """Get the velocity vector of a body in the observer frame"""
         sinf, cosf = self._get_true_anomaly(t)
diff --git a/tests/orbits/keplerian_test.py b/tests/orbits/keplerian_test.py
index dc2e9819..1b60b0d4 100644
--- a/tests/orbits/keplerian_test.py
+++ b/tests/orbits/keplerian_test.py
@@ -696,3 +696,74 @@ def test_jacobians():
         orbit.jacobians["b"]["cos_incl"].eval(),
         theano.grad(orbit.cos_incl, bv).eval(),
     )
+
+
+def test_relative_angles():
+
+    # test separation and position angle with Earth and Sun
+    p_earth = 365.256
+    t = np.linspace(0, 1000, 1000)
+    m_earth = 1.0 * 3.00273e-6  # units m_sun
+    orbit_earth = KeplerianOrbit(
+        m_star=1.0,
+        r_star=1.0,
+        t0=0.5,
+        period=p_earth,
+        ecc=0.0167,
+        omega=np.radians(102.9),
+        Omega=np.radians(0.0),
+        incl=np.radians(45.0),
+        m_planet=m_earth,
+    )
+
+    rho_star_earth, theta_star_earth = theano.function(
+        [], orbit_earth.get_star_relative_angles(t, parallax=0.1)
+    )()
+    rho_earth, theta_earth = theano.function(
+        [], orbit_earth.get_relative_angles(t, parallax=0.1)
+    )()
+
+    rho_star_earth_diff = np.max(rho_star_earth) - np.min(rho_star_earth)
+    rho_earth_diff = np.max(rho_earth) - np.min(rho_earth)
+
+    # make sure amplitude of separation is correct for star and planet motion
+    assert np.isclose(rho_earth_diff, 3.0813126e-02)
+    assert np.isclose(rho_star_earth_diff, 9.2523221e-08)
+
+    # make sure planet and star position angle are the same relative to each other
+    assert np.allclose(theta_earth, theta_star_earth)
+
+    ########################################
+    ########################################
+    # test separation and position angle with Jupiter and Sun
+    p_jup = 4327.631
+    t = np.linspace(0, 10000, 10000)
+    m_jup = 317.83 * 3.00273e-6  # units m_sun
+    orbit_jup = KeplerianOrbit(
+        m_star=1.0,
+        r_star=1.0,
+        t0=0.5,
+        period=p_jup,
+        ecc=0.0484,
+        omega=np.radians(274.3) - 2 * np.pi,
+        Omega=np.radians(100.4),
+        incl=np.radians(45.0),
+        m_planet=m_jup,
+    )
+
+    rho_star_jup, theta_star_jup = theano.function(
+        [], orbit_jup.get_star_relative_angles(t, parallax=0.1)
+    )()
+    rho_jup, theta_jup = theano.function(
+        [], orbit_jup.get_relative_angles(t, parallax=0.1)
+    )()
+
+    rho_star_jup_diff = np.max(rho_star_jup) - np.min(rho_star_jup)
+    rho_jup_diff = np.max(rho_jup) - np.min(rho_jup)
+
+    # make sure amplitude of separation is correct for star and planet motion
+    assert np.isclose(rho_jup_diff, 1.7190731e-01)
+    assert np.isclose(rho_star_jup_diff, 1.6390463e-04)
+
+    # make sure planet and star position angle are the same relative to each other
+    assert np.allclose(theta_jup, theta_star_jup)