From b47c85789cb53699e3c46aafbb79b46c5ff687b8 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 16 Mar 2026 19:38:02 +0000
Subject: [PATCH 1/5] Initial plan


From a1ca1bcc5091443865222eeaa834c881a545d74b Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 16 Mar 2026 20:07:55 +0000
Subject: [PATCH 2/5] Fix calibration to use fixed source Stokes state

Co-authored-by: benjym <3380296+benjym@users.noreply.github.com>
---
 docs/user/configuration.md    |   7 +-
 photoelastimetry/calibrate.py | 117 ++++++++++++----------------------
 tests/test_calibrate.py       |  96 +++++++++++++++++++---------
 3 files changed, 113 insertions(+), 107 deletions(-)

diff --git a/docs/user/configuration.md b/docs/user/configuration.md
index 7e5f598..80ae97f 100644
--- a/docs/user/configuration.md
+++ b/docs/user/configuration.md
@@ -142,10 +142,13 @@ Common optional keys:
 
 - `dark_frame_file` + `blank_frame_file` (both or neither)
 - `fit.max_points`, `fit.seed`, `fit.loss`, `fit.f_scale`, `fit.max_nfev`
-- `fit.initial_C`, `fit.initial_S_i_hat`
-- `fit.s3_identifiability_threshold`, `fit.prior_weight`, `fit.c_relative_bounds`
+- `fit.initial_C`, `fit.S_i_hat` (or legacy `fit.initial_S_i_hat`)
+- `fit.prior_weight`, `fit.c_relative_bounds`
 - `output_profile`, `output_report`, `output_diagnostics`
 
+`fit.S_i_hat` is treated as a fixed source state during calibration; the
+optimiser only fits `C`.
+
 Load step input notes:
 
 - Each `load_steps[i].image_file` may be either:
diff --git a/photoelastimetry/calibrate.py b/photoelastimetry/calibrate.py
index 96f0f21..dd550f6 100644
--- a/photoelastimetry/calibrate.py
+++ b/photoelastimetry/calibrate.py
@@ -1,9 +1,9 @@
 """
 Calibration workflows for photoelastimetry experiments.
 
-This module fits per-wavelength stress-optic coefficients (C), incoming
-polarisation state (S_i_hat), and optional detector blank correction from a
-multi-load calibration sequence.
+This module fits per-wavelength stress-optic coefficients (C) for a fixed
+incoming polarisation state (S_i_hat), and optional detector blank correction
+from a multi-load calibration sequence.
 
 Supported calibration methods:
 - ``brazilian_disk``: diametrically-loaded disk with analytical stress field.
@@ -357,7 +357,7 @@ def validate_calibration_config(config):
     if n_no_load < 1:
         warnings.warn(
             "Calibration usually benefits from at least one no-load step (load ≈ 0); "
-            "continuing without it will use a default initial S_i_hat guess unless one is provided in fit.initial_S_i_hat.",
+            "continuing without it will use a default fixed S_i_hat unless one is provided in fit.S_i_hat or fit.initial_S_i_hat.",
             stacklevel=2,
         )
     if n_loaded < 3:
@@ -373,7 +373,6 @@ def validate_calibration_config(config):
     fit_cfg.setdefault("f_scale", 0.05)
     fit_cfg.setdefault("max_nfev", 300)
     fit_cfg.setdefault("seed", 0)
-    fit_cfg.setdefault("s3_identifiability_threshold", 0.02)
     fit_cfg.setdefault("prior_weight", 0.0)
 
     output_profile = config.get("output_profile", "calibration_profile.json5")
@@ -740,7 +739,7 @@ def _prepare_sampling_points(roi_mask, max_points, seed):
 
 
 def _initial_s_i_hat_from_noload(measured_noload):
-    """Estimate initial incoming Stokes state from no-load measurements."""
+    """Estimate a fallback incoming Stokes state from no-load measurements."""
     s1 = float(np.median(measured_noload[..., 0]))
     s2 = float(np.median(measured_noload[..., 1]))
     magnitude_sq = s1**2 + s2**2
@@ -851,7 +850,7 @@ def _build_dataset(config):
     if len(no_load_measurements) == 0:
         warnings.warn(
             "No no-load measurements were found after processing load steps; "
-            "using a default initial S_i_hat guess of [1, 0, 0].",
+            "using a default fixed S_i_hat of [1, 0, 0].",
             stacklevel=2,
         )
         initial_s_i_hat = _default_initial_s_i_hat()
@@ -909,7 +908,18 @@ def _decode_params(params, n_channels, fixed_s3=None):
     return C, s_i_hat
 
 
-def calibration_residuals(params, dataset, fixed_s3=None):
+def _resolve_fixed_s_i_hat(fit_config, dataset):
+    """Resolve the fixed incoming Stokes state used during calibration."""
+    source = fit_config.get("S_i_hat", fit_config.get("initial_S_i_hat", dataset["initial_s_i_hat"]))
+    s_i_hat = np.asarray(source, dtype=float)
+    if s_i_hat.size == 2:
+        s_i_hat = np.append(s_i_hat, 0.0)
+    if s_i_hat.size != 3:
+        raise ValueError(f"fit.S_i_hat must have length 2 or 3, got {s_i_hat.size}.")
+    return _normalise_s_i_hat(s_i_hat)
+
+
+def calibration_residuals(params, dataset, fixed_s3=None, fixed_s_i_hat=None):
     """
     Compute calibration residuals for least-squares fitting.
 
@@ -921,6 +931,9 @@ def calibration_residuals(params, dataset, fixed_s3=None):
         Dataset dictionary from `_build_dataset`.
     fixed_s3 : float, optional
         If provided, S3 is fixed and only S1/S2 are optimised.
+    fixed_s_i_hat : array-like, optional
+        If provided, use this full fixed source state instead of decoding
+        Stokes components from `params`.
 
     Returns
     -------
@@ -929,7 +942,11 @@ def calibration_residuals(params, dataset, fixed_s3=None):
     """
     wavelengths = dataset["wavelengths"]
     n_channels = wavelengths.size
-    C, s_i_hat = _decode_params(params, n_channels, fixed_s3=fixed_s3)
+    C = np.asarray(params[:n_channels], dtype=float)
+    if fixed_s_i_hat is None:
+        _, s_i_hat = _decode_params(params, n_channels, fixed_s3=fixed_s3)
+    else:
+        s_i_hat = _normalise_s_i_hat(fixed_s_i_hat)
 
     residual_chunks = []
     for measured, sigma_xx, sigma_yy, sigma_xy in zip(
@@ -959,7 +976,7 @@ def calibration_residuals(params, dataset, fixed_s3=None):
 
 def fit_calibration_parameters(dataset, fit_config):
     """
-    Fit C and S_i_hat from calibration dataset.
+    Fit C from calibration dataset with fixed S_i_hat.
 
     Parameters
     ----------
@@ -971,7 +988,7 @@ def fit_calibration_parameters(dataset, fit_config):
     Returns
     -------
     dict
-        Fit results containing estimated C, S_i_hat, metrics, and optimizer state.
+        Fit results containing estimated C, fixed S_i_hat, metrics, and optimizer state.
     """
     wavelengths = dataset["wavelengths"]
     n_channels = wavelengths.size
@@ -979,13 +996,7 @@ def fit_calibration_parameters(dataset, fit_config):
     init_c = _as_float_array(
         fit_config.get("initial_C", np.full(n_channels, 3e-9)), expected_length=n_channels, name="initial_C"
     )
-    init_s = np.asarray(fit_config.get("initial_S_i_hat", dataset["initial_s_i_hat"]), dtype=float)
-    if init_s.size == 2:
-        init_s = np.append(init_s, 0.0)
-    if init_s.size != 3:
-        raise ValueError(f"fit.initial_S_i_hat must have length 2 or 3, got {init_s.size}.")
-
-    init_s = _normalise_s_i_hat(init_s)
+    fixed_s_i_hat = _resolve_fixed_s_i_hat(fit_config, dataset)
 
     c_relative_bounds = fit_config.get("c_relative_bounds")
     if c_relative_bounds is not None:
@@ -1001,70 +1012,26 @@ def fit_calibration_parameters(dataset, fit_config):
         c_lower = np.full(n_channels, 1e-15)
         c_upper = np.full(n_channels, 1e-4)
 
-    x0_full = np.concatenate([init_c, init_s], axis=0)
-    lb_full = np.concatenate([c_lower, np.full(3, -1.0)])
-    ub_full = np.concatenate([c_upper, np.full(3, 1.0)])
+    x0 = init_c.copy()
     prior_weight = float(fit_config.get("prior_weight", 0.0))
 
-    def residual_full(params):
-        residual = calibration_residuals(params, dataset, fixed_s3=None)
+    def residual_c_only(params):
+        residual = calibration_residuals(params, dataset, fixed_s_i_hat=fixed_s_i_hat)
         if prior_weight > 0:
-            prior = prior_weight * (params - x0_full)
+            prior = prior_weight * (params - x0)
             return np.concatenate([residual, prior], axis=0)
         return residual
 
-    result_full = least_squares(
-        residual_full,
-        x0_full,
-        bounds=(lb_full, ub_full),
+    result = least_squares(
+        residual_c_only,
+        x0,
+        bounds=(c_lower, c_upper),
         loss=fit_config["loss"],
         f_scale=float(fit_config["f_scale"]),
         max_nfev=int(fit_config["max_nfev"]),
     )
-
-    col_norms = (
-        np.linalg.norm(result_full.jac, axis=0) if result_full.jac.size > 0 else np.zeros_like(x0_full)
-    )
-    reference_norm = np.median(col_norms[:-1]) if col_norms.size > 1 else 0.0
-    s3_ratio = float(col_norms[-1] / max(reference_norm, 1e-12)) if col_norms.size > 0 else 0.0
-
-    fallback = False
-    fallback_reason = None
-    threshold = float(fit_config.get("s3_identifiability_threshold", 0.02))
-
-    if (not result_full.success) or (s3_ratio < threshold):
-        fallback = True
-        fallback_reason = "s3_identifiability" if s3_ratio < threshold else "full_fit_failed"
-
-        if fallback_reason == "s3_identifiability":
-            # If S3 is unidentifiable, restart fixed-S3 fit from the configured
-            # initial guess to avoid inheriting unstable full-model updates.
-            x0_fix = np.concatenate([init_c, init_s[:2]], axis=0)
-        else:
-            c_guess, s_guess = _decode_params(result_full.x, n_channels, fixed_s3=None)
-            x0_fix = np.concatenate([c_guess, s_guess[:2]], axis=0)
-        lb_fix = np.concatenate([c_lower, np.full(2, -1.0)])
-        ub_fix = np.concatenate([c_upper, np.full(2, 1.0)])
-
-        def residual_fixed(params):
-            residual = calibration_residuals(params, dataset, fixed_s3=0.0)
-            if prior_weight > 0:
-                prior = prior_weight * (params - x0_fix)
-                return np.concatenate([residual, prior], axis=0)
-            return residual
-
-        result = least_squares(
-            residual_fixed,
-            x0_fix,
-            bounds=(lb_fix, ub_fix),
-            loss=fit_config["loss"],
-            f_scale=float(fit_config["f_scale"]),
-            max_nfev=int(fit_config["max_nfev"]),
-        )
-        C_fit, S_fit = _decode_params(result.x, n_channels, fixed_s3=0.0)
-    else:
-        result = result_full
-        C_fit, S_fit = _decode_params(result.x, n_channels, fixed_s3=None)
+    C_fit = np.asarray(result.x[:n_channels], dtype=float)
+    S_fit = fixed_s_i_hat
 
     residual = result.fun
     fit_metrics = {
@@ -1077,9 +1044,9 @@ def residual_fixed(params):
         "n_residuals": int(residual.size),
         "n_samples": int(dataset["sample_y"].size),
         "n_load_steps": int(dataset["loads"].size),
-        "fallback_used": bool(fallback),
-        "fallback_reason": fallback_reason,
-        "s3_identifiability_ratio": float(s3_ratio),
+        "fallback_used": False,
+        "fallback_reason": None,
+        "s3_identifiability_ratio": None,
     }
 
     return {
diff --git a/tests/test_calibrate.py b/tests/test_calibrate.py
index fc17091..5c27514 100644
--- a/tests/test_calibrate.py
+++ b/tests/test_calibrate.py
@@ -9,12 +9,14 @@
 from photoelastimetry.image import simulate_four_step_polarimetry
 
 
-def _make_synthetic_calibration_case(tmp_path, noisy=False, inconsistent_shapes=False):
+def _make_synthetic_calibration_case(
+    tmp_path, noise_scale=0.0, inconsistent_shapes=False, source_state=None, load_scale=1.0
+):
     height, width = 34, 34
     thickness = 0.01
     wavelengths = np.array([650e-9, 550e-9, 450e-9], dtype=float)
     c_true = np.array([2.2e-9, 2.5e-9, 2.8e-9], dtype=float)
-    s_true = np.array([0.82, 0.18, 0.0], dtype=float)
+    s_true = np.array([0.82, 0.18, 0.0] if source_state is None else source_state, dtype=float)
 
     radius_m = 0.008
     pixels_per_meter = 1800.0
@@ -51,10 +53,10 @@ def _make_synthetic_calibration_case(tmp_path, noisy=False, inconsistent_shapes=
         (offset + 1.0 / scale)[np.newaxis, np.newaxis, :, :], (height, width, 3, 4)
     ).copy()
 
-    if noisy:
+    if noise_scale > 0.0:
         rng = np.random.default_rng(42)
-        dark += rng.normal(scale=1e-6, size=dark.shape)
-        blank += rng.normal(scale=1e-6, size=blank.shape)
+        dark += rng.normal(scale=0.02 * noise_scale, size=dark.shape)
+        blank += rng.normal(scale=0.02 * noise_scale, size=blank.shape)
 
     dark_file = tmp_path / "dark.npy"
     blank_file = tmp_path / "blank.npy"
@@ -88,15 +90,15 @@ def _make_synthetic_calibration_case(tmp_path, noisy=False, inconsistent_shapes=
         corrected[~disk_mask, :, :] = 1.0
 
         raw = corrected / scale[np.newaxis, np.newaxis, :, :] + offset[np.newaxis, np.newaxis, :, :]
-        if noisy:
-            raw += rng.normal(scale=5e-5, size=raw.shape)
+        if noise_scale > 0.0:
+            raw += rng.normal(scale=noise_scale, size=raw.shape)
 
         if inconsistent_shapes and idx == len(loads) - 1:
             raw = raw[:-1, :, :, :]
 
         image_file = tmp_path / f"load_{idx}.npy"
         io.save_image(str(image_file), raw.astype(np.float32))
-        load_steps.append({"load": float(load), "image_file": str(image_file)})
+        load_steps.append({"load": float(load_scale * load), "image_file": str(image_file)})
 
     config = {
         "method": "brazilian_disk",
@@ -119,10 +121,10 @@ def _make_synthetic_calibration_case(tmp_path, noisy=False, inconsistent_shapes=
             "loss": "soft_l1",
             "f_scale": 0.05,
             "max_nfev": 250,
-            "prior_weight": 10000.0 if noisy else 0.0,
-            "c_relative_bounds": [0.5, 2.0] if noisy else None,
+            "prior_weight": 10000.0 if noise_scale > 0.0 else 0.0,
+            "c_relative_bounds": [0.5, 2.0] if noise_scale > 0.0 else None,
             "initial_C": [2.0e-9, 2.0e-9, 2.0e-9],
-            "initial_S_i_hat": [0.8, 0.1, 0.0],
+            "S_i_hat": s_true.tolist(),
         },
         "output_profile": str(tmp_path / "calibration_profile.json5"),
         "output_report": str(tmp_path / "calibration_report.md"),
@@ -211,13 +213,13 @@ def test_build_dataset_accepts_raw_recording_directory_load_steps(tmp_path):
     assert dataset["measured_steps"][0].shape[-2:] == (3, 2)
 
 
-def _make_synthetic_coupon_case(tmp_path):
+def _make_synthetic_coupon_case(tmp_path, noise_scale=0.0, source_state=None, load_scale=1.0):
     height, width = 30, 42
     thickness = 0.01
     coupon_width_m = 0.012
     wavelengths = np.array([650e-9, 550e-9, 450e-9], dtype=float)
     c_true = np.array([2.4e-9, 2.7e-9, 3.1e-9], dtype=float)
-    s_true = np.array([0.88, 0.12, 0.0], dtype=float)
+    s_true = np.array([0.88, 0.12, 0.0] if source_state is None else source_state, dtype=float)
 
     # [x0, x1, y0, y1]
     gauge_roi = np.array([10, 32, 8, 22], dtype=int)
@@ -251,10 +253,17 @@ def _make_synthetic_coupon_case(tmp_path):
 
     dark_file = tmp_path / "coupon_dark.npy"
     blank_file = tmp_path / "coupon_blank.npy"
+
+    if noise_scale > 0.0:
+        rng = np.random.default_rng(314)
+        dark += rng.normal(scale=0.02 * noise_scale, size=dark.shape)
+        blank += rng.normal(scale=0.02 * noise_scale, size=blank.shape)
+
     io.save_image(str(dark_file), dark.astype(np.float32))
     io.save_image(str(blank_file), blank.astype(np.float32))
 
     load_steps = []
+    rng = np.random.default_rng(2718)
     for idx, load in enumerate(loads):
         sigma = load / (thickness * coupon_width_m)
         sigma_xx = np.full((height, width), sigma, dtype=float)
@@ -273,10 +282,12 @@ def _make_synthetic_coupon_case(tmp_path):
 
         corrected[~roi_mask, :, :] = 1.0
         raw = corrected / scale[np.newaxis, np.newaxis, :, :] + offset[np.newaxis, np.newaxis, :, :]
+        if noise_scale > 0.0:
+            raw += rng.normal(scale=noise_scale, size=raw.shape)
 
         image_file = tmp_path / f"coupon_load_{idx}.npy"
         io.save_image(str(image_file), raw.astype(np.float32))
-        load_steps.append({"load": float(load), "image_file": str(image_file)})
+        load_steps.append({"load": float(load_scale * load), "image_file": str(image_file)})
 
     config = {
         "method": "coupon_test",
@@ -300,9 +311,9 @@ def _make_synthetic_coupon_case(tmp_path):
             "f_scale": 0.05,
             "max_nfev": 200,
             "initial_C": [2.2e-9, 2.6e-9, 3.0e-9],
-            "initial_S_i_hat": [0.85, 0.1, 0.0],
-            "c_relative_bounds": [0.8, 1.2],
-            "prior_weight": 8000.0,
+            "S_i_hat": s_true.tolist(),
+            "c_relative_bounds": [0.8, 1.2] if noise_scale > 0.0 else None,
+            "prior_weight": 8000.0 if noise_scale > 0.0 else 0.0,
         },
         "output_profile": str(tmp_path / "coupon_calibration_profile.json5"),
         "output_report": str(tmp_path / "coupon_calibration_report.md"),
@@ -312,28 +323,39 @@ def _make_synthetic_coupon_case(tmp_path):
 
 
 def test_calibration_residuals_near_zero_for_noiseless_data(tmp_path):
-    config, c_true, s_true = _make_synthetic_calibration_case(tmp_path, noisy=False)
+    config, c_true, s_true = _make_synthetic_calibration_case(tmp_path)
     validated = calibrate.validate_calibration_config(config)
     dataset = calibrate._build_dataset(validated)
 
-    params = np.concatenate([c_true, s_true], axis=0)
-    residual = calibrate.calibration_residuals(params, dataset)
+    residual = calibrate.calibration_residuals(c_true, dataset, fixed_s_i_hat=s_true)
 
     assert residual.ndim == 1
     assert np.sqrt(np.mean(residual**2)) < 1e-6
 
 
-def test_fit_recovers_c_and_s_i_hat_on_noisy_data(tmp_path):
-    config, c_true, s_true = _make_synthetic_calibration_case(tmp_path, noisy=True)
+@pytest.mark.parametrize(
+    ("noise_scale", "source_state"),
+    [
+        (0.0, [0.82, 0.18, 0.0]),
+        (5e-5, [0.82, 0.18, 0.0]),
+        (0.0, [0.2, 0.1, np.sqrt(0.95)]),
+        (5e-5, [0.2, 0.1, np.sqrt(0.95)]),
+    ],
+)
+def test_fit_recovers_c_with_fixed_s_i_hat_from_synthetic_disk_data(tmp_path, noise_scale, source_state):
+    config, c_true, s_true = _make_synthetic_calibration_case(
+        tmp_path, noise_scale=noise_scale, source_state=source_state
+    )
     result = calibrate.run_calibration(config)
     profile = result["profile"]
 
     c_fit = np.asarray(profile["C"], dtype=float)
     s_fit = np.asarray(profile["S_i_hat"], dtype=float)
 
-    assert np.allclose(c_fit, c_true, rtol=0.2, atol=0.0)
-    assert np.allclose(s_fit[:2], s_true[:2], atol=0.12)
+    assert np.allclose(c_fit, c_true, rtol=0.25, atol=0.0)
+    assert np.allclose(s_fit, calibrate._normalise_s_i_hat(s_true), atol=1e-12)
     assert profile["fit_metrics"]["success"]
+    assert profile["fit_metrics"]["fallback_used"] is False
 
 
 def test_blank_correction_coefficients_and_application():
@@ -368,7 +390,7 @@ def test_blank_correction_coefficients_and_application():
 
 
 def test_validation_errors_for_bad_configuration(tmp_path):
-    config, _, _ = _make_synthetic_calibration_case(tmp_path, noisy=False)
+    config, _, _ = _make_synthetic_calibration_case(tmp_path)
 
     no_noload = dict(config)
     no_noload["load_steps"] = [dict(step, load=abs(step["load"]) + 1.0) for step in config["load_steps"]]
@@ -395,22 +417,23 @@ def test_validation_errors_for_bad_configuration(tmp_path):
 
 
 def test_build_dataset_without_noload_uses_default_initial_guess(tmp_path):
-    config, _, _ = _make_synthetic_calibration_case(tmp_path, noisy=False)
+    config, _, _ = _make_synthetic_calibration_case(tmp_path)
     config["load_steps"] = [
         dict(step, load=float(index + 1)) for index, step in enumerate(config["load_steps"][:3])
     ]
+    config["fit"].pop("S_i_hat")
 
     with pytest.warns(UserWarning, match="no-load"):
         validated = calibrate.validate_calibration_config(config)
 
-    with pytest.warns(UserWarning, match="default initial S_i_hat"):
+    with pytest.warns(UserWarning, match="default fixed S_i_hat"):
         dataset = calibrate._build_dataset(validated)
 
     assert np.allclose(dataset["initial_s_i_hat"], np.array([1.0, 0.0, 0.0]))
 
 
 def test_validation_error_for_inconsistent_image_shapes(tmp_path):
-    config, _, _ = _make_synthetic_calibration_case(tmp_path, noisy=False, inconsistent_shapes=True)
+    config, _, _ = _make_synthetic_calibration_case(tmp_path, inconsistent_shapes=True)
     validated = calibrate.validate_calibration_config(config)
 
     with pytest.raises(ValueError, match="share shape"):
@@ -418,7 +441,7 @@ def test_validation_error_for_inconsistent_image_shapes(tmp_path):
 
 
 def test_end_to_end_writes_profile_report_and_diagnostics(tmp_path):
-    config, _, _ = _make_synthetic_calibration_case(tmp_path, noisy=False)
+    config, _, _ = _make_synthetic_calibration_case(tmp_path)
     result = calibrate.run_calibration(config)
 
     assert Path(result["profile_file"]).exists()
@@ -459,7 +482,7 @@ def test_coupon_test_end_to_end_fit_and_profile(tmp_path):
 
     assert profile["method"] == "coupon_test"
     assert np.allclose(np.asarray(profile["C"], dtype=float), c_true, rtol=0.2, atol=0.0)
-    assert np.allclose(np.asarray(profile["S_i_hat"], dtype=float)[:2], s_true[:2], atol=0.12)
+    assert np.allclose(np.asarray(profile["S_i_hat"], dtype=float), s_true, atol=1e-12)
     assert profile["fit_metrics"]["success"]
 
     loaded = calibrate.load_calibration_profile(result["profile_file"])
@@ -474,3 +497,16 @@ def test_coupon_validation_requires_coupon_geometry_fields(tmp_path):
 
     with pytest.raises(ValueError, match="coupon_width_m"):
         calibrate.validate_calibration_config(bad)
+
+
+@pytest.mark.parametrize("noise_scale", [0.0, 5e-5])
+def test_coupon_calibration_recovers_c_with_fixed_circular_source(tmp_path, noise_scale):
+    config, c_true, s_true = _make_synthetic_coupon_case(
+        tmp_path, noise_scale=noise_scale, source_state=[0.0, 0.0, 1.0], load_scale=0.5
+    )
+    result = calibrate.run_calibration(config)
+    profile = result["profile"]
+
+    assert np.allclose(np.asarray(profile["C"], dtype=float), c_true, rtol=0.2, atol=0.0)
+    assert np.allclose(np.asarray(profile["S_i_hat"], dtype=float), s_true, atol=1e-12)
+    assert profile["fit_metrics"]["success"]

From f60b833d16d47c92a1970a982fd91ddc28f8dab6 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 16 Mar 2026 20:10:34 +0000
Subject: [PATCH 3/5] Add fixed source-state calibration coverage

Co-authored-by: benjym <3380296+benjym@users.noreply.github.com>
---
 photoelastimetry/calibrate.py | 7 ++++++-
 tests/test_calibrate.py       | 5 ++---
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/photoelastimetry/calibrate.py b/photoelastimetry/calibrate.py
index dd550f6..a1258d8 100644
--- a/photoelastimetry/calibrate.py
+++ b/photoelastimetry/calibrate.py
@@ -910,7 +910,12 @@ def _decode_params(params, n_channels, fixed_s3=None):
 
 def _resolve_fixed_s_i_hat(fit_config, dataset):
     """Resolve the fixed incoming Stokes state used during calibration."""
-    source = fit_config.get("S_i_hat", fit_config.get("initial_S_i_hat", dataset["initial_s_i_hat"]))
+    if "S_i_hat" in fit_config:
+        source = fit_config["S_i_hat"]
+    elif "initial_S_i_hat" in fit_config:
+        source = fit_config["initial_S_i_hat"]
+    else:
+        source = dataset["initial_s_i_hat"]
     s_i_hat = np.asarray(source, dtype=float)
     if s_i_hat.size == 2:
         s_i_hat = np.append(s_i_hat, 0.0)
diff --git a/tests/test_calibrate.py b/tests/test_calibrate.py
index 5c27514..b2cd0db 100644
--- a/tests/test_calibrate.py
+++ b/tests/test_calibrate.py
@@ -254,8 +254,8 @@ def _make_synthetic_coupon_case(tmp_path, noise_scale=0.0, source_state=None, lo
     dark_file = tmp_path / "coupon_dark.npy"
     blank_file = tmp_path / "coupon_blank.npy"
 
+    rng = np.random.default_rng(314)
     if noise_scale > 0.0:
-        rng = np.random.default_rng(314)
         dark += rng.normal(scale=0.02 * noise_scale, size=dark.shape)
         blank += rng.normal(scale=0.02 * noise_scale, size=blank.shape)
 
@@ -263,7 +263,6 @@ def _make_synthetic_coupon_case(tmp_path, noise_scale=0.0, source_state=None, lo
     io.save_image(str(blank_file), blank.astype(np.float32))
 
     load_steps = []
-    rng = np.random.default_rng(2718)
     for idx, load in enumerate(loads):
         sigma = load / (thickness * coupon_width_m)
         sigma_xx = np.full((height, width), sigma, dtype=float)
@@ -507,6 +506,6 @@ def test_coupon_calibration_recovers_c_with_fixed_circular_source(tmp_path, nois
     result = calibrate.run_calibration(config)
     profile = result["profile"]
 
-    assert np.allclose(np.asarray(profile["C"], dtype=float), c_true, rtol=0.2, atol=0.0)
+    assert np.allclose(np.asarray(profile["C"], dtype=float), c_true, rtol=0.25, atol=0.0)
     assert np.allclose(np.asarray(profile["S_i_hat"], dtype=float), s_true, atol=1e-12)
     assert profile["fit_metrics"]["success"]

From f2ba7280ab48b8368e884bc75b5b6e16a6a1a2ff Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 16 Mar 2026 20:17:48 +0000
Subject: [PATCH 4/5] Require explicit calibration source state

Co-authored-by: benjym <3380296+benjym@users.noreply.github.com>
---
 docs/user/configuration.md    |   9 ++-
 photoelastimetry/calibrate.py | 104 +++++++---------------------------
 tests/test_calibrate.py       |  23 ++++++--
 3 files changed, 45 insertions(+), 91 deletions(-)

diff --git a/docs/user/configuration.md b/docs/user/configuration.md
index 80ae97f..21f3944 100644
--- a/docs/user/configuration.md
+++ b/docs/user/configuration.md
@@ -142,12 +142,15 @@ Common optional keys:
 
 - `dark_frame_file` + `blank_frame_file` (both or neither)
 - `fit.max_points`, `fit.seed`, `fit.loss`, `fit.f_scale`, `fit.max_nfev`
-- `fit.initial_C`, `fit.S_i_hat` (or legacy `fit.initial_S_i_hat`)
-- `fit.prior_weight`, `fit.c_relative_bounds`
+- `fit.initial_C`, `fit.prior_weight`, `fit.c_relative_bounds`
 - `output_profile`, `output_report`, `output_diagnostics`
 
+Required fit key:
+
+- `fit.S_i_hat` as an explicit 3-component fixed source state `[S1_hat, S2_hat, S3_hat]`
+
 `fit.S_i_hat` is treated as a fixed source state during calibration; the
-optimiser only fits `C`.
+optimiser only fits `C`. Calibration does not provide a default `S_i_hat`.
 
 Load step input notes:
 
diff --git a/photoelastimetry/calibrate.py b/photoelastimetry/calibrate.py
index a1258d8..379cf69 100644
--- a/photoelastimetry/calibrate.py
+++ b/photoelastimetry/calibrate.py
@@ -1,9 +1,9 @@
 """
 Calibration workflows for photoelastimetry experiments.
 
-This module fits per-wavelength stress-optic coefficients (C) for a fixed
-incoming polarisation state (S_i_hat), and optional detector blank correction
-from a multi-load calibration sequence.
+This module fits per-wavelength stress-optic coefficients (C) for an explicit
+fixed incoming polarisation state (S_i_hat), and optional detector blank
+correction from a multi-load calibration sequence.
 
 Supported calibration methods:
 - ``brazilian_disk``: diametrically-loaded disk with analytical stress field.
@@ -357,7 +357,7 @@ def validate_calibration_config(config):
     if n_no_load < 1:
         warnings.warn(
             "Calibration usually benefits from at least one no-load step (load ≈ 0); "
-            "continuing without it will use a default fixed S_i_hat unless one is provided in fit.S_i_hat or fit.initial_S_i_hat.",
+            "continuing without it may reduce robustness.",
             stacklevel=2,
         )
     if n_loaded < 3:
@@ -368,12 +368,19 @@ def validate_calibration_config(config):
         )
 
     fit_cfg = dict(config.get("fit", {}))
+    if "S_i_hat" not in fit_cfg:
+        raise ValueError("Calibration config must include fit.S_i_hat.")
+    fit_s_i_hat = np.asarray(fit_cfg["S_i_hat"], dtype=float)
+    if fit_s_i_hat.shape != (3,):
+        raise ValueError(f"fit.S_i_hat must have shape (3,), got {fit_s_i_hat.shape}.")
+
     fit_cfg.setdefault("max_points", 6000)
     fit_cfg.setdefault("loss", "soft_l1")
     fit_cfg.setdefault("f_scale", 0.05)
     fit_cfg.setdefault("max_nfev", 300)
     fit_cfg.setdefault("seed", 0)
     fit_cfg.setdefault("prior_weight", 0.0)
+    fit_cfg["S_i_hat"] = fit_s_i_hat.tolist()
 
     output_profile = config.get("output_profile", "calibration_profile.json5")
     output_report = config.get("output_report", "calibration_report.md")
@@ -738,27 +745,6 @@ def _prepare_sampling_points(roi_mask, max_points, seed):
     return y_idx, x_idx
 
 
-def _initial_s_i_hat_from_noload(measured_noload):
-    """Estimate a fallback incoming Stokes state from no-load measurements."""
-    s1 = float(np.median(measured_noload[..., 0]))
-    s2 = float(np.median(measured_noload[..., 1]))
-    magnitude_sq = s1**2 + s2**2
-    s3 = float(np.sqrt(max(0.0, 1.0 - magnitude_sq)))
-    s_i = np.array([s1, s2, s3], dtype=float)
-
-    norm = np.linalg.norm(s_i)
-    if norm < 1e-12:
-        return np.array([1.0, 0.0, 0.0], dtype=float)
-    if norm > 1.0:
-        s_i = s_i / norm
-    return s_i
-
-
-def _default_initial_s_i_hat():
-    """Return a conservative default incoming Stokes state guess."""
-    return np.array([1.0, 0.0, 0.0], dtype=float)
-
-
 def _build_dataset(config):
     """Load calibration images and construct regression dataset."""
     steps = config["load_steps"]
@@ -813,8 +799,6 @@ def _build_dataset(config):
     sigma_yy_steps = []
     sigma_xy_steps = []
 
-    no_load_tolerance = config["load_zero_tolerance"]
-    no_load_measurements = []
     diagnostic_step_index = int(np.argmax(np.abs(loads)))
     diagnostic_image = corrected_images[diagnostic_step_index]
 
@@ -844,20 +828,6 @@ def _build_dataset(config):
             sigma_yy_steps.append(sigma_yy)
             sigma_xy_steps.append(sigma_xy)
 
-        if abs(load) <= no_load_tolerance:
-            no_load_measurements.append(measured)
-
-    if len(no_load_measurements) == 0:
-        warnings.warn(
-            "No no-load measurements were found after processing load steps; "
-            "using a default fixed S_i_hat of [1, 0, 0].",
-            stacklevel=2,
-        )
-        initial_s_i_hat = _default_initial_s_i_hat()
-    else:
-        noload = np.concatenate(no_load_measurements, axis=0)
-        initial_s_i_hat = _initial_s_i_hat_from_noload(noload)
-
     dataset = {
         "method": config["method"],
         "wavelengths": config["wavelengths"],
@@ -875,7 +845,6 @@ def _build_dataset(config):
         # Backward-compatible key retained for previous diagnostics consumers.
         "disk_mask": model_mask,
         "blank_correction": blank_correction,
-        "initial_s_i_hat": initial_s_i_hat,
         "X": X,
         "Y": Y,
         "geometry": config["geometry"],
@@ -897,48 +866,26 @@ def _normalise_s_i_hat(s_i_hat):
     return s_i
 
 
-def _decode_params(params, n_channels, fixed_s3=None):
-    """Decode parameter vector into C and S_i_hat."""
-    C = np.asarray(params[:n_channels], dtype=float)
-    if fixed_s3 is None:
-        s_i_hat = np.asarray(params[n_channels : n_channels + 3], dtype=float)
-    else:
-        s_i_hat = np.array([params[n_channels], params[n_channels + 1], fixed_s3], dtype=float)
-    s_i_hat = _normalise_s_i_hat(s_i_hat)
-    return C, s_i_hat
-
-
-def _resolve_fixed_s_i_hat(fit_config, dataset):
+def _resolve_fixed_s_i_hat(fit_config):
     """Resolve the fixed incoming Stokes state used during calibration."""
-    if "S_i_hat" in fit_config:
-        source = fit_config["S_i_hat"]
-    elif "initial_S_i_hat" in fit_config:
-        source = fit_config["initial_S_i_hat"]
-    else:
-        source = dataset["initial_s_i_hat"]
-    s_i_hat = np.asarray(source, dtype=float)
-    if s_i_hat.size == 2:
-        s_i_hat = np.append(s_i_hat, 0.0)
-    if s_i_hat.size != 3:
-        raise ValueError(f"fit.S_i_hat must have length 2 or 3, got {s_i_hat.size}.")
+    s_i_hat = np.asarray(fit_config["S_i_hat"], dtype=float)
+    if s_i_hat.shape != (3,):
+        raise ValueError(f"fit.S_i_hat must have shape (3,), got {s_i_hat.shape}.")
     return _normalise_s_i_hat(s_i_hat)
 
 
-def calibration_residuals(params, dataset, fixed_s3=None, fixed_s_i_hat=None):
+def calibration_residuals(params, dataset, fixed_s_i_hat):
     """
     Compute calibration residuals for least-squares fitting.
 
     Parameters
     ----------
     params : ndarray
-        Parameter vector containing C values and S_i_hat components.
+        Parameter vector containing C values.
     dataset : dict
         Dataset dictionary from `_build_dataset`.
-    fixed_s3 : float, optional
-        If provided, S3 is fixed and only S1/S2 are optimised.
-    fixed_s_i_hat : array-like, optional
-        If provided, use this full fixed source state instead of decoding
-        Stokes components from `params`.
+    fixed_s_i_hat : array-like
+        Fixed source state used for the forward model.
 
     Returns
     -------
@@ -946,12 +893,8 @@ def calibration_residuals(params, dataset, fixed_s3=None, fixed_s_i_hat=None):
         1D residual vector.
     """
     wavelengths = dataset["wavelengths"]
-    n_channels = wavelengths.size
-    C = np.asarray(params[:n_channels], dtype=float)
-    if fixed_s_i_hat is None:
-        _, s_i_hat = _decode_params(params, n_channels, fixed_s3=fixed_s3)
-    else:
-        s_i_hat = _normalise_s_i_hat(fixed_s_i_hat)
+    C = np.asarray(params[: wavelengths.size], dtype=float)
+    s_i_hat = _normalise_s_i_hat(fixed_s_i_hat)
 
     residual_chunks = []
     for measured, sigma_xx, sigma_yy, sigma_xy in zip(
@@ -1001,7 +944,7 @@ def fit_calibration_parameters(dataset, fit_config):
     init_c = _as_float_array(
         fit_config.get("initial_C", np.full(n_channels, 3e-9)), expected_length=n_channels, name="initial_C"
     )
-    fixed_s_i_hat = _resolve_fixed_s_i_hat(fit_config, dataset)
+    fixed_s_i_hat = _resolve_fixed_s_i_hat(fit_config)
 
     c_relative_bounds = fit_config.get("c_relative_bounds")
     if c_relative_bounds is not None:
@@ -1049,9 +992,6 @@ def residual_c_only(params):
         "n_residuals": int(residual.size),
         "n_samples": int(dataset["sample_y"].size),
         "n_load_steps": int(dataset["loads"].size),
-        "fallback_used": False,
-        "fallback_reason": None,
-        "s3_identifiability_ratio": None,
     }
 
     return {
diff --git a/tests/test_calibrate.py b/tests/test_calibrate.py
index b2cd0db..1b993bf 100644
--- a/tests/test_calibrate.py
+++ b/tests/test_calibrate.py
@@ -203,6 +203,7 @@ def test_build_dataset_accepts_raw_recording_directory_load_steps(tmp_path):
             {"image_file": str(recording_dir), "load": 2.0},
             {"image_file": str(recording_dir), "load": 3.0},
         ],
+        "fit": {"S_i_hat": [1.0, 0.0, 0.0]},
     }
 
     validated = calibrate.validate_calibration_config(config)
@@ -354,7 +355,6 @@ def test_fit_recovers_c_with_fixed_s_i_hat_from_synthetic_disk_data(tmp_path, no
     assert np.allclose(c_fit, c_true, rtol=0.25, atol=0.0)
     assert np.allclose(s_fit, calibrate._normalise_s_i_hat(s_true), atol=1e-12)
     assert profile["fit_metrics"]["success"]
-    assert profile["fit_metrics"]["fallback_used"] is False
 
 
 def test_blank_correction_coefficients_and_application():
@@ -414,21 +414,32 @@ def test_validation_errors_for_bad_configuration(tmp_path):
     with pytest.raises(ValueError, match="at least one load_step"):
         calibrate.validate_calibration_config(empty_steps)
 
+    missing_s_i_hat = dict(config)
+    missing_s_i_hat["fit"] = dict(config["fit"])
+    missing_s_i_hat["fit"].pop("S_i_hat")
+    with pytest.raises(ValueError, match="fit.S_i_hat"):
+        calibrate.validate_calibration_config(missing_s_i_hat)
 
-def test_build_dataset_without_noload_uses_default_initial_guess(tmp_path):
+    legacy_s_i_hat = dict(config)
+    legacy_s_i_hat["fit"] = dict(config["fit"])
+    legacy_s_i_hat["fit"].pop("S_i_hat")
+    legacy_s_i_hat["fit"]["initial_S_i_hat"] = [0.8, 0.1, 0.0]
+    with pytest.raises(ValueError, match="fit.S_i_hat"):
+        calibrate.validate_calibration_config(legacy_s_i_hat)
+
+
+def test_validation_without_noload_still_accepts_explicit_s_i_hat(tmp_path):
     config, _, _ = _make_synthetic_calibration_case(tmp_path)
     config["load_steps"] = [
         dict(step, load=float(index + 1)) for index, step in enumerate(config["load_steps"][:3])
     ]
-    config["fit"].pop("S_i_hat")
 
     with pytest.warns(UserWarning, match="no-load"):
         validated = calibrate.validate_calibration_config(config)
 
-    with pytest.warns(UserWarning, match="default fixed S_i_hat"):
-        dataset = calibrate._build_dataset(validated)
+    dataset = calibrate._build_dataset(validated)
 
-    assert np.allclose(dataset["initial_s_i_hat"], np.array([1.0, 0.0, 0.0]))
+    assert dataset["loads"].shape == (3,)
 
 
 def test_validation_error_for_inconsistent_image_shapes(tmp_path):

From c1aa426ca4812d2862c808a9f328590aef198c3c Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 16 Mar 2026 20:22:47 +0000
Subject: [PATCH 5/5] Remove calibration source-state fallback paths

Co-authored-by: benjym <3380296+benjym@users.noreply.github.com>
---
 photoelastimetry/calibrate.py | 8 ++++++--
 tests/test_calibrate.py       | 5 ++++-
 2 files changed, 10 insertions(+), 3 deletions(-)

diff --git a/photoelastimetry/calibrate.py b/photoelastimetry/calibrate.py
index 379cf69..115fdfb 100644
--- a/photoelastimetry/calibrate.py
+++ b/photoelastimetry/calibrate.py
@@ -368,8 +368,12 @@ def validate_calibration_config(config):
         )
 
     fit_cfg = dict(config.get("fit", {}))
+    if "initial_S_i_hat" in fit_cfg:
+        raise ValueError("fit.initial_S_i_hat is no longer supported; use fit.S_i_hat.")
     if "S_i_hat" not in fit_cfg:
-        raise ValueError("Calibration config must include fit.S_i_hat.")
+        raise ValueError(
+            "Calibration config must include fit.S_i_hat as a 3-component array " "[S1_hat, S2_hat, S3_hat]."
+        )
     fit_s_i_hat = np.asarray(fit_cfg["S_i_hat"], dtype=float)
     if fit_s_i_hat.shape != (3,):
         raise ValueError(f"fit.S_i_hat must have shape (3,), got {fit_s_i_hat.shape}.")
@@ -881,7 +885,7 @@ def calibration_residuals(params, dataset, fixed_s_i_hat):
     Parameters
     ----------
     params : ndarray
-        Parameter vector containing C values.
+        1D parameter vector containing one C value per wavelength channel.
     dataset : dict
         Dataset dictionary from `_build_dataset`.
     fixed_s_i_hat : array-like
diff --git a/tests/test_calibrate.py b/tests/test_calibrate.py
index 1b993bf..76a354f 100644
--- a/tests/test_calibrate.py
+++ b/tests/test_calibrate.py
@@ -424,7 +424,7 @@ def test_validation_errors_for_bad_configuration(tmp_path):
     legacy_s_i_hat["fit"] = dict(config["fit"])
     legacy_s_i_hat["fit"].pop("S_i_hat")
     legacy_s_i_hat["fit"]["initial_S_i_hat"] = [0.8, 0.1, 0.0]
-    with pytest.raises(ValueError, match="fit.S_i_hat"):
+    with pytest.raises(ValueError, match="initial_S_i_hat"):
         calibrate.validate_calibration_config(legacy_s_i_hat)
 
 
@@ -440,6 +440,9 @@ def test_validation_without_noload_still_accepts_explicit_s_i_hat(tmp_path):
     dataset = calibrate._build_dataset(validated)
 
     assert dataset["loads"].shape == (3,)
+    assert np.allclose(
+        np.asarray(validated["fit"]["S_i_hat"], dtype=float), np.asarray(config["fit"]["S_i_hat"])
+    )
 
 
 def test_validation_error_for_inconsistent_image_shapes(tmp_path):