diff --git a/CHANGES.md b/CHANGES.md
index af7b8ec1..5dd28b26 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -7,6 +7,11 @@ issue tracker.  The initial project development roadmap is documented here:
 
 ## Unreleased
 
+### Fixed
+
+- `ReferenceConstraint.evaluate()` computes real residuals; verified post-solve. (#304)
+
+
 ## Release v0.11.2
 
 Released 2026-06-22
diff --git a/src/ad_hoc_diffractometer/forward.py b/src/ad_hoc_diffractometer/forward.py
index d40157f1..3bc6acda 100644
--- a/src/ad_hoc_diffractometer/forward.py
+++ b/src/ad_hoc_diffractometer/forward.py
@@ -4234,6 +4234,7 @@ def _validate_solutions(
     from .display import precision_atol
     from .mode import BisectConstraint
     from .mode import ConstraintViolation
+    from .mode import ReferenceConstraint
     from .mode import SampleConstraint
 
     if not solutions:
@@ -4243,10 +4244,15 @@ def _validate_solutions(
 
     for i, sol in enumerate(solutions):
         for c in mode.constraints:
-            # Only validate constraints that have a numeric residual we can check.
-            # BisectConstraint and SampleConstraint are the checkable ones;
-            # DetectorConstraint is applied by the solver directly (ttheta_deg);
-            # ReferenceConstraint is not yet implemented.
+            # Validate constraints that expose a numeric per-solution residual.
+            # BisectConstraint and SampleConstraint are checked against the
+            # motor angles directly; DetectorConstraint is applied by the
+            # solver (ttheta_deg) so it never appears here.
+            # ReferenceConstraint pseudo-angle residuals
+            # (incidence/emergence/specular/omega/naz) are checked too — except
+            # "psi", which is enforced upstream as a validation filter in
+            # _solve_psi_mode (its motor-frame residual is subject to ±360°
+            # azimuthal wraparound and is not a reliable per-solution check).
             if isinstance(c, BisectConstraint | SampleConstraint):
                 try:
                     residual = c.evaluate(sol, geometry)
@@ -4272,6 +4278,29 @@ def _validate_solutions(
                         residual=residual,
                         tolerance=atol,
                     )
+            elif isinstance(c, ReferenceConstraint):
+                # "psi" is enforced by the validation filter, not here.
+                if c.name == "psi":
+                    continue
+                # Skip when the required reference vector is unset (evaluate()
+                # would raise ValueError); is_implemented() already gated the
+                # solve, but a mode could in principle reach here without one.
+                if not c.has_reference_vector(geometry):
+                    continue
+                try:
+                    residual = c.evaluate(sol, geometry)
+                except (KeyError, ValueError):
+                    continue
+                if abs(residual) > atol:
+                    raise ConstraintViolation(
+                        solution_index=i,
+                        constraint_repr=(
+                            f"ReferenceConstraint({c.name!r}, {c.value!r}): "
+                            f"residual {residual!s}° exceeds tolerance"
+                        ),
+                        residual=residual,
+                        tolerance=atol,
+                    )
 
 
 def _check_limits(
diff --git a/src/ad_hoc_diffractometer/mode.py b/src/ad_hoc_diffractometer/mode.py
index 37de8597..916ac8f8 100644
--- a/src/ad_hoc_diffractometer/mode.py
+++ b/src/ad_hoc_diffractometer/mode.py
@@ -918,17 +918,68 @@ def evaluate(
         geometry: AdHocDiffractometer,
     ) -> float:
         """
-        Return the constraint residual (zero when satisfied).
+        Return the constraint residual in degrees (zero when satisfied).
+
+        The residual is the difference between the physical pseudo-angle
+        computed from the supplied motor angles (via the standalone
+        functions in :mod:`~ad_hoc_diffractometer.reference`) and the
+        target :attr:`value`:
+
+        - ``"incidence"`` — ``incidence_angle - value``
+        - ``"emergence"`` — ``emergence_angle - value``
+        - ``"specular"`` — ``incidence_angle - emergence_angle``
+          (the relational condition α_i = α_f; :attr:`value` is the
+          boolean flag ``True`` and does not enter the residual)
+        - ``"psi"`` — ``psi_angle - value``
+        - ``"naz"`` — ``naz_angle - value``
+        - ``"omega"`` — ``omega_pseudo - value``
+
+        Parameters
+        ----------
+        angles : dict[str, float]
+            Current motor angles in degrees, keyed by stage name.
+        geometry : AdHocDiffractometer
+            The diffractometer.  The required reference vector
+            (``surface_normal`` for incidence/emergence/specular,
+            ``azimuth`` for psi/naz) must be set, except for ``"omega"``
+            which needs no reference vector.
+
+        Returns
+        -------
+        float
+            Residual in degrees.  Zero means the constraint is satisfied.
 
-        Not yet implemented — all reference constraints require the
-        reference vector infrastructure (Issue J).  Raises
-        ``NotImplementedError`` when called.
+        Raises
+        ------
+        ValueError
+            If the required reference vector is not set on the geometry
+            (raised by the underlying
+            :mod:`~ad_hoc_diffractometer.reference` function).
         """
-        raise NotImplementedError(
-            f"ReferenceConstraint('{self._name}') solver is not yet implemented. "
-            "Reference constraint solvers require the reference vector "
-            "infrastructure (Issue J)."
-        )
+        # Local import to avoid module-load ordering issues (reference.py
+        # imports geometry types lazily; mirror the pattern used elsewhere
+        # in this package).
+        from . import reference
+
+        if self._name == "incidence":
+            return reference.incidence_angle(geometry, angles=angles) - float(
+                self._value
+            )
+        if self._name == "emergence":
+            return reference.emergence_angle(geometry, angles=angles) - float(
+                self._value
+            )
+        if self._name == "specular":
+            # Relational: incidence == emergence.  value is the flag True.
+            return reference.incidence_angle(
+                geometry, angles=angles
+            ) - reference.emergence_angle(geometry, angles=angles)
+        if self._name == "psi":
+            return reference.psi_angle(geometry, angles=angles) - float(self._value)
+        if self._name == "naz":
+            return reference.naz_angle(geometry, angles=angles) - float(self._value)
+        # omega
+        return reference.omega_pseudo(geometry, angles=angles) - float(self._value)
 
     def is_satisfied(
         self,
@@ -936,10 +987,8 @@ def is_satisfied(
         geometry: AdHocDiffractometer,
         tol: float = 1e-6,
     ) -> bool:
-        """Return True when the constraint is satisfied (not yet implemented)."""
-        raise NotImplementedError(
-            f"ReferenceConstraint('{self._name}') is not yet implemented."
-        )
+        """Return True when ``|evaluate(angles, geometry)| < tol``."""
+        return abs(self.evaluate(angles, geometry)) < tol
 
     def has_reference_vector(self, geometry: AdHocDiffractometer) -> bool:
         """
diff --git a/tests/test_mode.py b/tests/test_mode.py
index 88c549a8..1eaaf497 100644
--- a/tests/test_mode.py
+++ b/tests/test_mode.py
@@ -797,39 +797,145 @@ def test_reference_constraint_to_dict_from_dict_a_eq_b():
 
 
 @pytest.mark.parametrize(
-    "name, value, surface_normal, expected",
+    "name, value, surface_normal, azimuth, expected",
     [
         # incidence/emergence/specular: implemented when surface_normal is set
-        pytest.param("incidence", 0.0, (0, 0, 1), True, id="incidence-with-sn"),
-        pytest.param("incidence", 0.0, None, False, id="incidence-no-sn"),
-        pytest.param("emergence", 0.0, (0, 0, 1), True, id="emergence-with-sn"),
-        pytest.param("specular", True, (0, 0, 1), True, id="specular-with-sn"),
-        pytest.param("specular", True, None, False, id="specular-no-sn"),
-        # psi/naz: never implemented (no forward solver yet)
-        pytest.param("psi", 0.0, (0, 0, 1), False, id="psi-not-implemented"),
-        pytest.param("naz", 0.0, (0, 0, 1), False, id="naz-not-implemented"),
+        pytest.param("incidence", 0.0, (0, 0, 1), None, True, id="incidence-with-sn"),
+        pytest.param("incidence", 0.0, None, None, False, id="incidence-no-sn"),
+        pytest.param("emergence", 0.0, (0, 0, 1), None, True, id="emergence-with-sn"),
+        pytest.param("specular", True, (0, 0, 1), None, True, id="specular-with-sn"),
+        pytest.param("specular", True, None, None, False, id="specular-no-sn"),
+        # psi: implemented (validation filter) when azimuth is set
+        pytest.param("psi", 0.0, None, (0, 0, 1), True, id="psi-with-azimuth"),
+        pytest.param("psi", 0.0, None, None, False, id="psi-no-azimuth"),
+        # naz: no forward solver yet
+        pytest.param("naz", 0.0, (0, 0, 1), None, False, id="naz-not-implemented"),
     ],
 )
-def test_reference_constraint_is_implemented(name, value, surface_normal, expected):
-    """ReferenceConstraint.is_implemented() reflects surface_normal presence and solver."""
+def test_reference_constraint_is_implemented(
+    name, value, surface_normal, azimuth, expected
+):
+    """ReferenceConstraint.is_implemented() reflects reference vector and solver."""
     g = _psic_like()
     g._surface_normal = surface_normal  # noqa: SLF001
+    g._azimuth = azimuth  # noqa: SLF001
     rc = ReferenceConstraint(name, value)
     assert rc.is_implemented(g) is expected
 
 
-def test_reference_constraint_evaluate_raises():
-    g = _psic_like()
-    rc = ReferenceConstraint("psi", 90.0)
-    with pytest.raises(NotImplementedError, match="not yet implemented"):
-        rc.evaluate({}, g)
+def _psic_reference_geometry(**kwargs):
+    """Real psic geometry with wavelength, cubic UB, and a reference vector."""
+    import ad_hoc_diffractometer as ahd
+    from ad_hoc_diffractometer import ub_identity
 
+    g = ahd.make_geometry("psic")
+    g.wavelength = 1.54
+    g.sample.lattice = ahd.Lattice(a=4.0)
+    ub_identity(g.sample)
+    for key, val in kwargs.items():
+        setattr(g, key, val)
+    return g
 
-def test_reference_constraint_is_satisfied_raises():
-    g = _psic_like()
-    rc = ReferenceConstraint("psi", 90.0)
-    with pytest.raises(NotImplementedError):
-        rc.is_satisfied({}, g)
+
+@pytest.mark.parametrize(
+    "name, value, kwargs, context",
+    [
+        pytest.param(
+            "incidence",
+            0.0,
+            {"surface_normal": (0, 0, 1)},
+            does_not_raise(),
+            id="incidence-residual",
+        ),
+        pytest.param(
+            "emergence",
+            0.0,
+            {"surface_normal": (0, 0, 1)},
+            does_not_raise(),
+            id="emergence-residual",
+        ),
+        pytest.param(
+            "specular",
+            True,
+            {"surface_normal": (0, 0, 1)},
+            does_not_raise(),
+            id="specular-residual",
+        ),
+        pytest.param(
+            "omega",
+            0.0,
+            {},
+            does_not_raise(),
+            id="omega-residual",
+        ),
+        pytest.param(
+            "incidence",
+            0.0,
+            {},
+            pytest.raises(ValueError, match=re.escape("surface_normal must be set")),
+            id="incidence-no-surface-normal-raises",
+        ),
+    ],
+)
+def test_reference_constraint_evaluate(name, value, kwargs, context):
+    """ReferenceConstraint.evaluate() returns a real pseudo-angle residual.
+
+    A solution returned by forward() satisfies its reference constraint, so
+    the residual evaluated at that solution must be ~0.  evaluate() raises
+    ValueError when the required reference vector is not set.
+    """
+    g = _psic_reference_geometry(**kwargs)
+    rc = ReferenceConstraint(name, value)
+    with context:
+        # Use the geometry's current (zeroed) angles for the raising case;
+        # for the satisfied cases use a real forward solution.
+        angles = {s.name: s.angle for s in g._stages.values()}  # noqa: SLF001
+        if "surface_normal" in kwargs or name == "omega":
+            g.mode_name = {
+                "incidence": "fixed_incidence_vertical",
+                "emergence": "fixed_emergence_vertical",
+                "specular": "specular_vertical",
+                "omega": "fixed_omega_vertical",
+            }[name]
+            sols = g.forward(1, 0, 1)
+            angles = sols[0]
+        residual = rc.evaluate(angles, g)
+        assert residual == pytest.approx(0.0, abs=1e-6)
+
+
+def test_reference_constraint_is_satisfied():
+    """is_satisfied() is True at a satisfying solution, False for a wrong target."""
+    g = _psic_reference_geometry(surface_normal=(0, 0, 1))
+    g.mode_name = "fixed_incidence_vertical"
+    sol = g.forward(1, 0, 1)[0]
+    assert ReferenceConstraint("incidence", 0.0).is_satisfied(sol, g) is True
+    assert ReferenceConstraint("incidence", 99.0).is_satisfied(sol, g) is False
+
+
+def test_reference_constraint_evaluate_psi_branch():
+    """evaluate() for the 'psi' branch returns psi_angle - target."""
+    import ad_hoc_diffractometer.reference as reference
+
+    g = _psic_reference_geometry(azimuth=(0, 0, 1))
+    angles = {s.name: s.angle for s in g._stages.values()}  # noqa: SLF001
+    angles["delta"] = 30.0  # ensure a non-degenerate Q for psi
+    angles["eta"] = 15.0
+    expected = reference.psi_angle(g, angles=angles) - 12.0
+    assert ReferenceConstraint("psi", 12.0).evaluate(angles, g) == pytest.approx(
+        expected
+    )
+
+
+def test_reference_constraint_evaluate_naz_branch():
+    """evaluate() for the 'naz' branch returns naz_angle - target."""
+    import ad_hoc_diffractometer.reference as reference
+
+    g = _psic_reference_geometry(surface_normal=(1, 0, 0))
+    angles = {s.name: s.angle for s in g._stages.values()}  # noqa: SLF001
+    expected = reference.naz_angle(g, angles=angles) - 7.0
+    assert ReferenceConstraint("naz", 7.0).evaluate(angles, g) == pytest.approx(
+        expected
+    )
 
 
 # ---------------------------------------------------------------------------
diff --git a/tests/test_regression_issue_304.py b/tests/test_regression_issue_304.py
new file mode 100644
index 00000000..7c5a4505
--- /dev/null
+++ b/tests/test_regression_issue_304.py
@@ -0,0 +1,248 @@
+# Copyright (c) 2026-2026 UChicago Argonne, LLC
+# SPDX-License-Identifier: LicenseRef-UChicago-Argonne-LLC-License
+"""
+Regression tests for issue #304.
+
+Issue #304 reported that ``ReferenceConstraint.evaluate()`` /
+``.is_satisfied()`` were stubs raising ``NotImplementedError``, with code
+comments pointing at the long-closed "Issue J" (#157).  In fact the
+forward solvers for every reference-constraint mode
+(``incidence`` / ``emergence`` / ``specular`` / ``psi`` / ``omega``)
+were already implemented in :mod:`ad_hoc_diffractometer.forward`; the
+dispatcher routes those modes to dedicated ``_solve_*`` functions and
+never called the stubbed methods.
+
+The fix makes ``ReferenceConstraint.evaluate()`` /
+``.is_satisfied()`` compute real pseudo-angle residuals using
+:mod:`ad_hoc_diffractometer.reference`, and wires the reference
+constraints (except ``psi``, which is enforced upstream as a validation
+filter) into the post-solve verification loop in
+:func:`ad_hoc_diffractometer.forward._validate_solutions`.
+
+These tests span the ``mode``, ``forward``, and ``reference`` modules,
+so they live in a cross-module regression file per the project's testing
+conventions.
+
+The tests below confirm:
+
+* every affected reference-constraint mode solves end-to-end through
+  ``forward()`` once its reference vector is set (no spurious
+  ``ConstraintViolation`` from the verification loop), and
+* the verification loop rejects a solution whose reference residual
+  exceeds tolerance.
+"""
+
+from __future__ import annotations
+
+from contextlib import nullcontext as does_not_raise
+
+import pytest
+
+import ad_hoc_diffractometer as ahd
+from ad_hoc_diffractometer import ConstraintViolation
+from ad_hoc_diffractometer import ReferenceConstraint
+from ad_hoc_diffractometer.forward import _validate_solutions
+
+WAVELENGTH = 1.5406  # Cu Kα
+
+
+def _setup_cubic(name, a=4.0, **kwargs):
+    """Return a cubic geometry ready for forward()."""
+    g = ahd.make_geometry(name)
+    g.wavelength = WAVELENGTH
+    g.sample.lattice = ahd.Lattice(a=a)
+    ahd.ub_identity(g.sample)
+    for key, val in kwargs.items():
+        setattr(g, key, val)
+    return g
+
+
+# ---------------------------------------------------------------------------
+# Every affected reference-constraint mode solves end-to-end.
+# ---------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize(
+    "geometry, mode_name, ref_attr, ref_val, context",
+    [
+        pytest.param(
+            "psic",
+            "fixed_incidence_vertical",
+            "surface_normal",
+            (0, 0, 1),
+            does_not_raise(),
+            id="psic-incidence-vertical",
+        ),
+        pytest.param(
+            "psic",
+            "fixed_emergence_vertical",
+            "surface_normal",
+            (0, 0, 1),
+            does_not_raise(),
+            id="psic-emergence-vertical",
+        ),
+        pytest.param(
+            "psic",
+            "specular_vertical",
+            "surface_normal",
+            (0, 0, 1),
+            does_not_raise(),
+            id="psic-specular-vertical",
+        ),
+        pytest.param(
+            "psic",
+            "specular_horizontal",
+            "surface_normal",
+            (0, 0, 1),
+            does_not_raise(),
+            id="psic-specular-horizontal",
+        ),
+        pytest.param(
+            "psic",
+            "fixed_omega_vertical",
+            "surface_normal",
+            (0, 0, 1),
+            does_not_raise(),
+            id="psic-omega-vertical",
+        ),
+        pytest.param(
+            "sixc",
+            "fixed_incidence_zaxis",
+            "surface_normal",
+            (0, 0, 1),
+            does_not_raise(),
+            id="sixc-incidence-zaxis",
+        ),
+        pytest.param(
+            "sixc",
+            "specular_zaxis",
+            "surface_normal",
+            (0, 0, 1),
+            does_not_raise(),
+            id="sixc-specular-zaxis",
+        ),
+        pytest.param(
+            "zaxis",
+            "zaxis",
+            "surface_normal",
+            (0, 0, 1),
+            does_not_raise(),
+            id="zaxis-zaxis",
+        ),
+        pytest.param(
+            "zaxis",
+            "reflectivity",
+            "surface_normal",
+            (0, 0, 1),
+            does_not_raise(),
+            id="zaxis-reflectivity",
+        ),
+        pytest.param(
+            "s2d2",
+            "reflectivity",
+            "surface_normal",
+            (0, 0, 1),
+            does_not_raise(),
+            id="s2d2-reflectivity",
+        ),
+    ],
+)
+def test_reference_mode_solves_without_spurious_violation(
+    geometry, mode_name, ref_attr, ref_val, context
+):
+    """Reference-constraint modes solve end-to-end with verification active.
+
+    The post-solve verification loop now evaluates the reference-constraint
+    residual; a correct solver must produce solutions that pass it without
+    raising :class:`ConstraintViolation`.
+    """
+    with context:
+        g = _setup_cubic(geometry, **{ref_attr: ref_val})
+        g.mode_name = mode_name
+        assert g.mode.is_implemented(g)
+        sols = g.forward(1, 0, 1)
+        assert len(sols) > 0
+        # Re-running the verification loop explicitly must not raise.
+        _validate_solutions(sols, g.mode, g)
+        # Every returned solution satisfies the reference constraint.
+        rc = next(c for c in g.mode.constraints if c.category == "reference")
+        if rc.name != "psi" and rc.has_reference_vector(g):
+            for sol in sols:
+                assert abs(rc.evaluate(sol, g)) < 1e-6
+
+
+# ---------------------------------------------------------------------------
+# The verification loop rejects a bad reference solution.
+# ---------------------------------------------------------------------------
+
+
+def test_verification_rejects_bad_reference_solution():
+    """A solution violating a reference constraint raises ConstraintViolation."""
+    g = _setup_cubic("psic", surface_normal=(0, 0, 1))
+    g.mode_name = "fixed_incidence_vertical"
+    good = g.forward(1, 0, 1)[0]
+
+    # A mode demanding incidence=45° is not met by the incidence=0 solution.
+    bad_mode = ahd.ConstraintSet(
+        [ReferenceConstraint("incidence", 45.0)],
+        extras={"n_hat": ahd.REQUIRED, "incidence": None},
+    )
+    with pytest.raises(ConstraintViolation):
+        _validate_solutions([good], bad_mode, g)
+
+
+# ---------------------------------------------------------------------------
+# psi is skipped by the verification loop (enforced upstream as a filter).
+# ---------------------------------------------------------------------------
+
+
+def test_verification_skips_reference_without_vector():
+    """A reference constraint whose vector is unset is skipped, not raised."""
+    g = _setup_cubic("psic")  # no surface_normal set
+    assert g.surface_normal is None
+    mode = ahd.ConstraintSet(
+        [ReferenceConstraint("incidence", 0.0)],
+        extras={"n_hat": ahd.REQUIRED, "incidence": None},
+    )
+    # Any motor-angle dict works; the guard short-circuits before evaluate().
+    angles = {s.name: s.angle for s in g._stages.values()}  # noqa: SLF001
+    _validate_solutions([angles], mode, g)  # must not raise
+
+
+def test_verification_swallows_reference_evaluate_error(monkeypatch):
+    """If evaluate() raises despite a set vector, the loop skips that solution."""
+    g = _setup_cubic("psic", surface_normal=(0, 0, 1))
+    mode = ahd.ConstraintSet(
+        [ReferenceConstraint("incidence", 0.0)],
+        extras={"n_hat": ahd.REQUIRED, "incidence": None},
+    )
+    angles = {s.name: s.angle for s in g._stages.values()}  # noqa: SLF001
+
+    def _boom(*_args, **_kwargs):
+        raise ValueError("synthetic incidence failure for issue #304 coverage")
+
+    monkeypatch.setattr("ad_hoc_diffractometer.reference.incidence_angle", _boom)
+    _validate_solutions([angles], mode, g)  # must not raise
+
+
+def test_psi_constraint_skipped_by_verification_loop():
+    """psi residuals (subject to ±360° wrap) do not trip the verification loop."""
+    g = _setup_cubic("fourcv", azimuth=(0, 0, 1))
+    import numpy as np
+
+    from ad_hoc_diffractometer.forward import _compute_natural_psi
+
+    Q_phi = g.sample.UB @ np.array([1.0, 1.0, 1.0])
+    natural = _compute_natural_psi(g, Q_phi)
+    assert natural is not None
+
+    g.modes["fixed_psi"] = ahd.ConstraintSet(
+        [ReferenceConstraint("psi", natural)],
+        computed=g.modes["fixed_psi"].computed,
+        extras={"n_hat": ahd.REQUIRED, "psi": None},
+    )
+    g.mode_name = "fixed_psi"
+    sols = g.forward(1, 1, 1)
+    assert len(sols) > 0
+    # Must not raise even though a naive psi residual could be ~360°.
+    _validate_solutions(sols, g.mode, g)