diff --git a/docs/bibliography.bib b/docs/bibliography.bib
index bb8e39886..7211107de 100644
--- a/docs/bibliography.bib
+++ b/docs/bibliography.bib
@@ -48,3 +48,16 @@ @misc{kahan:2000:CPR
   note =         "Lecture notes for Math 128.",
   URL =          "http://www.cs.berkeley.edu/~wkahan/Math128/Cross.pdf",
 }
+
+@article{STAHN201644,
+    title = {Focusing neutron reflectometry: Implementation and experience on the TOF-reflectometer Amor},
+    journal = {Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment},
+    volume = {821},
+    pages = {44-54},
+    year = {2016},
+    issn = {0168-9002},
+    doi = {https://doi.org/10.1016/j.nima.2016.03.007},
+    url = {https://www.sciencedirect.com/science/article/pii/S0168900216300250},
+    author = {J. Stahn and A. Glavic},
+    keywords = {Reflectometry, Neutron, Focusing},
+}
diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py
index 367f48174..65019e3e0 100644
--- a/src/scippneutron/conversion/beamline.py
+++ b/src/scippneutron/conversion/beamline.py
@@ -104,6 +104,8 @@
 
 These definitions assume that gravity can be neglected.
 See :func:`scattering_angles_with_gravity` for definitions that account for gravity.
+And :func:`scattering_angle_in_yz_plane` for the definition used in reflectometry ---
+which also includes gravity.
 """
 
 from typing import TypedDict
@@ -475,21 +477,22 @@ def scattering_angles_with_gravity(
 
     .. math::
 
-        \mathsf{tan}(2\theta) &= \frac{\sqrt{x_d^2 + y_d^{\prime\, 2}}}{z} \\
-        \mathsf{tan}(\phi) &= \frac{y'_d}{x}
+        \mathsf{tan}(2\theta) &= \frac{\sqrt{x_d^2 + y_d^{\prime\, 2}}}{z_d} \\
+        \mathsf{tan}(\phi) &= \frac{y'_d}{x_d}
 
     Attention
     ---------
-    The above equation for :math:`y'_d` contains :math:`L_2^{\prime\, 2} = |b'_2|`
-    which in turn depends on :math:`y'_d`.
-    Solving this equation for :math:`y'_d` is too difficult.
-    Instead, we approximate :math:`L'_2 \approx L_2`.
-    The impact of this approximation on :math:`2\theta` is of the order of
-    :math:`10^{-6}` or less for beamlines at ESS.
-    This is within the expected statistical uncertainties and can be ignored.
-
-    See the file ``docs/user-guide/auxiliary/gravity_approx/gravity_approx.typ``
-    for details.
+        The above equation for :math:`y'_d` contains :math:`L_2^{\prime\, 2} = |b'_2|`
+        which in turn depends on :math:`y'_d`.
+        Solving this equation for :math:`y'_d` is too difficult.
+        Instead, we approximate :math:`L'_2 \approx L_2`.
+        The impact of this approximation on :math:`2\theta` is of the order of
+        :math:`10^{-6}` or less for beamlines at ESS.
+        This is within the expected statistical uncertainties and can be ignored.
+
+        See `two_theta gravity correction
+        <../../user-guide/algorithms-background/two_theta-gravity-correction.rst>`_
+        for details.
 
     Parameters
     ----------
@@ -507,16 +510,19 @@ def scattering_angles_with_gravity(
     :
         A dict containing the polar scattering angle ``'two_theta'`` and
         the azimuthal angle ``'phi'``.
+
+    See also
+    --------
+    scattering_angle_in_yz_plane:
+        Ignores the ``x`` component when computing ``theta``.
+        This is used in reflectometry.
     """
-    match beam_aligned_unit_vectors(incident_beam=incident_beam, gravity=gravity):
-        case {
-            'beam_aligned_unit_x': ex,
-            'beam_aligned_unit_y': ey,
-            'beam_aligned_unit_z': ez,
-        }:
-            pass
-        case _:
-            raise RuntimeError('Unexpected return value of beam_aligned_unit_vectors')
+    unit_vectors = beam_aligned_unit_vectors(
+        incident_beam=incident_beam, gravity=gravity
+    )
+    ex = unit_vectors['beam_aligned_unit_x']
+    ey = unit_vectors['beam_aligned_unit_y']
+    ez = unit_vectors['beam_aligned_unit_z']
 
     y = _drop_due_to_gravity(
         distance=sc.norm(scattered_beam), wavelength=wavelength, gravity=gravity
@@ -536,3 +542,83 @@ def scattering_angles_with_gravity(
     two_theta_ = sc.atan2(y=y, x=z, out=y)
 
     return {'two_theta': two_theta_, 'phi': phi}
+
+
+def scattering_angle_in_yz_plane(
+    incident_beam: sc.Variable,
+    scattered_beam: sc.Variable,
+    wavelength: sc.Variable,
+    gravity: sc.Variable,
+) -> sc.Variable:
+    r"""Compute polar scattering angles in the y-z plane using gravity.
+
+    Note
+    ----
+        This function uses the reflectometry definition of the polar scattering angle.
+        Other techniques define the angle w.r.t. the incident beam.
+        See :func:`scattering_angles_with_gravity` for those use cases.
+
+    With the definitions given in :func:`scattering_angles_with_gravity`,
+    and ignoring :math:`x_d`, we get
+
+    .. math::
+
+        \mathsf{tan}(\gamma) = \frac{|y_d^{\prime}|}{z_d}
+
+    with
+
+    .. math::
+
+        y'_d = y_d + \frac{|g| m_n^2}{2 h^2} L_2^{\prime\, 2} \lambda^2
+
+    The angle :math:`\gamma` is defined as in Fig. 5 of :cite:`STAHN201644`.
+
+    Attention
+    ---------
+        The above equation for :math:`y'_d` contains :math:`L_2^{\prime\, 2} = |b'_2|`
+        which in turn depends on :math:`y'_d`.
+        Solving this equation for :math:`y'_d` is too difficult.
+        Instead, we approximate :math:`L'_2 \approx L_2`.
+        The impact of this approximation on :math:`\gamma` is of the order of
+        :math:`10^{-6}` or less for beamlines at ESS.
+        This is within the expected statistical uncertainties and can be ignored.
+
+        See `two_theta gravity correction
+        <../../user-guide/algorithms-background/two_theta-gravity-correction.rst>`_
+        for details.
+
+    Parameters
+    ----------
+    incident_beam:
+        Beam from source to sample. Expects ``dtype=vector3``.
+    scattered_beam:
+        Beam from sample to detector. Expects ``dtype=vector3``.
+    wavelength:
+        Wavelength of neutrons.
+    gravity:
+        Gravity vector.
+
+    Returns
+    -------
+    :
+        The polar scattering angle :math:`\gamma`.
+
+    See also
+    --------
+    scattering_angles_with_gravity:
+        Includes the ``x`` component when computing ``theta``.
+        This is used in techniques other than reflectometry.
+    """
+    unit_vectors = beam_aligned_unit_vectors(
+        incident_beam=incident_beam, gravity=gravity
+    )
+    ey = unit_vectors['beam_aligned_unit_y']
+    ez = unit_vectors['beam_aligned_unit_z']
+
+    y = _drop_due_to_gravity(
+        distance=sc.norm(scattered_beam), wavelength=wavelength, gravity=gravity
+    )
+    y += sc.dot(scattered_beam, ey).to(dtype=elem_dtype(wavelength), copy=False)
+    y = sc.abs(y, out=y)
+    z = sc.dot(scattered_beam, ez).to(dtype=elem_dtype(y), copy=False)
+    return sc.atan2(y=y, x=z, out=y)
diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py
index bbeabbdaa..0039974ae 100644
--- a/tests/conversion/beamline_conversions_test.py
+++ b/tests/conversion/beamline_conversions_test.py
@@ -544,6 +544,320 @@ def test_scattering_angles_with_gravity_supports_mismatching_units():
     sc.testing.assert_allclose(res['phi'], expected['phi'])
 
 
+def test_scattering_angle_in_yz_plane_requires_gravity_orthogonal_to_incident_beam():
+    incident_beam = sc.vector([0.564, 1.2, -10.4], unit='m')
+    scattered_beam = sc.vectors(
+        dims=['beam'], values=[[13, 24, 35], [51, -42, 33]], unit='m'
+    )
+    wavelength = sc.array(dims=['wavelength'], values=[1.2, 1.6, 1.8], unit='Å')
+    gravity = sc.vector([0, 0, sc.constants.g.value], unit=sc.constants.g.unit)
+
+    with pytest.raises(
+        ValueError, match='`gravity` and `incident_beam` must be orthogonal'
+    ):
+        beamline.scattering_angle_in_yz_plane(
+            incident_beam=incident_beam,
+            scattered_beam=scattered_beam,
+            wavelength=wavelength,
+            gravity=gravity,
+        )
+
+
+def test_scattering_angle_in_yz_plane_small_gravity():
+    # This case is unphysical but tests the consistency with `two_theta`.
+    # Note that the scattered beam must be in the x-z plane for `two_theta`
+    # and `scattering_angle_from_sample` to compute the same angle.
+    incident_beam = sc.vector([0.0, 0.0, 10.4], unit='m')
+    scattered_beam = sc.vectors(
+        dims=['beam'],
+        values=[[0, 24, 0], [0, -42, 0]],
+        unit='m',
+    )
+    wavelength = sc.array(dims=['wavelength'], values=[1.2, 1.6, 1.8], unit='Å')
+    gravity = sc.vector([0, -1e-11, 0], unit=sc.constants.g.unit)
+
+    res = beamline.scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+    expected = beamline.two_theta(
+        incident_beam=incident_beam, scattered_beam=scattered_beam
+    ).broadcast(dims=['wavelength', 'beam'], shape=[3, 2])
+    sc.testing.assert_allclose(res, expected)
+
+
+@pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi])
+def test_scattering_angle_in_yz_plane_reproduces_polar_angle(
+    polar: float,
+):
+    # This case is unphysical but tests that the function reproduces
+    # the expected angles using a rotated vector.
+
+    gravity = sc.vector([0.0, -1e-11, 0.0], unit='cm/s^2')
+    incident_beam = sc.vector([0.0, 0.0, 968.0], unit='cm')
+
+    # With this definition, the x-axis has azimuthal=0.
+    rot1 = sc.spatial.rotations_from_rotvecs(sc.vector([-polar, 0, 0], unit='rad'))
+    scattered_beam = rot1 * incident_beam
+
+    wavelength = sc.scalar(1e-6, unit='Å')
+
+    res = beamline.scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+
+    sc.testing.assert_allclose(res, sc.scalar(polar, unit='rad'))
+
+
+@pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi])
+@pytest.mark.parametrize('x', [0.5, 11.4, -9.7])
+def test_scattering_angle_in_yz_plane_does_not_depend_on_x(polar: float, x: float):
+    # This case is unphysical but tests that the function reproduces
+    # the expected angles using a rotated vector.
+
+    gravity = sc.vector([0.0, -1e-11, 0.0], unit='cm/s^2')
+    incident_beam = sc.vector([0.0, 0.0, 968.0], unit='cm')
+
+    # With this definition, the x-axis has azimuthal=0.
+    rot1 = sc.spatial.rotations_from_rotvecs(sc.vector([-polar, 0, 0], unit='rad'))
+    scattered_beam_ref = rot1 * incident_beam
+    scattered_beam_shift = scattered_beam_ref + sc.vector([x, 0.0, 0.0], unit='cm')
+
+    wavelength = sc.scalar(1e-6, unit='Å')
+
+    res_shift = beamline.scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam_shift,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+    res_ref = beamline.scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam_ref,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+
+    sc.testing.assert_allclose(res_shift, res_ref)
+
+
+def test_scattering_angle_in_yz_plane_drops_in_expected_direction():
+    wavelength = sc.scalar(1.6, unit='Å')
+    gravity = sc.vector([0.0, -sc.constants.g.value, 0.0], unit=sc.constants.g.unit)
+    incident_beam = sc.vector([0.0, 0.0, 41.1], unit='m')
+    scattered_beam = sc.vectors(
+        dims=['det'], values=[[0.0, 2.5, 8.6], [0.0, -1.7, 6.9]], unit='m'
+    )
+
+    with_gravity = beamline.scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+    without_gravity = beamline.two_theta(
+        incident_beam=incident_beam, scattered_beam=scattered_beam
+    )
+
+    # The neutron was detected above the incident beam.
+    # So using straight paths, it looks like it scattered at a
+    # smaller angle (detected at a lower y) than in the real case with gravity.
+    assert sc.all(with_gravity[0] > without_gravity[0]).value
+    # The neutron was detected below the incident beam.
+    # So the opposite of the above comment applies.
+    assert sc.all(with_gravity[1] < without_gravity[1]).value
+
+
+def test_scattering_angle_in_yz_plane_beams_aligned_with_lab_coords():
+    wavelength = sc.array(dims=['wavelength'], values=[1.6, 0.9, 0.7], unit='Å')
+    gravity = sc.vector([0.0, -sc.constants.g.value, 0.0], unit=sc.constants.g.unit)
+    incident_beam = sc.vector([0.0, 0.0, 41.1], unit='m')
+    scattered_beam = sc.vectors(
+        dims=['det'], values=[[0.0, 2.5, 3.6], [0.0, -1.7, 2.9]], unit='m'
+    )
+
+    original_wavelength = wavelength.copy()
+    original_gravity = gravity.copy()
+    original_incident_beam = incident_beam.copy()
+    original_scattered_beam = scattered_beam.copy()
+
+    res = beamline.scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+
+    L2 = sc.norm(scattered_beam)
+    x = scattered_beam.fields.x
+    y = scattered_beam.fields.y
+    z = scattered_beam.fields.z
+    drop = (
+        L2**2
+        * sc.norm(gravity)
+        * wavelength**2
+        * sc.constants.m_n**2
+        / (2 * sc.constants.h**2)
+    )
+    drop = drop.to(unit=y.unit)
+    true_y = y + drop
+    true_scattered_beam = sc.spatial.as_vectors(x, true_y, z)
+    expected_theta = sc.asin(abs(true_y) / sc.norm(true_scattered_beam))
+    expected_theta = expected_theta.transpose(res.dims)
+
+    sc.testing.assert_allclose(res, expected_theta)
+
+    sc.testing.assert_identical(wavelength, original_wavelength)
+    sc.testing.assert_identical(gravity, original_gravity)
+    sc.testing.assert_identical(incident_beam, original_incident_beam)
+    sc.testing.assert_identical(scattered_beam, original_scattered_beam)
+
+
+def _reference_scattering_angle_in_yz_plane(
+    incident_beam: sc.Variable,
+    scattered_beam: sc.Variable,
+    gravity: sc.Variable,
+    wavelength: sc.Variable,
+) -> sc.Variable:
+    # This is a simplified, independently checked implementation.
+    e_z = incident_beam / sc.norm(incident_beam)
+    e_y = -gravity / sc.norm(gravity)
+
+    y = sc.dot(scattered_beam, e_y)
+    z = sc.dot(scattered_beam, e_z)
+
+    L2 = sc.norm(scattered_beam)
+    drop = (
+        L2**2
+        * wavelength**2
+        * sc.norm(gravity)
+        * (sc.constants.m_n**2 / (2 * sc.constants.h**2))
+    )
+    dropped_y = y + drop.to(unit=y.unit)
+
+    return sc.atan2(y=abs(dropped_y), x=z)
+
+
+def test_scattering_angle_in_yz_plane_beams_unaligned_with_lab_coords():
+    wavelength = sc.array(dims=['wavelength'], values=[1.6, 0.9, 0.7], unit='Å')
+    # Gravity and incident_beam are not aligned with the coordinate system
+    # but orthogonal to each other.
+    gravity = sc.vector([-0.3, -9.81, 0.01167883211678832], unit='m/s^2')
+    incident_beam = sc.vector([1.6, 0.0, 41.1], unit='m')
+    scattered_beam = sc.vectors(
+        dims=['det'], values=[[1.8, 2.5, 3.6], [-0.4, -1.7, 2.9]], unit='m'
+    )
+
+    original_wavelength = wavelength.copy()
+    original_gravity = gravity.copy()
+    original_incident_beam = incident_beam.copy()
+    original_scattered_beam = scattered_beam.copy()
+
+    res = beamline.scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+    expected = _reference_scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+
+    sc.testing.assert_allclose(res, expected.transpose(res.dims))
+
+    sc.testing.assert_identical(wavelength, original_wavelength)
+    sc.testing.assert_identical(gravity, original_gravity)
+    sc.testing.assert_identical(incident_beam, original_incident_beam)
+    sc.testing.assert_identical(scattered_beam, original_scattered_beam)
+
+
+def test_scattering_angle_in_yz_plane_binned_data():
+    wavelength = sc.array(dims=['wavelength'], values=[1.6, 0.9, 0.7], unit='Å')
+    wavelength = sc.bins(
+        dim='wavelength',
+        data=wavelength,
+        begin=sc.array(dims=['det'], values=[0, 2], unit=None),
+        end=sc.array(dims=['det'], values=[2, 3], unit=None),
+    )
+    gravity = sc.vector([-0.3, -9.81, 0.01167883211678832], unit='m/s^2')
+    incident_beam = sc.vector([1.6, 0.0, 41.1], unit='m')
+    scattered_beam = sc.vectors(
+        dims=['det'], values=[[1.8, 2.5, 3.6], [-0.4, -1.7, 2.9]], unit='m'
+    )
+
+    original_wavelength = wavelength.copy()
+    original_gravity = gravity.copy()
+    original_incident_beam = incident_beam.copy()
+    original_scattered_beam = scattered_beam.copy()
+
+    res = beamline.scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+    expected = _reference_scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+
+    sc.testing.assert_allclose(res, expected.transpose(res.dims))
+    sc.testing.assert_identical(wavelength, original_wavelength)
+    sc.testing.assert_identical(gravity, original_gravity)
+    sc.testing.assert_identical(incident_beam, original_incident_beam)
+    sc.testing.assert_identical(scattered_beam, original_scattered_beam)
+
+
+def test_scattering_angle_in_yz_plane_uses_wavelength_dtype():
+    wavelength = sc.array(
+        dims=['wavelength'], values=[1.6, 0.7], unit='Å', dtype='float32'
+    )
+    gravity = sc.vector([0.0, -9.81, 0.0], unit='m/s^2')
+    incident_beam = sc.vector([0.0, 0.0, 41.1], unit='m')
+    scattered_beam = sc.vector([1.8, 2.5, 3.6], unit='m')
+
+    res = beamline.scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+    assert res.dtype == 'float32'
+
+
+def test_scattering_angle_in_yz_plane_supports_mismatching_units():
+    wavelength = sc.array(dims=['wavelength'], values=[1.6, 0.7], unit='Å')
+    gravity = sc.vector([0.0, -9810, 0.0], unit='m/ms^2')
+    incident_beam = sc.vector([0.0, 0.0, 410], unit='cm')
+    scattered_beam = sc.vector([180, 1800, 2400], unit='mm')
+
+    res = beamline.scattering_angle_in_yz_plane(
+        incident_beam=incident_beam,
+        scattered_beam=scattered_beam,
+        wavelength=wavelength,
+        gravity=gravity,
+    )
+
+    expected = _reference_scattering_angle_in_yz_plane(
+        incident_beam=incident_beam.to(unit='m'),
+        scattered_beam=scattered_beam.to(unit='m'),
+        wavelength=wavelength,
+        gravity=gravity.to(unit='m/s^2'),
+    )
+
+    sc.testing.assert_allclose(res, expected)
+
+
 def test_beam_aligned_unit_vectors_axis_aligned_inputs():
     res = beamline.beam_aligned_unit_vectors(
         incident_beam=sc.vector([0.0, 0.0, 2.1], unit='mm'),