fix(core): correct Crucible-confirmed QS + spectral physics defects

src/magnetics/_slcontour/fit.py

@@ -88,7 +88,7 @@ def form_basis_function(
         return (
             (np.exp(1j * m * np.deg2rad(y2)) - np.exp(1j * m * np.deg2rad(y1)))
             * (np.exp(1j * n * np.deg2rad(x2)) - np.exp(1j * n * np.deg2rad(x1)))
-        ) / (np.deg2rad(dx * dy) * n * m)
+        ) / (np.deg2rad(dx) * np.deg2rad(dy) * (1j * n) * (1j * m))
 
     # ── gaussian-point ────────────────────────────────────────────────────────
     if fit_basis == "gaussian-point":
@@ -352,7 +352,9 @@ def _printv(*a):
     w_inv = 1.0 / w_a
     w_inv[~valid] = 0.0
     w_inv = np.hstack((w_inv, [0.0] * max(Vh_a.shape[0] - w_a.shape[0], 0)))
-    fit_sigma2 = np.sum((Vh_a.T * w_inv) ** 2, axis=0)
+    # diag(cov) = diag(V W^-2 V^H) = Σ_k V[j,k]^2 w_inv[k]^2 — sum over the singular
+    # index k (axis=1 of Vh_a.T, shape [n_cols, k]), NOT over the coefficient index.
+    fit_sigma2 = np.sum((Vh_a.T * w_inv) ** 2, axis=1)
     fit_sigmas = np.sqrt(fit_sigma2)  # (n_columns,)
 
     # ── reform per-mode complex (sinusoidal) or real (Gaussian) coefficients ──

src/magnetics/core/qs_bridge.py

@@ -55,14 +55,6 @@ def _mode_label(n: int | float, m: int | float) -> str:
     return f"n={n}" if m == 0 else f"m/n={m}/{n}"
 
 
-def _quality_status(K: float) -> str:
-    if K > 20:
-        return "error"
-    if K > 10:
-        return "warn"
-    return "ok"
-
-
 def _reconstruct_grid(
     fit_ds, phi_grid: np.ndarray, theta_grid: np.ndarray, t_idx: int
 ) -> np.ndarray:
@@ -257,8 +249,12 @@ def fit_to_fit_quality_node(fit_ds) -> dict:
     return contracts.metrics(
         title="fit quality",
         fields=[
-            {"label": "K (raw)", "value": f"{K:.1f}", "status": _quality_status(K)},
-            {"label": "K (eff)", "value": f"{eff_cn:.1f}", "status": _quality_status(eff_cn)},
+            {"label": "K (raw)", "value": f"{K:.1f}", "status": contracts.quality_for_k(K)},
+            {
+                "label": "K (eff)",
+                "value": f"{eff_cn:.1f}",
+                "status": contracts.quality_for_k(eff_cn),
+            },
             {"label": "K cutoff", "value": f"{fit_cond:.0f}"},
             {"label": "χ² (mean)", "value": f"{mean_chi2:.3f}"},
             {"label": "channels", "value": n_ch},

src/magnetics/core/quasistationary.py

@@ -170,7 +170,10 @@ def fit(
 
     # Per-coeff uncertainty from SVD pseudo-inverse
     w_inv = np.where(valid, 1.0 / np.where(w_a != 0, w_a, 1.0), 0.0)
-    fit_sigmas = np.sqrt(np.sum((Vh_a.T * w_inv) ** 2, axis=0))  # [n_cols]
+    # diag(cov) = diag(V W^-2 V^H) = Σ_k V[j,k]^2 w_inv[k]^2 — sum over the singular
+    # index k (axis=1 of Vh_a.T, shape [n_cols, k]). axis=1 also yields length n_cols
+    # for an underdetermined fit (k = n_ch < n_cols), so the reform loop can't overrun.
+    fit_sigmas = np.sqrt(np.sum((Vh_a.T * w_inv) ** 2, axis=1))  # [n_cols]
 
     # ── reform complex coefficients (one complex number per mode) ─────────────
     coeffs_c: list[np.ndarray] = []
@@ -250,7 +253,9 @@ def reconstruct_grid(
     z = np.zeros((len(theta_grid), len(phi_grid)))
     for i, (n, m) in enumerate(zip(result.ns, result.ms)):
         c = result.coeffs[i, t_idx]
-        # outer product: [n_theta, n_phi]
+        # The fit's forward model is A@x = x_r·cos ψ + x_i·sin ψ = Re(conj(c)·e^{iψ})
+        # with c = x_r + i·x_i and ψ = mθ + nφ, so the reconstruction must use conj(c)
+        # (matching qs_bridge._reconstruct_grid); using c mirrors the map toroidally.
         basis = np.exp(1j * m * theta_rad)[:, None] * np.exp(1j * n * phi_rad)[None, :]
-        z += (c * basis).real
+        z += (np.conj(c) * basis).real
     return z

src/magnetics/core/spectral.py

@@ -800,9 +800,10 @@ def fit_toroidal_mode(
     best_n, best_R, best_c = candidates[best_j], float(rmag[best_j]), float(inter[best_j])
 
     # Uncertainty (eigspec-style): propagate per-probe phase σ from the cross-spectral
-    # statistics into (a) the intercept's 1σ via inverse-variance combination and
-    # (b) a posterior over candidate n from each hypothesis's weighted χ². Falls back
-    # to unit σ when phase_error is absent so older callers still get sane numbers.
+    # statistics into (a) the intercept's 1σ by propagating σ through the SAME weights w
+    # that produced the intercept, and (b) a posterior over candidate n from each
+    # hypothesis's weighted χ². Falls back to unit σ when phase_error is absent so older
+    # callers still get sane numbers.
     sigma = np.asarray(
         mode_result.phase_error
         if mode_result.phase_error is not None
@@ -825,7 +826,10 @@ def _wrap180(x):
     post = np.exp(-(chi2 - chi2.min()) / 2.0)
     post /= post.sum()
     n_confidence = float(post[best_j])
-    phase_sigma = float(np.sqrt(1.0 / np.sum(1.0 / sigma**2)))
+    # 1σ of the ACTUALLY-returned intercept — the w-weighted circular mean, not the
+    # inverse-variance-optimal one: Var(c) = Σ(w_k·σ_k)² / (Σ w_k)². (Σw > 0: w is
+    # replaced by ones above if every weight was zero.)
+    phase_sigma = float(np.sqrt(np.sum((w * sigma) ** 2)) / w.sum())
 
     return ToroidalFitResult(
         kind="toroidal_fit",

tests/test_crucible_regressions.py

@@ -0,0 +1,110 @@
+"""Regression tests for the 2026-07-01 Crucible review fixes.
+
+Each test names one confirmed defect and fails against the pre-fix behavior:
+  P3  SVD per-coefficient σ summed the wrong axis (wrong QS error bars)
+  P4  toroidal phase_sigma reported the BLUE variance, not the returned estimator's
+  P6  the 2D sinusoidal-integral basis was mis-normalized by -180/π for m≠0 modes
+  P7  reconstruct_grid used c instead of conj(c) (toroidal mirror); underdetermined
+      fits raised IndexError
+"""
+
+from __future__ import annotations
+
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+
+from magnetics._slcontour import fit as slfit
+from magnetics.core import quasistationary
+from magnetics.core.spectral import fit_toroidal_mode
+
+
+# ── P6: 2D sinusoidal-integral basis normalization ──────────────────────────
+def test_slcontour_2d_sinusoidal_integral_basis_is_separable_product():
+    """The (n≠0, m≠0) integral basis must equal the product of the two 1D branches.
+    Pre-fix it divided by deg2rad(dx*dy)*n*m (π/180 once, no i²), off by -180/π."""
+    x1 = np.array([10.0, 100.0])
+    x2 = x1 + 20.0
+    y1 = np.array([5.0, 60.0])
+    y2 = y1 + 15.0
+    n, m = 2, 3
+    fmn_2d = slfit.form_basis_function(n, m, x1, x2, y1, y2, fit_basis="sinusoidal-integral")
+    fmn_n = slfit.form_basis_function(n, 0, x1, x2, y1, y2, fit_basis="sinusoidal-integral")
+    fmn_m = slfit.form_basis_function(0, m, x1, x2, y1, y2, fit_basis="sinusoidal-integral")
+    assert np.allclose(fmn_2d, fmn_n * fmn_m)
+
+
+# ── P3: SVD per-coefficient uncertainty axis ────────────────────────────────
+def test_svd_covariance_diagonal_sums_the_singular_axis():
+    """diag((AᵀA)⁻¹) = Σ_k V[j,k]² / w_k² sums over the singular index (axis=1 of
+    Vh.T). The old axis=0 collapsed to 1/w_k² by singular position — not the
+    covariance diagonal. Guards the reduction in fit.py:355 and quasistationary.py:173."""
+    rng = np.random.default_rng(0)
+    a = rng.standard_normal((10, 4))  # non-orthogonal ⇒ condition number ≠ 1
+    _, w, vh = np.linalg.svd(a, full_matrices=False)
+    w_inv = 1.0 / w
+    ref = np.diag(np.linalg.inv(a.T @ a))
+    assert np.allclose(np.sum((vh.T * w_inv) ** 2, axis=1), ref)
+    assert not np.allclose(np.sum((vh.T * w_inv) ** 2, axis=0), ref)
+
+
+# ── P7: underdetermined fit must not IndexError (resolved by the axis fix) ───
+def test_quasistationary_fit_underdetermined_sizes_sigmas_without_indexerror():
+    phi = np.linspace(0.0, 300.0, 6)
+    theta = np.linspace(0.0, 150.0, 6)
+    signal = np.cos(np.deg2rad(phi))[:, None] + 0.3 * np.sin(np.deg2rad(theta))[:, None]
+    res = quasistationary.fit(
+        np.array([0.0]),
+        signal,
+        phi,
+        phi,
+        theta,
+        theta,
+        ns=(1, 2, 3),
+        ms=(0, 1, 2),  # 9 modes → ~18 basis columns > 6 sensors
+        fit_basis="sinusoidal-point",
+    )
+    assert res.sigmas.shape[0] == res.ns.shape[0]  # one σ per mode; overran before the fix
+    assert np.all(np.isfinite(np.abs(res.sigmas)))
+
+
+# ── P7: reconstruct_grid must not mirror a phase-shifted mode ────────────────
+def test_reconstruct_grid_roundtrips_phase_shifted_mode():
+    """cos(φ − 40°) fit as n=1, m=0 then reconstructed must reproduce cos(φ − 40°).
+    Pre-fix (c instead of conj(c)) it produced cos(φ + 40°) — a toroidal mirror."""
+    phi = np.linspace(0.0, 315.0, 8)
+    theta0 = np.zeros_like(phi)
+    signal = np.cos(np.deg2rad(phi - 40.0))[:, None]  # [n_ch, 1]
+    res = quasistationary.fit(
+        np.array([0.0]),
+        signal,
+        phi,
+        phi,
+        theta0,
+        theta0,
+        ns=(1,),
+        ms=(0,),
+        fit_basis="sinusoidal-point",
+    )
+    grid = np.linspace(0.0, 350.0, 36)
+    z = quasistationary.reconstruct_grid(res, grid, np.array([0.0]), 0)  # [1, n_phi]
+    assert np.allclose(z[0], np.cos(np.deg2rad(grid - 40.0)), atol=1e-6)
+
+
+# ── P4: toroidal phase_sigma reflects the actually-returned (weighted) estimator ─
+def test_toroidal_phase_sigma_matches_the_weighted_estimator_not_the_blue():
+    mode = SimpleNamespace(
+        frequency=3000.0,
+        toroidal_angle=np.array([0.0, 120.0, 240.0]),
+        phase=np.array([0.0, 240.0, 120.0]),
+        amplitude=np.array([100.0, 1.0, 1.0]),  # one loud, noisy probe dominates the weight
+        coherence=np.array([1.0, 1.0, 1.0]),
+        phase_error=np.array([30.0, 3.0, 3.0]),
+    )
+    fit = fit_toroidal_mode(mode)  # weights="amplitude"
+    w, sigma = mode.amplitude, mode.phase_error
+    expected = float(np.sqrt(np.sum((w * sigma) ** 2)) / w.sum())
+    blue = float(np.sqrt(1.0 / np.sum(1.0 / sigma**2)))  # the old, over-optimistic value
+    assert fit.phase_sigma == pytest.approx(expected)
+    assert fit.phase_sigma > 5 * blue

tests/test_qs_bridge.py

@@ -45,6 +45,16 @@ def test_dropping_a_required_variable_raises_not_silently_misserves(fit_ds):
         qs_bridge.fit_to_qs_fit_node(broken)
 
 
+def test_fit_quality_statuses_use_the_contract_vocabulary(fit_ds):
+    """The traffic-light status must be one of the GUI's Quality values
+    (good/warn/bad = contracts.quality_for_k), not the old 'ok'/'error' strings
+    which NodeView's QCOLOR can't color."""
+    node = qs_bridge.fit_to_fit_quality_node(fit_ds)
+    statuses = [f["status"] for f in node["fields"] if "status" in f]
+    assert statuses  # the K (raw)/K (eff) rows carry a status
+    assert all(s in {"good", "warn", "bad"} for s in statuses)
+
+
 def test_amplitude_sigma_is_finite(fit_ds):
     node = qs_bridge.fit_to_amplitude_node(fit_ds)
     sigma = np.asarray(node["meta"]["sigma"], dtype=float)