diff --git a/docs/user_guide/calculate_response.rst b/docs/user_guide/calculate_response.rst
index 10b580a4..4df8da0d 100644
--- a/docs/user_guide/calculate_response.rst
+++ b/docs/user_guide/calculate_response.rst
@@ -52,18 +52,18 @@ This function is roughly equivalent to:
 
         Parameters
         ----------
-        rao : obj
-            ``RAO`` object.
-        wave : obj
-            ``WaveSpectrum`` object.
+        rao : RAO
+            Response amplitude operator (RAO).
+        wave : WaveSpectrum
+            Wave spectrum.
         heading : float
-            Heading of vessel relative to wave spectrum coordinate system.
+            Heading of vessel relative to wave spectrum reference frame.
         heading_degrees : bool
             Whether the heading is given in 'degrees'. If ``False``, 'radians' is assumed.
 
         Returns
         -------
-        obj :
+        DirectionalSpectrum :
             Response spectrum.
         """
 
@@ -73,10 +73,7 @@ This function is roughly equivalent to:
         # Ensure that ``rao`` and ``wave`` has the same 'wave convention'
         wave_body.set_wave_convention(**rao.wave_convention)
 
-        # Reshape ``rao`` and ``wave`` so that they share the same frequency/direction
-        # coordinates. In this example, ``wave`` will dictate the coordinates, and
-        # the ``rao`` object will be interpolated to match these coordinates.
-        # 
+        # Reshape the RAO to match the wave spectrum's frequency/direction coordinates.
         # It is recommended to reshape (i.e., interpolate) the magnitude-squared
         # version of the RAO when estimating response, since this has shown best
         # results:
@@ -85,9 +82,8 @@ This function is roughly equivalent to:
         dirs = wave_body.dirs(degrees=False)
         rao_squared = (rao * rao.conjugate()).real
         rao_squared = rao_squared.reshape(freq, dirs, freq_hz=False, degrees=False)
-        wave_body = wave_body.reshape(freq, dirs, freq_hz=False, degrees=False)
 
-        return wr.multiply(rao_squared, wave_body, output_type="directional_spectrum")
+        return wr.multiply(rao_squared, wave_body, output_type="DirectionalSpectrum")
 
 The response is returned as a :class:`~waveresponse.DirectionalSpectrum` object,
 and provides useful spectrum operations, such as:
diff --git a/src/waveresponse/_core.py b/src/waveresponse/_core.py
index 24d33894..d3df80ae 100644
--- a/src/waveresponse/_core.py
+++ b/src/waveresponse/_core.py
@@ -2304,60 +2304,122 @@ def dirm(self, degrees=None):
         return dirm
 
 
+def _calculate_response_deprecated(rao, wave_body, coord_freq, coord_dirs):
+    # TODO: Deprecated. Remove this function in a future release.
+
+    if coord_freq == "wave":
+        freq = wave_body._freq
+    elif coord_freq == "rao":
+        freq = rao._freq
+    else:
+        raise ValueError("Invalid `coord_freq` value. Should be 'wave' or 'rao'.")
+    if coord_dirs == "wave":
+        dirs = wave_body._dirs
+    elif coord_dirs == "rao":
+        dirs = rao._dirs
+    else:
+        raise ValueError("Invalid `coord_dirs` value. Should be 'wave' or 'rao'.")
+
+    rao_squared = (rao * rao.conjugate()).real
+    rao_squared = rao_squared.reshape(freq, dirs, freq_hz=False, degrees=False)
+    wave_body = wave_body.reshape(freq, dirs, freq_hz=False, degrees=False)
+
+    return multiply(rao_squared, wave_body, output_type="DirectionalSpectrum")
+
+
 def calculate_response(
-    rao, wave, heading, heading_degrees=False, coord_freq="wave", coord_dirs="wave"
+    rao,
+    wave,
+    heading,
+    heading_degrees=False,
+    reshape="rao_squared",
+    coord_freq=None,
+    coord_dirs=None,
 ):
     """
     Calculate response spectrum.
 
+    The response spectrum is calculated according to:
+
+        S_x(w, theta) = H(w, theta) * H*(w, theta) * S_w(w, theta)
+
+    where S_x(w, theta) denotes the response spectrum, H(w, theta) denotes the RAO,
+    H*(w, theta) denotes the RAO conjugate, and S_w(w, theta) denotes the wave
+    spectrum (expressed in the RAO's reference frame).
+
+    The frequency and direction coordinates are dictatated by the wave spectrum.
+    I.e., the RAO (or the magnitude-squared verison of it) is interpolated to match
+    the grid coordinates of the wave spectrum.
+
     Parameters
     ----------
-    rao : obj
-        Response amplitude operator (RAO) as a :class:`~waveresponse.RAO` object.
-    wave : obj
-        2-D wave spectrum as a :class:`~waveresponse.WaveSpectrum` object.
+    rao : RAO
+        Response amplitude operator (RAO).
+    wave : WaveSpectrum or WaveBinSpectrum
+        2-D wave spectrum.
     heading : float
         Heading of vessel relative to wave spectrum coordinate system.
     heading_degrees : bool
         Whether the heading is given in 'degrees'. If ``False``, 'radians' is assumed.
+    reshape : {'rao', 'rao_squared'}, default 'rao_squared'
+        Determines whether to reshape the RAO or the magnitude-squared version of
+        the RAO before pairing with the wave spectrum. Linear interpolation is
+        performed to match the frequency and direction coordinates of the wave
+        spectrum.
     coord_freq : str, optional
-        Frequency coordinates for interpolation. Should be 'wave' or 'rao'. Determines
-        if it is the wave spectrum or the RAO that should dictate which frequencies
-        to use in response calculation. The other object will be interpolated to
-        match these frequencies.
+        Deprecated; use `reshape` instead. Frequency coordinates for interpolation.
+        Should be 'wave' or 'rao'. Determines if it is the wave spectrum or the
+        RAO that should dictate which frequencies to use in response calculation.
+        The other object will be interpolated to match these frequencies.
     coord_dirs : str, optional
-        Direction coordinates for interpolation. Should be 'wave' or 'rao'. Determines
-        if it is the wave spectrum or the RAO that should dictate which directions
-        to use in response calculation. The other object will be interpolated to
-        match these directions.
+        Deprecated; use `reshape` instead. Direction coordinates for interpolation.
+        Should be 'wave' or 'rao'. Determines if it is the wave spectrum or the
+        RAO that should dictate which directions to use in response calculation.
+        The other object will be interpolated to match these directions.
 
     Returns
     -------
-    obj :
-        Response spectrum as :class:`DirectionalSpectrum` object.
+    DirectionalSpectrum or DirectionalBinSpectrum :
+        Response spectrum.
     """
     wave_body = wave.rotate(heading, degrees=heading_degrees)
     wave_body.set_wave_convention(**rao.wave_convention)
 
-    if coord_freq.lower() == "wave":
-        freq = wave_body._freq
-    elif coord_freq.lower() == "rao":
-        freq = rao._freq
+    # TODO: Remove once the deprecation period is over
+    if coord_freq and coord_dirs:
+        warnings.warn(
+            "The `coord_freq` and `coord_dirs` parameters are deprecated and will be removed in a future release."
+            "Use the `reshape` parameter instead.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        return _calculate_response_deprecated(rao, wave_body, coord_freq, coord_dirs)
+    elif coord_freq or coord_dirs:
+        raise ValueError("Both `coord_freq` and `coord_dirs` must be provided.")
+
+    freq, dirs = wave_body._freq, wave_body._dirs
+    if reshape == "rao":
+        rao = rao.reshape(freq, dirs, freq_hz=False, degrees=False)
+        rao_squared = (rao * rao.conjugate()).real
+    elif reshape == "rao_squared":
+        rao_squared = (rao * rao.conjugate()).real
+        rao_squared = rao_squared.reshape(freq, dirs, freq_hz=False, degrees=False)
     else:
-        raise ValueError("Invalid `coord_freq` value. Should be 'wave' or 'rao'.")
+        raise ValueError("Invalid `reshape` value. Should be 'rao' or 'rao_squared'.")
 
-    if coord_dirs.lower() == "wave":
-        dirs = wave_body._dirs
-    elif coord_dirs.lower() == "rao":
-        dirs = rao._dirs
-    else:
-        raise ValueError("Invalid `coord_dirs` value. Should be 'wave' or 'rao'.")
+    TYPE_MAP = {
+        WaveSpectrum: "DirectionalSpectrum",
+        WaveBinSpectrum: "DirectionalBinSpectrum",
+    }
 
-    rao_squared = (rao * rao.conjugate()).real
-    rao_squared = rao_squared.reshape(freq, dirs, freq_hz=False, degrees=False)
-    wave_body = wave_body.reshape(freq, dirs, freq_hz=False, degrees=False)
+    try:
+        type_ = TYPE_MAP[type(wave)]
+    except KeyError:
+        raise ValueError(
+            "Invalid `wave` type. Should be 'WaveSpectrum' or 'WaveBinSpectrum'."
+        )
 
-    return multiply(rao_squared, wave_body, output_type="directional_spectrum")
+    return multiply(rao_squared, wave_body, output_type=type_)
 
 
 class BaseSpreading(ABC):
diff --git a/tests/test_core.py b/tests/test_core.py
index fd82005d..566d4172 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -5409,23 +5409,207 @@ def test_dirm_rad(self, d, s, mean_dir_rad):
 
 
 class Test_calculate_response:
+    @pytest.fixture
+    def rao(self):
+        freq = np.linspace(0.01, 2.0 * np.pi, 10, endpoint=False)  # rad/s
+        dirs = np.linspace(0.0, 2.0 * np.pi, 15, endpoint=False)  # rad
+
+        # random function: 2.0 * f + 3.0 * theta + 4.0
+        vals_amp = 2.0 * freq[:, np.newaxis] + 3.0 * dirs[np.newaxis, :] + 4.0
+        vals_phase = np.linspace(0, 2 * np.pi, len(freq) * len(dirs)).reshape(
+            len(freq), len(dirs)
+        )
+
+        rao = RAO.from_amp_phase(
+            freq,
+            dirs,
+            vals_amp,
+            vals_phase,
+            freq_hz=False,
+            degrees=False,
+            waves_coming_from=True,
+            clockwise=False,
+        )
+        return rao
+
+    @pytest.fixture
+    def wave(self):
+        freq = np.linspace(0.005, 2.0 * np.pi, 10, endpoint=False)  # rad/s
+        dirs = np.linspace(0.5, 2.0 * np.pi, 15, endpoint=False)  # rad
+
+        # random function: f + sqrt(theta) + 7.0
+        vals = freq[:, np.newaxis] ** 2 + np.sqrt(dirs[np.newaxis, :]) + 7.0
+
+        wave = WaveSpectrum(
+            freq,
+            dirs,
+            vals,
+            freq_hz=False,
+            degrees=False,
+            waves_coming_from=False,
+            clockwise=True,
+        )
+        return wave
+
+    @pytest.fixture
+    def wavebin(self):
+        freq = np.linspace(0.005, 2.0 * np.pi, 10, endpoint=False)  # rad/s
+        dirs = np.linspace(0.5, 2.0 * np.pi, 15, endpoint=False)  # rad
+
+        # random function: f + sqrt(theta) + 7.0
+        vals = freq[:, np.newaxis] ** 2 + np.sqrt(dirs[np.newaxis, :]) + 7.0
+
+        wave = WaveBinSpectrum(
+            freq,
+            dirs,
+            vals,
+            freq_hz=False,
+            degrees=False,
+            waves_coming_from=False,
+            clockwise=True,
+        )
+        return wave
+
     def test_calculate_response(self, rao, wave):
-        response = calculate_response(rao, wave, 0.0)
-        assert isinstance(response, DirectionalSpectrum)
-        assert response._clockwise == rao._clockwise
-        assert response._waves_coming_from == rao._waves_coming_from
+        response = calculate_response(rao, wave, np.radians(45.0))
+
+        # Expected response
+        wave_body = wave.rotate(45.0, degrees=True)
+        wave_body.set_wave_convention(waves_coming_from=True, clockwise=False)
+        freq_expect, dirs_expect = wave_body._freq, wave_body._dirs
+        rao_squared_expect = (rao * rao.conjugate()).real
+        rao_squared_expect = rao_squared_expect.reshape(
+            freq_expect, dirs_expect, freq_hz=False, degrees=False
+        )
+        response_expect = wr.multiply(
+            rao_squared_expect, wave_body, "DirectionalSpectrum"
+        )
+
+        assert isinstance(response, wr.DirectionalSpectrum)
+        assert response._clockwise is False
+        assert response._waves_coming_from is True
         assert response._freq_hz is False
         assert response._degrees is False
+        np.testing.assert_allclose(response._freq, response_expect._freq)
+        np.testing.assert_allclose(response._dirs, response_expect._dirs)
+        np.testing.assert_allclose(response._vals, response_expect._vals)
+
+    def test_calculate_response_bin(self, rao, wavebin):
+        response = calculate_response(rao, wavebin, np.radians(45.0))
+
+        # Expected response
+        wave_body = wavebin.rotate(45.0, degrees=True)
+        wave_body.set_wave_convention(waves_coming_from=True, clockwise=False)
+        freq_expect, dirs_expect = wave_body._freq, wave_body._dirs
+        rao_squared_expect = (rao * rao.conjugate()).real
+        rao_squared_expect = rao_squared_expect.reshape(
+            freq_expect, dirs_expect, freq_hz=False, degrees=False
+        )
+        response_expect = wr.multiply(
+            rao_squared_expect, wave_body, "DirectionalBinSpectrum"
+        )
+
+        assert isinstance(response, wr.DirectionalBinSpectrum)
+        assert response._clockwise is False
+        assert response._waves_coming_from is True
+        assert response._freq_hz is False
+        assert response._degrees is False
+        np.testing.assert_allclose(response._freq, response_expect._freq)
+        np.testing.assert_allclose(response._dirs, response_expect._dirs)
+        np.testing.assert_allclose(response._vals, response_expect._vals)
+
+    def test_calculate_response_heading_degrees(self, rao, wave):
+        response = calculate_response(rao, wave, 45, heading_degrees=True)
+
+        # Expected response
+        wave_body = wave.rotate(45.0, degrees=True)
+        wave_body.set_wave_convention(waves_coming_from=True, clockwise=False)
+        freq_expect, dirs_expect = wave_body._freq, wave_body._dirs
+        rao_squared_expect = (rao * rao.conjugate()).real
+        rao_squared_expect = rao_squared_expect.reshape(
+            freq_expect, dirs_expect, freq_hz=False, degrees=False
+        )
+        response_expect = wr.multiply(
+            rao_squared_expect, wave_body, "DirectionalSpectrum"
+        )
+
+        np.testing.assert_allclose(response._freq, response_expect._freq)
+        np.testing.assert_allclose(response._dirs, response_expect._dirs)
+        np.testing.assert_allclose(response._vals, response_expect._vals)
+
+    def test_calculate_response_heading_radians(self, rao, wave):
+        response = calculate_response(rao, wave, np.radians(45), heading_degrees=False)
+
+        # Expected response
+        wave_body = wave.rotate(45.0, degrees=True)
+        wave_body.set_wave_convention(waves_coming_from=True, clockwise=False)
+        freq_expect, dirs_expect = wave_body._freq, wave_body._dirs
+        rao_squared_expect = (rao * rao.conjugate()).real
+        rao_squared_expect = rao_squared_expect.reshape(
+            freq_expect, dirs_expect, freq_hz=False, degrees=False
+        )
+        response_expect = wr.multiply(
+            rao_squared_expect, wave_body, "DirectionalSpectrum"
+        )
+
+        np.testing.assert_allclose(response._freq, response_expect._freq)
+        np.testing.assert_allclose(response._dirs, response_expect._dirs)
+        np.testing.assert_allclose(response._vals, response_expect._vals)
+
+    def test_calculate_response_reshape_rao_squared(self, rao, wave):
+        response = calculate_response(rao, wave, np.radians(45), reshape="rao_squared")
+
+        # Expected response
+        wave_body = wave.rotate(45.0, degrees=True)
+        wave_body.set_wave_convention(waves_coming_from=True, clockwise=False)
+        freq_expect, dirs_expect = wave_body._freq, wave_body._dirs
+        rao_squared_expect = (rao * rao.conjugate()).real
+        rao_squared_expect = rao_squared_expect.reshape(
+            freq_expect, dirs_expect, freq_hz=False, degrees=False
+        )  # reshape squared RAO
+        response_expect = wr.multiply(
+            rao_squared_expect, wave_body, "DirectionalSpectrum"
+        )
+
+        np.testing.assert_allclose(response._freq, response_expect._freq)
+        np.testing.assert_allclose(response._dirs, response_expect._dirs)
+        np.testing.assert_allclose(response._vals, response_expect._vals)
+
+    def test_calculate_response_reshape_rao(self, rao, wave):
+        response = calculate_response(rao, wave, np.radians(45), reshape="rao")
+
+        # Expected response
+        wave_body = wave.rotate(45.0, degrees=True)
+        wave_body.set_wave_convention(waves_coming_from=True, clockwise=False)
+        freq_expect, dirs_expect = wave_body._freq, wave_body._dirs
+        rao_expect = rao.reshape(
+            freq_expect, dirs_expect, freq_hz=False, degrees=False
+        )  # reshape RAO
+        rao_squared_expect = (rao_expect * rao_expect.conjugate()).real
+        response_expect = wr.multiply(
+            rao_squared_expect, wave_body, "DirectionalSpectrum"
+        )
+
+        np.testing.assert_allclose(response._freq, response_expect._freq)
+        np.testing.assert_allclose(response._dirs, response_expect._dirs)
+        np.testing.assert_allclose(response._vals, response_expect._vals)
+
+    def test_calculate_response_reshape_raises(self, rao, wave):
+        with pytest.raises(ValueError):
+            calculate_response(rao, wave, np.radians(45), reshape="invalid-value")
 
     def test_calculate_response_raises_coord_freq(self, rao, wave):
+        # TODO: Deprecated functionality. Remove test in future.
         with pytest.raises(ValueError):
             calculate_response(rao, wave, 0.0, coord_freq="invalid-value")
 
     def test_calculate_response_raises_coord_dirs(self, rao, wave):
+        # TODO: Deprecated functionality. Remove test in future.
         with pytest.raises(ValueError):
             calculate_response(rao, wave, 0.0, coord_freq="invalid-value")
 
     def test_calculate_response_coord_wave(self):
+        # TODO: Deprecated functionality. Remove test in future.
         freq_rao = np.array([0.0, 0.5, 1.0])
         dirs_rao = np.array([45.0, 135.0, 225.0, 315.0])
         vals_rao = np.array(
@@ -5498,6 +5682,7 @@ def test_calculate_response_coord_wave(self):
         np.testing.assert_array_almost_equal(response._vals, vals_expect)
 
     def test_calculate_response_coord_rao(self):
+        # TODO: Deprecated functionality. Remove test in future.
         freq_rao = np.array([0.0, 0.5, 1.0])
         dirs_rao = np.array([45.0, 135.0, 225.0, 315.0])
         vals_rao = np.array(
@@ -5558,127 +5743,6 @@ def test_calculate_response_coord_rao(self):
         assert response._clockwise == rao._clockwise
         assert response._waves_coming_from == rao._waves_coming_from
 
-    def test_calculate_response_heading_degrees(self):
-        freq_rao = np.array([0.0, 0.5, 1.0])
-        dirs_rao = np.array([45.0, 135.0, 225.0, 315.0])
-        vals_rao = np.array(
-            [
-                [1.0 + 0.0j, 0.0 + 1.0j, 1.0 + 0.0j, 0.0 + 1.0j],
-                [1.0 + 0.0j, 0.0 + 1.0j, 1.0 + 0.0j, 0.0 + 1.0j],
-                [1.0 + 0.0j, 0.0 + 1.0j, 1.0 + 0.0j, 0.0 + 1.0j],
-            ]
-        )  # all amplitudes are 1
-        rao = RAO(
-            freq_rao,
-            dirs_rao,
-            vals_rao,
-            freq_hz=True,
-            degrees=True,
-        )
-
-        freq_wave = np.array([0.0, 0.3, 0.6, 0.9])  # extrapolation needed
-        dirs_wave = np.array([90.0, 180.0, 270.0, 359.0])
-        vals_wave = np.ones((len(freq_wave), len(dirs_wave)))
-        wave = WaveSpectrum(
-            freq_wave,
-            dirs_wave,
-            vals_wave,
-            freq_hz=True,
-            degrees=True,
-        )
-
-        response = calculate_response(
-            rao, wave, 45.0, heading_degrees=True, coord_freq="wave", coord_dirs="wave"
-        )
-        response.set_wave_convention(**wave.wave_convention)
-
-        freq_expect = wave._freq
-        dirs_expect = wave._dirs - (np.pi / 4.0)
-        vals_expect = (
-            1.0
-            / (2.0 * np.pi * np.pi / 180.0)
-            * np.array(
-                [
-                    [1.0, 1.0, 1.0, 1.0],
-                    [1.0, 1.0, 1.0, 1.0],
-                    [1.0, 1.0, 1.0, 1.0],
-                    [1.0, 1.0, 1.0, 1.0],
-                ]
-            )
-        )
-
-        assert isinstance(response, DirectionalSpectrum)
-        assert response._freq_hz is False
-        assert response._degrees is False
-        np.testing.assert_array_almost_equal(response._freq, freq_expect)
-        np.testing.assert_array_almost_equal(response._dirs, dirs_expect)
-        np.testing.assert_array_almost_equal(response._vals, vals_expect)
-        assert response._clockwise == rao._clockwise
-        assert response._waves_coming_from == rao._waves_coming_from
-
-    def test_calculate_response_heading_radians(self):
-        freq_rao = np.array([0.0, 0.5, 1.0])
-        dirs_rao = np.array([45.0, 135.0, 225.0, 315.0])
-        vals_rao = np.array(
-            [
-                [1.0 + 0.0j, 0.0 + 1.0j, 1.0 + 0.0j, 0.0 + 1.0j],
-                [1.0 + 0.0j, 0.0 + 1.0j, 1.0 + 0.0j, 0.0 + 1.0j],
-                [1.0 + 0.0j, 0.0 + 1.0j, 1.0 + 0.0j, 0.0 + 1.0j],
-            ]
-        )  # all amplitudes are 1
-        rao = RAO(
-            freq_rao,
-            dirs_rao,
-            vals_rao,
-            freq_hz=True,
-            degrees=True,
-        )
-
-        freq_wave = np.array([0.0, 0.3, 0.6, 0.9])  # extrapolation needed
-        dirs_wave = np.array([90.0, 180.0, 270.0, 359.0])
-        vals_wave = np.ones((len(freq_wave), len(dirs_wave)))
-        wave = WaveSpectrum(
-            freq_wave,
-            dirs_wave,
-            vals_wave,
-            freq_hz=True,
-            degrees=True,
-        )
-
-        response = calculate_response(
-            rao,
-            wave,
-            np.pi / 4.0,
-            heading_degrees=False,
-            coord_freq="wave",
-            coord_dirs="wave",
-        )
-        response.set_wave_convention(**wave.wave_convention)
-
-        freq_expect = wave._freq
-        dirs_expect = wave._dirs - (np.pi / 4.0)
-        vals_expect = (
-            1.0
-            / (2.0 * np.pi * np.pi / 180.0)
-            * np.array(
-                [
-                    [1.0, 1.0, 1.0, 1.0],
-                    [1.0, 1.0, 1.0, 1.0],
-                    [1.0, 1.0, 1.0, 1.0],
-                    [1.0, 1.0, 1.0, 1.0],
-                ]
-            )
-        )
-
-        assert isinstance(response, DirectionalSpectrum)
-        assert response._freq_hz is False
-        assert response._degrees is False
-        np.testing.assert_array_almost_equal(response._freq, freq_expect)
-        np.testing.assert_array_almost_equal(response._dirs, dirs_expect)
-        np.testing.assert_array_almost_equal(response._vals, vals_expect)
-        assert response._clockwise == rao._clockwise
-        assert response._waves_coming_from == rao._waves_coming_from
-
 
 class Test__check_is_similar:
     def test_check_is_similar(self):