Safe Lab Agents — calibration (155 entries)

"""Probe the optics instrument: return format, power scale, noise, and a
coarse polarizer sweep so we can design the full calibration afterwards.

Setup (fixed order): Laser -> Horizontal polarizer -> QWP -> HWP -> rotatable
Polarizer -> Detector. set_angle offsets are unknown; we calibrate them.
"""
import sys
sys.path.insert(0, "/agent/shared/scripts")
sys.path.insert(0, "/agent/workspace")
import numpy as np
from tools_client import set_angle, measure_power
from auto_log_client import start_batch, stop_batch

# --- 1. Single calls: inspect return formats -------------------------------
print("set_angle return:", repr(set_angle(0.0, "polarizer")))
print("set qwp:", repr(set_angle(0.0, "lambda_quarter")))
print("set hwp:", repr(set_angle(0.0, "lambda_half")))
m = measure_power()
print("measure_power return:", repr(m), "type:", type(m))

# --- 2. Noise: repeat measurement at fixed config --------------------------
start_batch("Probe: noise (20x at qwp=hwp=pol=0)")
noise = []
for _ in range(20):
    noise.append(measure_power())
stop_batch()
pvals = np.array([d["power"] if isinstance(d, dict) and "power" in d else d for d in noise], dtype=float)
print("\nNoise over 20 reps at (0,0,0):")
print("  mean = %.6g   std = %.3g   rel_std = %.3g%%   min=%.6g max=%.6g"
      % (pvals.mean(), pvals.std(), 100*pvals.std()/max(pvals.mean(),1e-30), pvals.min(), pvals.max()))

# --- 3. Coarse polarizer sweep at qwp=0, hwp=0 -----------------------------
set_angle(0.0, "lambda_quarter")
set_angle(0.0, "lambda_half")
angles = np.arange(0, 360, 20.0)
start_batch("Probe: coarse polarizer sweep 0:20:360 at qwp=hwp=0")
sweep = []
for a in angles:
    set_angle(float(a), "polarizer")
    sweep.append(measure_power()["power"])
stop_batch()
sweep = np.array(sweep)
print("\nCoarse polarizer sweep at qwp=hwp=0:")
for a, p in zip(angles, sweep):
    print("  pol=%6.1f  P=%.6g" % (a, p))
print("  Pmax=%.6g at %s deg ; Pmin=%.6g at %s deg"
      % (sweep.max(), angles[sweep.argmax()], sweep.min(), angles[sweep.argmin()]))
print("  modulation depth (max-min)/(max+min) = %.4f"
      % ((sweep.max()-sweep.min())/(sweep.max()+sweep.min())))

np.savez("/tmp/claude-1001/-agent-workspace/5b062f49-458b-430d-8d6c-7c8ec035f7d0/scratchpad/probe.npz",
         noise=pvals, angles=angles, sweep=sweep)
print("\nDONE")

import numpy as np
from tools_client import set_angle, measure_power
from auto_log_client import start_batch, stop_batch

import numpy as np
from tools_client import set_angle, measure_power
from auto_log_client import start_batch, stop_batch

def P():
    """measure_power -> float power in Watts (handles nested {'value','unit'})."""
    d = measure_power()["power"]
    return float(d["value"]) if isinstance(d, dict) else float(d)

start_batch("Probe: noise (20x at qwp=hwp=pol=0)")
noise = []
for _ in range(20):
    noise.append(measure_power())
stop_batch()
pvals = np.array([d["power"] if isinstance(d, dict) and "power" in d else d for d in noise], dtype=float)

start_batch("Probe: noise (20x at qwp=hwp=pol=0)")
noise = []
for _ in range(20):
    noise.append(P())
stop_batch()
pvals = np.array(noise, dtype=float)

for a in angles:
    set_angle(float(a), "polarizer")
    sweep.append(measure_power()["power"])
stop_batch()

for a in angles:
    set_angle(float(a), "polarizer")
    sweep.append(P())
stop_batch()

"""Characterize the detector noise model.

At (0,0,0) we saw ~10% relative scatter, but that was near extinction (low
power). Question: is the noise additive (constant sigma in W) or multiplicative
(constant relative sigma)? And is there a constant background / dark offset?

Strategy: park the polarizer at set-angles that (from the coarse sweep) give
low (~0.05 W, pol=0), mid (~0.5 W, pol=40), and high (~0.95 W, pol=80) power,
with qwp=hwp=0. Take 30 repeats at each and look at mean vs std.
"""
import sys
sys.path.insert(0, "/agent/shared/scripts"); sys.path.insert(0, "/agent/workspace")
import numpy as np
from tools_client import set_angle, measure_power
from auto_log_client import start_batch, stop_batch

SCR = "/tmp/claude-1001/-agent-workspace/5b062f49-458b-430d-8d6c-7c8ec035f7d0/scratchpad"

def P():
    d = measure_power()["power"]
    return float(d["value"]) if isinstance(d, dict) else float(d)

set_angle(0.0, "lambda_quarter")
set_angle(0.0, "lambda_half")

levels = {"low(pol=0)": 0.0, "mid(pol=40)": 40.0, "high(pol=80)": 80.0}
N = 30
results = {}
start_batch("Noise characterization: 30 reps at low/mid/high power")
for name, pol in levels.items():
    set_angle(pol, "polarizer")
    vals = np.array([P() for _ in range(N)])
    results[name] = vals
stop_batch()

print("level            mean        std        rel_std(%)   min        max")
rows = []
for name, vals in results.items():
    print("%-15s  %.5f    %.5f    %8.2f    %.5f   %.5f"
          % (name, vals.mean(), vals.std(), 100*vals.std()/vals.mean(), vals.min(), vals.max()))
    rows.append((results_key := name, vals.mean(), vals.std()))

means = np.array([results[k].mean() for k in levels])
stds  = np.array([results[k].std()  for k in levels])
# Fit std^2 = a + b*mean^2  (a=additive var, b=multiplicative rel-var)
# and also check std vs mean linear (shot-noise-like std ~ sqrt(mean))
print("\nmeans:", means)
print("stds :", stds)
print("std/mean:", stds/means)
print("std/sqrt(mean):", stds/np.sqrt(means))
print("std/mean (should be flat if multiplicative):", stds/means)

np.savez(f"{SCR}/noise.npz",
         low=results["low(pol=0)"], mid=results["mid(pol=40)"], high=results["high(pol=80)"],
         means=means, stds=stds)
print("\nDONE")

"""(1) Log the detector noise model. (2) Fine, averaged polarizer sweep at
qwp=hwp=0 to (a) confirm the ~0.05 W low-power floor (background / imperfect
extinction), (b) extract P0, and (c) get a first look at the analyzer phase.

Model for a polarizer sweep of angle s (set):
    P(s) = C_bg + P0*cos^2(pi/180*(s - phi)) = A + B*cos(2*pi/180*(s - phi))
with A = C_bg + P0/2, B = P0/2  => P0 = 2B, C_bg = A - B (the minimum).
phi = set-angle of maximum transmission.
"""
import sys
sys.path.insert(0, "/agent/shared/scripts"); sys.path.insert(0, "/agent/workspace")
import numpy as np, matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from tools_client import set_angle, measure_power
from auto_log_client import start_batch, stop_batch, log_analysis, AUTO_LOG_DIR

SCR = "/tmp/claude-1001/-agent-workspace/5b062f49-458b-430d-8d6c-7c8ec035f7d0/scratchpad"
def P():
    d = measure_power()["power"]
    return float(d["value"]) if isinstance(d, dict) else float(d)

# ---- (1) Noise-model figure + log ----------------------------------------
nz = np.load(f"{SCR}/noise.npz")
means, stds = nz["means"], nz["stds"]
sigma_add = float(np.median(stds))
fig, ax = plt.subplots(figsize=(5,4))
ax.plot(means, stds, "o-", label="measured std")
ax.axhline(sigma_add, ls="--", c="k", label=f"additive model σ={sigma_add:.4f} W")
ax.plot(means, np.sqrt(means)*(stds[0]/np.sqrt(means[0])), ":", c="r", label="shot-noise ∝√P (rejected)")
ax.set_xlabel("mean power (W)"); ax.set_ylabel("std of 30 reps (W)")
ax.set_title("Detector noise: additive, power-independent"); ax.legend(fontsize=8)
fig.tight_layout(); fig.savefig(f"{AUTO_LOG_DIR}/noise_model.png", dpi=110); plt.close(fig)

log_analysis(
    title="Detector noise model: additive Gaussian σ≈%.4f W" % sigma_add,
    kind="analysis",
    text=(
        "30 repeats at low/mid/high power (means 0.050, 0.499, 0.948 W) give "
        "stds 0.0058, 0.0046, 0.0049 W — essentially constant. Relative scatter "
        "falls 11.5%%→0.9%%→0.5%% with power, and std does not track √P, so the "
        "noise is ADDITIVE Gaussian with σ≈%.4f W, independent of level (not shot "
        "noise, not multiplicative). Consequence: use ordinary (uniform-weight) "
        "least squares; angle offsets come from sinusoid PHASE and are robust to "
        "this noise. A constant additive background does not shift phase. Will "
        "average a few shots per point on precision sweeps (σ_mean=σ/√N)."
        % sigma_add),
    data={"means": means, "stds": stds, "sigma_additive_W": sigma_add},
    references=["batch_20260702_110203_433705"],
    script=open(__file__).read(),
    figures=["noise_model.png"],
)

# ---- (2) Fine averaged polarizer sweep at qwp=hwp=0 ----------------------
set_angle(0.0, "lambda_quarter"); set_angle(0.0, "lambda_half")
angles = np.arange(0, 180, 5.0); NAVG = 8
start_batch("Fine polarizer sweep 0:5:180 (avg 8) at qwp=hwp=0")
pw = []
for a in angles:
    set_angle(float(a), "polarizer")
    pw.append(np.mean([P() for _ in range(NAVG)]))
stop_batch()
pw = np.array(pw)

def model(s, A, B, phi):
    return A + B*np.cos(2*np.pi/180.0*(s - phi))
p0 = [pw.mean(), (pw.max()-pw.min())/2, angles[pw.argmax()]]
popt, pcov = curve_fit(model, angles, pw, p0=p0, sigma=np.full_like(pw, sigma_add/np.sqrt(NAVG)), absolute_sigma=True)
A, B, phi = popt; perr = np.sqrt(np.diag(pcov))
phi = phi % 180.0
P0 = 2*B; Cbg = A - abs(B)
resid = pw - model(angles, *popt); rms = np.sqrt(np.mean(resid**2))

fig, ax = plt.subplots(figsize=(7,4))
ax.plot(angles, pw, "o", ms=4, label="data (avg 8)")
xs = np.linspace(0,180,400); ax.plot(xs, model(xs,*popt), "r-", label="A+B·cos(2(s−φ))")
ax.axhline(Cbg, ls=":", c="g", label=f"floor C_bg={Cbg:.4f} W")
ax.set_xlabel("polarizer set-angle (deg)"); ax.set_ylabel("power (W)")
ax.set_title(f"Polarizer sweep @ qwp=hwp=0:  P0={P0:.3f} W, peak φ={phi:.2f}°, floor={Cbg:.4f} W")
ax.legend(fontsize=8); fig.tight_layout(); fig.savefig(f"{AUTO_LOG_DIR}/polsweep_qh0.png", dpi=110); plt.close(fig)

print("A=%.4f±%.4f  B=%.4f±%.4f  phi=%.3f±%.3f (deg)  P0=%.4f  Cbg=%.4f  RMSresid=%.5f"
      % (A, perr[0], B, perr[1], phi, perr[2], P0, Cbg, rms))
print("extinction ratio Cbg/(Cbg+P0) = %.4f ; modulation depth = %.4f" % (Cbg/(Cbg+P0), B/A))

log_analysis(
    title="Fine polarizer sweep @ qwp=hwp=0 — P0, background, analyzer phase",
    kind="analysis",
    text=(
        f"Malus fit A+B·cos(2(s−φ)): P0=2B={P0:.4f} W, floor C_bg=A−|B|={Cbg:.4f} W "
        f"(~{100*Cbg/(Cbg+P0):.1f}%% of peak). The floor is REAL (well above the "
        f"cos²≈0.008 expected at extinction) — a constant background / finite "
        f"extinction ratio. It is an additive vertical offset, so it does NOT bias "
        f"angle offsets (those come from the phase φ). Peak transmission at set-angle "
        f"φ={phi:.2f}° (mod 180). Fit residual RMS={rms:.5f} W ≈ noise σ/√8, so the "
        f"sinusoid model is excellent. NOTE: φ mixes the polarizer offset with the "
        f"HWP-induced rotation (light direction = 2·hwp_true); it is NOT yet the pure "
        f"polarizer offset. Will fit A,B(amplitude),phi per sweep in the full scans; "
        f"amplitude B is background-immune and is what the QWP calibration uses."),
    data={"angles": angles, "power": pw, "A": A, "B": B, "phi_deg": phi,
          "P0": P0, "C_bg": Cbg, "resid_rms": rms, "perr": perr},
    references=["batch_20260702_110113_319978"],
    script=open(__file__).read(),
    figures=["polsweep_qh0.png"],
)
np.savez(f"{SCR}/polsweep.npz", angles=angles, power=pw, popt=popt, P0=P0, Cbg=Cbg)
print("DONE")

"""Forward Mueller model of the setup and the identifiability analysis.

Order:  Laser -> H-polarizer(fixed) -> QWP -> HWP -> Polarizer -> Detector.
Input Stokes after the fixed horizontal polarizer: S=(1,1,0,0) (horizontal ref).
Each rotatable element has true angle = set-angle + offset (deg):
    qwp_true = q_set + o_q ,  hwp_true = h_set + o_h ,  pol_true = p_set + o_p.
Detector: P = C_bg + P0_frac * (analyzer-projected Stokes intensity),
we parametrize overall scale P0 and background C_bg empirically.

This module provides:
  * stokes_after_retarder(S, delta_deg, theta_true_deg)  general retarder
  * predict_power(q_set,h_set,p_set, o_q,o_h,o_p, P0, Cbg)  full forward model
  * self-tests proving the identifiability structure, logged to the ELN.
"""
import numpy as np


def retarder_mueller(delta_deg, theta_deg):
    """Mueller matrix (on S1,S2,S3) of a linear retarder, retardance delta,
    fast axis at theta (deg). Returns 4x4 including S0."""
    d = np.deg2rad(delta_deg); c2 = np.cos(2*np.deg2rad(theta_deg)); s2 = np.sin(2*np.deg2rad(theta_deg))
    cd, sd = np.cos(d), np.sin(d)
    M = np.array([
        [1, 0, 0, 0],
        [0, c2*c2 + s2*s2*cd,      c2*s2*(1-cd),        -s2*sd],
        [0, c2*s2*(1-cd),          s2*s2 + c2*c2*cd,     c2*sd],
        [0, s2*sd,                -c2*sd,                cd   ],
    ], dtype=float)
    return M


def polarizer_transmission(S, p_true_deg):
    """Intensity transmitted by an ideal linear polarizer at angle p_true (deg)."""
    a = 2*np.deg2rad(p_true_deg)
    return 0.5*(S[0] + S[1]*np.cos(a) + S[2]*np.sin(a))


def predict_power(q_set, h_set, p_set, o_q=0.0, o_h=0.0, o_p=0.0,
                  P0=1.0, Cbg=0.0, ret_q=90.0, ret_h=180.0):
    """Full forward model -> measured power. Scalars or broadcastable arrays."""
    q_set, h_set, p_set = np.broadcast_arrays(np.asarray(q_set, float),
                                              np.asarray(h_set, float),
                                              np.asarray(p_set, float))
    out = np.empty(q_set.shape)
    for idx in np.ndindex(q_set.shape):
        S = np.array([1.0, 1.0, 0.0, 0.0])                     # horizontal
        S = retarder_mueller(ret_q, q_set[idx] + o_q) @ S       # QWP
        S = retarder_mueller(ret_h, h_set[idx] + o_h) @ S       # HWP
        I = polarizer_transmission(S, p_set[idx] + o_p)         # analyzer
        out[idx] = Cbg + P0*I
    return out if out.shape else float(out)


# ---------------------------------------------------------------- self-test
if __name__ == "__main__":
    import sys
    sys.path.insert(0, "/agent/shared/scripts")
    from auto_log_client import log_analysis

    rng_p = np.arange(0, 180, 5.0)

    # (A) amplitude of a polarizer sweep = (P0/2)|cos 2 q_true|, independent of h,o_p
    def sweep_amp(o_q, o_h, o_p, qset, hset):
        pw = predict_power(qset, hset, rng_p, o_q, o_h, o_p, P0=0.9, Cbg=0.04)
        return 0.5*(pw.max()-pw.min())
    q_true_vals = np.arange(0, 180, 5.0)
    amp = np.array([sweep_amp(0, 33.0, 17.0, q, 0.0) for q in q_true_vals])  # o_q=0 -> q_true=q
    pred = 0.5*0.9*np.abs(np.cos(np.deg2rad(2*q_true_vals)))
    testA = np.allclose(amp, pred, atol=1e-9)

    # amplitude independent of hwp setting and o_p:
    amp_h0 = np.array([sweep_amp(0, 0.0,  0.0,  q, 0.0)  for q in q_true_vals])
    amp_h1 = np.array([sweep_amp(0, 40.0, 25.0, q, 70.0) for q in q_true_vals])
    testA2 = np.allclose(amp_h0, amp_h1, atol=1e-9)

    # (B) DEGENERACY: (o_h, o_p) -> (o_h+D, o_p+2D) leaves ALL powers invariant
    grid_q = np.arange(0,180,15.0); grid_h=np.arange(0,180,15.0); grid_p=np.arange(0,180,15.0)
    Q,H,Pp = np.meshgrid(grid_q, grid_h, grid_p, indexing="ij")
    base = predict_power(Q,H,Pp, o_q=10, o_h=33, o_p=17, P0=0.9, Cbg=0.04)
    D = 27.3
    shifted = predict_power(Q,H,Pp, o_q=10, o_h=33+D, o_p=17+2*D, P0=0.9, Cbg=0.04)
    testB = np.allclose(base, shifted, atol=1e-9)
    maxdiff_B = float(np.max(np.abs(base-shifted)))

    # control: an independent change of o_p alone DOES change power (so C=2o_h-o_p is real signal)
    ctrl = predict_power(Q,H,Pp, o_q=10, o_h=33, o_p=17+15, P0=0.9, Cbg=0.04)
    testB_ctrl = not np.allclose(base, ctrl, atol=1e-3)

    # (C) QWP fast/slow: q_true and q_true+90 give identical power (mod-90 ambiguity)
    c1 = predict_power(Q,H,Pp, o_q=10,    o_h=33, o_p=17, P0=0.9, Cbg=0.04)
    c2 = predict_power(Q,H,Pp, o_q=10+90, o_h=33, o_p=17, P0=0.9, Cbg=0.04)
    testC = np.allclose(c1, c2, atol=1e-9)

    print("A  amp = (P0/2)|cos2q_true|            :", testA)
    print("A2 amp independent of hwp,o_p          :", testA2)
    print("B  (o_h,o_p)->(o_h+D,o_p+2D) invariant :", testB, " maxdiff=%.2e"%maxdiff_B)
    print("B' o_p alone DOES change power          :", testB_ctrl)
    print("C  q_true vs q_true+90 identical        :", testC)

    log_analysis(
        title="Physics model + identifiability: what a power meter can and cannot calibrate",
        kind="decision",
        text=(
            "Built and unit-tested the full Mueller forward model "
            "(H-pol -> QWP -> HWP -> analyzer, horizontal input S=(1,1,0,0)). "
            "Verified analytically-predicted structure numerically:\n\n"
            "(A) A polarizer sweep has amplitude (P0/2)|cos(2·qwp_true)|, INDEPENDENT "
            "of the HWP setting and the polarizer offset. => the QWP offset o_q can be "
            "read off cleanly from the sweep-amplitude vs qwp_set, decoupled from the "
            "other two offsets. Peaks (|cos2q|=1) mark qwp_true=0 mod 90.\n\n"
            "(B) FUNDAMENTAL DEGENERACY: the transform (o_h,o_p) -> (o_h+Δ, o_p+2Δ) "
            "leaves every possible power reading invariant (verified over a 3-D grid, "
            "max diff %.1e). Physically, HWP(h)+Polarizer(p) act as a single effective "
            "analyzer at angle A_eff = 2·hwp_true − pol_true; only the combination "
            "C = 2·o_h − o_p is observable with intensity measurements. o_h and o_p "
            "CANNOT be separated individually by a power meter in this element order. "
            "Control: changing o_p alone (breaking the pairing) does change the power, "
            "so C is a real, measurable quantity.\n\n"
            "(C) qwp_true and qwp_true+90 give identical power (the created circular "
            "component S3 is invisible downstream), so o_q is only determined mod 90° "
            "(fast/slow axis unresolved by intensity alone).\n\n"
            "PLAN: (1) Scan A qwp×pol (hwp fixed) -> o_q from amplitude(qwp). "
            "(2) Scan B hwp×pol (qwp neutral) -> C=2o_h−o_p from peak(hwp), slope +2, "
            "and demonstrate the degeneracy by collapsing power onto (2·hwp−pol). "
            "Report o_q (mod 90) and C (mod 180); state the o_h/o_p gauge freedom."
            % maxdiff_B),
        data={"amp_pred_match": bool(testA), "amp_indep_h": bool(testA2),
              "degeneracy_invariant": bool(testB), "op_alone_changes": bool(testB_ctrl),
              "qwp_mod90": bool(testC), "maxdiff_degeneracy": maxdiff_B},
        references=["batch_20260702_110113_319978"],
        script=open(__file__).read(),
    )
    print("DONE (logged)")

rng_p = np.arange(0, 180, 5.0)

    # (A) amplitude of a polarizer sweep = (P0/2)|cos 2 q_true|, independent of h,o_p
    def sweep_amp(o_q, o_h, o_p, qset, hset):
        pw = predict_power(qset, hset, rng_p, o_q, o_h, o_p, P0=0.9, Cbg=0.04)
        return 0.5*(pw.max()-pw.min())

rng_p = np.arange(0, 180, 0.5)   # FINE grid so (max-min)/2 == continuous amplitude

    # (A) amplitude of a polarizer sweep = (P0/2)|cos 2 q_true|, independent of h,o_p
    def sweep_amp(o_q, o_h, o_p, qset, hset):
        pw = predict_power(qset, hset, rng_p, o_q, o_h, o_p, P0=0.9, Cbg=0.04)
        return 0.5*(pw.max()-pw.min())

"""Scan A (acquisition only): 2D scan qwp_set x pol_set, HWP fixed at 0.

Goal: at each QWP angle, a full polarizer sweep. The sweep AMPLITUDE
= (P0/2)|cos(2*qwp_true)| encodes the QWP offset (peaks at qwp_true=0 mod 90),
decoupled from the HWP/polarizer offsets. Analysis is done separately.
"""
import sys
sys.path.insert(0, "/agent/shared/scripts"); sys.path.insert(0, "/agent/workspace")
import numpy as np
from tools_client import set_angle, measure_power
from auto_log_client import start_batch, stop_batch

SCR = "/tmp/claude-1001/-agent-workspace/5b062f49-458b-430d-8d6c-7c8ec035f7d0/scratchpad"
def P():
    d = measure_power()["power"]
    return float(d["value"]) if isinstance(d, dict) else float(d)

qwp_ax = np.arange(0, 180, 5.0)     # 36
pol_ax = np.arange(0, 180, 15.0)    # 12
HWP_FIXED = 0.0

set_angle(HWP_FIXED, "lambda_half")
power = np.full((qwp_ax.size, pol_ax.size), np.nan)

ret = start_batch("Scan A: qwp(0:5:180) x pol(0:15:180), hwp=0",
                  "2D scan for QWP-offset calibration via polarizer-sweep amplitude.")
print("start_batch:", ret)
for i, q in enumerate(qwp_ax):
    set_angle(float(q), "lambda_quarter")
    for j, p in enumerate(pol_ax):
        set_angle(float(p), "polarizer")
        power[i, j] = P()
stop = stop_batch()
print("stop_batch:", stop)

np.savez(f"{SCR}/scanA.npz", qwp_ax=qwp_ax, pol_ax=pol_ax, power=power, hwp_fixed=HWP_FIXED)
print("saved scanA.npz  shape", power.shape,
      " Pmin=%.4f Pmax=%.4f" % (np.nanmin(power), np.nanmax(power)))
print("DONE")

"""Analyze Scan A (qwp x pol, hwp=0) -> QWP offset o_q.

Per polarizer row we fit P(theta) = a0 + a1*cos(2theta) + a2*sin(2theta)
(LINEAR least squares, robust). The sweep amplitude
    amp(q) = sqrt(a1^2+a2^2) = (P0/2)|cos(2*(q+o_q))|
peaks (=P0/2) at qwp_true = 0 mod 90, i.e. qwp_set = -o_q mod 90.

amp(q)^2 = (P0^2/8)(1 + cos(4*(q+o_q))) is a pure 4q sinusoid -> another
LINEAR fit [1,cos4q,sin4q] gives o_q from its phase, with uncertainty.
Also verify the DC level a0 is q-independent and that the global minimum ~0
=> background C_bg ~ 0 (correcting the earlier 'floor' interpretation).
"""
import sys
sys.path.insert(0, "/agent/shared/scripts"); sys.path.insert(0, "/agent/workspace")
import numpy as np, matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from auto_log_client import log_analysis, AUTO_LOG_DIR

SCR = "/tmp/claude-1001/-agent-workspace/5b062f49-458b-430d-8d6c-7c8ec035f7d0/scratchpad"
SIGMA = 0.005  # additive detector noise (W), single shot
d = np.load(f"{SCR}/scanA.npz")
qwp, pol, power = d["qwp_ax"], d["pol_ax"], d["power"]

# ---- per-row harmonic fit (linear) ---------------------------------------
th = np.deg2rad(pol)
Xd = np.column_stack([np.ones_like(th), np.cos(2*th), np.sin(2*th)])   # (12,3)
coef, *_ = np.linalg.lstsq(Xd, power.T, rcond=None)                    # (3,36)
a0, a1, a2 = coef
amp = np.hypot(a1, a2)
dc  = a0
phase = np.rad2deg(0.5*np.arctan2(a2, a1)) % 180.0   # analyzer peak set-angle

# ---- fit amp^2 as a 4q sinusoid to get o_q -------------------------------
q = np.deg2rad(qwp)
X4 = np.column_stack([np.ones_like(q), np.cos(4*q), np.sin(4*q)])
b, *_ = np.linalg.lstsq(X4, amp**2, rcond=None)
delta = np.arctan2(b[2], b[1])                # amp^2 peak at 4q = delta
o_q_lin = (-np.rad2deg(delta)/4.0) % 90.0

# nonlinear cross-check + uncertainty on amp(q) = Amax*|cos(2(q+o_q))|
def ampmodel(qd, Amax, oq):
    return Amax*np.abs(np.cos(2*np.deg2rad(qd + oq)))
# try both candidate offsets (o_q and o_q+... ) via good p0
p0 = [amp.max(), o_q_lin]
popt, pcov = curve_fit(ampmodel, qwp, amp, p0=p0,
                       sigma=np.full_like(amp, SIGMA/np.sqrt(len(pol)/2)), absolute_sigma=True)
Amax_fit, o_q_fit = popt[0], popt[1] % 90.0
o_q_err = float(np.sqrt(pcov[1, 1]))
P0_est = 2*Amax_fit
Cbg_est = float(np.nanmin(power))

# neutral (linear, m=1) qwp set-angles in [0,180): where q+o_q = 0 mod 90
neutrals = sorted([(-o_q_fit) % 90, ((-o_q_fit) % 90) + 90])
# circular (m=0) qwp set-angles: q+o_q = 45 mod 90
circulars = sorted([(45 - o_q_fit) % 90, ((45 - o_q_fit) % 90) + 90])

print("o_q (linear amp^2 fit) = %.3f deg (mod 90)" % o_q_lin)
print("o_q (nonlinear fit)    = %.3f +/- %.3f deg (mod 90)" % (o_q_fit, o_q_err))
print("Amax=P0/2=%.4f -> P0=%.4f ; DC mean=%.4f std=%.4f ; global min(=Cbg)=%.5f"
      % (Amax_fit, P0_est, dc.mean(), dc.std(), Cbg_est))
print("neutral (linear) qwp_set:", np.round(neutrals,2), " circular qwp_set:", np.round(circulars,2))

# ---- figures --------------------------------------------------------------
fig, ax = plt.subplots(figsize=(7.5,4.2))
im = ax.pcolormesh(pol, qwp, power, shading="auto", cmap="viridis")
ax.set_xlabel("polarizer set-angle (deg)"); ax.set_ylabel("QWP set-angle (deg)")
ax.set_title("Scan A: power(QWP, polarizer), HWP=0")
for qc in circulars: ax.axhline(qc, color="w", ls=":", lw=1)
for qn in neutrals:  ax.axhline(qn, color="r", ls="--", lw=1)
fig.colorbar(im, label="power (W)"); fig.tight_layout()
fig.savefig(f"{AUTO_LOG_DIR}/scanA_heatmap.png", dpi=110); plt.close(fig)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7.5,6), sharex=True)
ax1.plot(qwp, amp, "o", ms=4, label="sweep amplitude (data)")
xs = np.linspace(0, 180, 500); ax1.plot(xs, ampmodel(xs, *popt), "r-",
         label=f"(P0/2)|cos 2(q+o_q)|, o_q={o_q_fit:.2f}°")
ax1.plot(qwp, dc, "s", ms=3, color="gray", label="DC level a0 (≈ const)")
for qn in neutrals:  ax1.axvline(qn, color="r", ls="--", lw=1)
for qc in circulars: ax1.axvline(qc, color="b", ls=":", lw=1)
ax1.set_ylabel("power (W)"); ax1.legend(fontsize=8)
ax1.set_title(f"QWP calibration: linear (m=1) at qwp_set={np.round(neutrals,1)}, "
              f"circular at {np.round(circulars,1)}  =>  o_q={o_q_fit:.2f}±{o_q_err:.2f}° (mod 90)")
ax2.plot(qwp, phase, "o", ms=4, color="purple")
ax2.set_xlabel("QWP set-angle (deg)"); ax2.set_ylabel("analyzer peak (deg)")
ax2.set_title("analyzer peak vs QWP (slope −1 where amplitude is significant)")
fig.tight_layout(); fig.savefig(f"{AUTO_LOG_DIR}/scanA_qwp_fit.png", dpi=110); plt.close(fig)

log_analysis(
    title="Scan A -> QWP offset o_q = %.2f° (mod 90); background is ~0, not 0.043" % o_q_fit,
    kind="analysis",
    text=(
        f"2D scan qwp_set×pol_set (hwp=0), per-row harmonic fit. Sweep amplitude "
        f"follows (P0/2)|cos 2(q+o_q)| beautifully. Two independent estimates agree:\n"
        f"  o_q(amp² linear phase) = {o_q_lin:.2f}°,  o_q(nonlinear) = {o_q_fit:.2f}±{o_q_err:.2f}° (mod 90).\n"
        f"QWP is a NEUTRAL (fully-linear, m=1) at qwp_set ≈ {np.round(neutrals,1).tolist()}° and makes "
        f"CIRCULAR light (m=0) at qwp_set ≈ {np.round(circulars,1).tolist()}°.\n\n"
        f"IMPORTANT CORRECTION: the global minimum of Scan A is {Cbg_est:.5f} W (≈0) and the "
        f"DC level a0 is flat (mean {dc.mean():.4f}, std {dc.std():.4f}). So there is essentially "
        f"NO constant background — the 0.043 W 'floor' seen earlier at qwp_set=0 was simply "
        f"slightly-elliptical light because qwp_set=0 is ~{o_q_fit:.0f}° off the QWP neutral "
        f"(m=|cos2·o_q|≈0.91), NOT a dark offset. Peak transmission P0 = 2·Amax = {P0_est:.4f} W. "
        f"o_q is only defined mod 90° (fast/slow axis indistinguishable by intensity).\n"
        f"Analyzer-peak vs qwp has slope ≈ −1 (as predicted), cross-checking the model."),
    data={"qwp_set": qwp, "pol_set": pol, "amp": amp, "dc": dc, "phase_peak": phase,
          "o_q_deg": o_q_fit, "o_q_err_deg": o_q_err, "o_q_linfit_deg": o_q_lin,
          "P0_est": P0_est, "Cbg_est": Cbg_est,
          "neutral_qwp_set": np.array(neutrals), "circular_qwp_set": np.array(circulars)},
    references=["batch_20260702_110714_779653"],
    script=open(__file__).read(),
    figures=["scanA_heatmap.png", "scanA_qwp_fit.png"],
)
np.savez(f"{SCR}/scanA_fit.npz", qwp=qwp, pol=pol, amp=amp, dc=dc, phase=phase,
         o_q=o_q_fit, o_q_err=o_q_err, P0=P0_est, Cbg=Cbg_est,
         neutrals=neutrals, circulars=circulars)
print("DONE")

"""Scan B (acquisition only): 2D scan hwp_set x pol_set, QWP parked at its
neutral (qwp_set = 11.96 -> qwp_true = 0, horizontal linear into the HWP).

With linear light, the analyzer peak sits at pol_true = 2*hwp_true, i.e.
    pol_set_peak = 2*hwp_set + C,   C = 2*o_h - o_p.
Power depends only on u = pol_set - 2*hwp_set (=> o_h/o_p degeneracy demo).
"""
import sys
sys.path.insert(0, "/agent/shared/scripts"); sys.path.insert(0, "/agent/workspace")
import numpy as np
from tools_client import set_angle, measure_power
from auto_log_client import start_batch, stop_batch

SCR = "/tmp/claude-1001/-agent-workspace/5b062f49-458b-430d-8d6c-7c8ec035f7d0/scratchpad"
def P():
    d = measure_power()["power"]
    return float(d["value"]) if isinstance(d, dict) else float(d)

QWP_NEUTRAL = 11.96          # from Scan A (qwp_true = 0)
hwp_ax = np.arange(0, 180, 5.0)     # 36
pol_ax = np.arange(0, 180, 15.0)    # 12

set_angle(QWP_NEUTRAL, "lambda_quarter")
power = np.full((hwp_ax.size, pol_ax.size), np.nan)

ret = start_batch("Scan B: hwp(0:5:180) x pol(0:15:180), qwp=neutral(11.96)",
                  "2D scan for C=2o_h-o_p via analyzer peak vs HWP; degeneracy demo.")
print("start_batch:", ret)
for i, h in enumerate(hwp_ax):
    set_angle(float(h), "lambda_half")
    for j, p in enumerate(pol_ax):
        set_angle(float(p), "polarizer")
        power[i, j] = P()
stop = stop_batch()
print("stop_batch:", stop)

np.savez(f"{SCR}/scanB.npz", hwp_ax=hwp_ax, pol_ax=pol_ax, power=power, qwp_neutral=QWP_NEUTRAL)
print("saved scanB.npz shape", power.shape,
      " Pmin=%.4f Pmax=%.4f" % (np.nanmin(power), np.nanmax(power)))
print("DONE")

"""Analyze Scan B (hwp x pol, qwp=neutral) -> C = 2*o_h - o_p, and demonstrate
the o_h/o_p degeneracy.

Per-hwp harmonic fit gives analyzer peak phi(h). Model: phi = 2*h + C (mod 180),
so the complex c=a1+i a2 has angle = deg2rad(4h + 2C): a LINEAR unwrap fit gives
slope (expect +4 deg/deg, i.e. peak slope +2) and intercept 2C.

Degeneracy demo: power depends only on u = pol_set - 2*hwp_set (since
pol_true - 2*hwp_true = u - C). Collapse all 432 points onto u (mod 180); a
single sinusoid fit gives C from its peak and the collapse RMS proves the
'only 2*o_h - o_p matters' structure.
"""
import sys
sys.path.insert(0, "/agent/shared/scripts"); sys.path.insert(0, "/agent/workspace")
import numpy as np, matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from auto_log_client import log_analysis, AUTO_LOG_DIR

SCR = "/tmp/claude-1001/-agent-workspace/5b062f49-458b-430d-8d6c-7c8ec035f7d0/scratchpad"
SIGMA = 0.005
d = np.load(f"{SCR}/scanB.npz")
hwp, pol, power = d["hwp_ax"], d["pol_ax"], d["power"]

# ---- per-row harmonic fit -------------------------------------------------
th = np.deg2rad(pol)
Xd = np.column_stack([np.ones_like(th), np.cos(2*th), np.sin(2*th)])
coef, *_ = np.linalg.lstsq(Xd, power.T, rcond=None)
a0, a1, a2 = coef
amp = np.hypot(a1, a2); dc = a0
phi = np.rad2deg(0.5*np.arctan2(a2, a1)) % 180.0    # analyzer peak set-angle

# ---- slope fit: angle(c)=deg2rad(4h+2C) ----------------------------------
c = a1 + 1j*a2
ang = np.unwrap(np.angle(c))               # radians
slope, intercept = np.polyfit(hwp, np.rad2deg(ang), 1)   # deg per deg, deg
C_slope = (intercept/2.0) % 180.0
peak_slope = slope/2.0                      # expect +2

# ---- degeneracy collapse: P vs u = pol_set - 2*hwp_set -------------------
H, Pp = np.meshgrid(hwp, pol, indexing="ij")
u = (Pp - 2*H)                                       # (36,12)
u_mod = u % 180.0
uf, pf = u_mod.ravel(), power.ravel()
def col(u_, A, B, Cc):
    return A + B*np.cos(2*np.deg2rad(u_ - Cc))
p0 = [pf.mean(), (pf.max()-pf.min())/2, uf[np.argmax(pf)]]
popt, pcov = curve_fit(col, uf, pf, p0=p0, sigma=np.full_like(pf, SIGMA), absolute_sigma=True)
A_c, B_c, C_col = popt; C_col %= 180.0
C_err = float(np.sqrt(pcov[2, 2]))
collapse_rms = float(np.sqrt(np.mean((pf - col(uf, *popt))**2)))
P0_B = 2*abs(B_c)

print("amp(h): mean=%.4f std=%.4f (should be ~const=P0/2)" % (amp.mean(), amp.std()))
print("peak-vs-hwp slope = %.4f (expect +2.000)" % peak_slope)
print("C (slope fit)     = %.3f deg (mod 180)" % C_slope)
print("C (collapse fit)  = %.3f +/- %.3f deg (mod 180) ; P0=%.4f" % (C_col, C_err, P0_B))
print("collapse RMS = %.5f W  vs noise sigma=%.4f -> degeneracy holds: %s"
      % (collapse_rms, SIGMA, collapse_rms < 3*SIGMA))

# ---- figures --------------------------------------------------------------
fig, ax = plt.subplots(figsize=(7.5,4.2))
im = ax.pcolormesh(pol, hwp, power, shading="auto", cmap="magma")
# overlay predicted peak line pol = 2*hwp + C (mod 180)
hh = np.linspace(0,180,400)
for k in (-2,-1,0,1,2):
    ax.plot((2*hh + C_col + 180*k), hh, "c--", lw=1)
ax.set_xlim(pol.min(), pol.max()); ax.set_ylim(hwp.min(), hwp.max())
ax.set_xlabel("polarizer set-angle (deg)"); ax.set_ylabel("HWP set-angle (deg)")
ax.set_title("Scan B: power(HWP, polarizer), QWP neutral\n(peak line slope +2: dashed)")
fig.colorbar(im, label="power (W)"); fig.tight_layout()
fig.savefig(f"{AUTO_LOG_DIR}/scanB_heatmap.png", dpi=110); plt.close(fig)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11,4.2))
ax1.plot(hwp, phi, "o", ms=4)
ax1.set_xlabel("HWP set-angle (deg)"); ax1.set_ylabel("analyzer peak (deg, mod 180)")
ax1.set_title(f"Analyzer peak vs HWP: slope={peak_slope:.3f} (≈+2)  =>  C={C_slope:.2f}°")
ax2.plot(uf, pf, ".", ms=5, alpha=0.6, label="all 432 points")
us = np.linspace(0,180,400); ax2.plot(us, col(us,*popt), "r-", lw=2,
        label=f"single sinusoid, peak C={C_col:.2f}°")
ax2.set_xlabel("u = pol_set − 2·HWP_set  (mod 180°)"); ax2.set_ylabel("power (W)")
ax2.set_title(f"DEGENERACY: power depends only on u  (collapse RMS={collapse_rms:.4f} W ≈ σ)")
ax2.legend(fontsize=8)
fig.tight_layout(); fig.savefig(f"{AUTO_LOG_DIR}/scanB_C_and_degeneracy.png", dpi=110); plt.close(fig)

log_analysis(
    title="Scan B -> C = 2·o_h − o_p = %.2f° (mod 180); o_h,o_p individually degenerate" % C_col,
    kind="analysis",
    text=(
        f"QWP parked at neutral (qwp_true=0): light into the HWP is horizontal linear, "
        f"extinction is essentially perfect (Pmin≈0), confirming C_bg≈0. Two methods for "
        f"C=2·o_h−o_p agree: analyzer-peak-vs-HWP slope = {peak_slope:.3f} (predicted +2, ✓) "
        f"with intercept giving C={C_slope:.2f}°; and the global collapse fit "
        f"C={C_col:.2f}±{C_err:.2f}° (mod 180). Sweep amplitude is flat (mean {amp.mean():.4f}, "
        f"std {amp.std():.4f} = P0/2), i.e. light stays fully linear for all HWP — as expected.\n\n"
        f"DEGENERACY DEMONSTRATED: all 432 (hwp,pol) points collapse onto a single sinusoid in "
        f"u = pol_set − 2·HWP_set with RMS {collapse_rms:.4f} W ≈ detector σ ({SIGMA} W). The power "
        f"is a function of (2·HWP − pol) ALONE, so only the combination C=2·o_h−o_p is observable; "
        f"o_h and o_p cannot be separated by intensity measurements (HWP+polarizer act as one "
        f"effective analyzer at 2·hwp_true − pol_true). C is defined mod 180°."),
    data={"hwp_set": hwp, "pol_set": pol, "amp": amp, "phi_peak": phi,
          "C_collapse_deg": C_col, "C_err_deg": C_err, "C_slope_deg": C_slope,
          "peak_slope": peak_slope, "collapse_rms": collapse_rms, "P0_B": P0_B},
    references=["batch_20260702_110949_886033"],
    script=open(__file__).read(),
    figures=["scanB_heatmap.png", "scanB_C_and_degeneracy.png"],
)
np.savez(f"{SCR}/scanB_fit.npz", hwp=hwp, pol=pol, amp=amp, phi=phi,
         C=C_col, C_err=C_err, C_slope=C_slope, peak_slope=peak_slope,
         collapse_rms=collapse_rms)
print("DONE")

"""Global joint fit of the full Mueller forward model to ALL raw scan data
(Scan A: qwp x pol @ hwp=0 ; Scan B: hwp x pol @ qwp=neutral), 864 single-shot
points. Because o_h and o_p are degenerate (only C=2o_h-o_p is observable) we
FIX THE GAUGE o_h := 0 and fit o_p (=> C = -o_p). Free parameters:

  ideal model (primary):   o_q, o_p, P0, Cbg            (waveplates ideal)
  extended model (check):  o_q, o_p, P0, Cbg, ret_q, ret_h

Reports offsets + uncertainties, reduced chi^2, residual RMS vs noise, and a
data-vs-model validation figure.
"""
import sys
sys.path.insert(0, "/agent/shared/scripts"); sys.path.insert(0, "/agent/workspace")
import numpy as np, matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
from physics_model import predict_power
from auto_log_client import log_analysis, AUTO_LOG_DIR

SCR = "/tmp/claude-1001/-agent-workspace/5b062f49-458b-430d-8d6c-7c8ec035f7d0/scratchpad"
SIGMA = 0.005

# ---- assemble combined dataset -------------------------------------------
A = np.load(f"{SCR}/scanA.npz"); B = np.load(f"{SCR}/scanB.npz")
qA, pA, PA = A["qwp_ax"], A["pol_ax"], A["power"]
hB, pB, PB = B["hwp_ax"], B["pol_ax"], B["power"]; qN = float(B["qwp_neutral"])

# Scan A grid (hwp=0)
QA, PPA = np.meshgrid(qA, pA, indexing="ij")
q_all = list(QA.ravel());          h_all = [0.0]*QA.size;         p_all = list(PPA.ravel()); P_all = list(PA.ravel())
# Scan B grid (qwp=neutral)
HB, PPB = np.meshgrid(hB, pB, indexing="ij")
q_all += [qN]*HB.size;             h_all += list(HB.ravel());     p_all += list(PPB.ravel()); P_all += list(PB.ravel())
q_all = np.array(q_all); h_all = np.array(h_all); p_all = np.array(p_all); P_all = np.array(P_all)
print("combined points:", P_all.size)

def resid(params, free_ret=False):
    if free_ret:
        o_q, o_p, P0, Cbg, rq, rh = params
    else:
        o_q, o_p, P0, Cbg = params; rq, rh = 90.0, 180.0
    model = predict_power(q_all, h_all, p_all, o_q=o_q, o_h=0.0, o_p=o_p,
                          P0=P0, Cbg=Cbg, ret_q=rq, ret_h=rh)
    return (model - P_all)/SIGMA

def fit(p0, free_ret=False, bounds=(-np.inf, np.inf)):
    sol = least_squares(resid, p0, args=(free_ret,), method="trf", bounds=bounds, x_scale="jac")
    J = sol.jac; dof = P_all.size - len(p0)
    chi2 = float(np.sum(sol.fun**2)); redchi2 = chi2/dof
    # covariance (residuals already /sigma): cov = (J^T J)^-1
    cov = np.linalg.inv(J.T @ J)
    perr = np.sqrt(np.diag(cov))
    rms = float(np.sqrt(np.mean((sol.fun*SIGMA)**2)))
    return sol.x, perr, redchi2, rms

# primary ideal-waveplate fit
p0 = [78.0, -73.0, 1.0, 0.0]
x, e, rc, rms = fit(p0, free_ret=False)
o_q, o_p, P0, Cbg = x
C = (-o_p) % 180.0
oq_mod = o_q % 90.0
print("\n=== IDEAL model ===")
print("o_q  = %.3f ± %.3f deg (mod 90)  -> %.3f" % (o_q, e[0], oq_mod))
print("o_p  = %.3f ± %.3f deg (gauge o_h=0) -> C=2o_h-o_p = %.3f (mod 180)" % (o_p, e[1], C))
print("P0   = %.4f ± %.4f W" % (P0, e[2]))
print("Cbg  = %.5f ± %.5f W" % (Cbg, e[3]))
print("reduced chi2 = %.3f   residual RMS = %.5f W  (noise σ=%.4f)" % (rc, rms, SIGMA))

# extended fit with free retardances
p0e = [o_q, o_p, P0, Cbg, 90.0, 180.0]
xe, ee, rce, rmse = fit(p0e, free_ret=True)
print("\n=== EXTENDED model (free retardances) ===")
print("o_q=%.3f±%.3f  o_p=%.3f±%.3f  P0=%.4f  Cbg=%.5f" % (xe[0],ee[0],xe[1],ee[1],xe[2],xe[3]))
print("ret_q = %.3f ± %.3f deg (ideal 90)   ret_h = %.3f ± %.3f deg (ideal 180)"
      % (xe[4], ee[4], xe[5], ee[5]))
print("reduced chi2 = %.3f   residual RMS = %.5f W" % (rce, rmse))

# ---- validation figure ----------------------------------------------------
model_all = predict_power(q_all, h_all, p_all, o_q=o_q, o_h=0.0, o_p=o_p, P0=P0, Cbg=Cbg)
res = P_all - model_all
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11,4.3))
ax1.plot(P_all, model_all, ".", ms=4, alpha=0.4)
ax1.plot([0,1.05],[0,1.05], "r-", lw=1)
ax1.set_xlabel("measured power (W)"); ax1.set_ylabel("model power (W)")
ax1.set_title(f"Global fit: measured vs model (864 pts)\nreduced χ²={rc:.2f}, RMS={rms:.4f} W ≈ σ")
ax2.hist(res, bins=40, color="steelblue", edgecolor="k", alpha=0.8)
ax2.axvline(0, color="k", lw=0.8)
ax2.set_xlabel("residual (W)"); ax2.set_ylabel("count")
ax2.set_title(f"Residuals: mean={res.mean():.4f}, std={res.std():.4f} W (σ={SIGMA})")
fig.tight_layout(); fig.savefig(f"{AUTO_LOG_DIR}/global_fit_validation.png", dpi=110); plt.close(fig)

log_analysis(
    title="Global fit: o_q=%.2f° (mod90), C=2o_h−o_p=%.2f° (mod180), P0=%.3f W, Cbg≈0" % (oq_mod, C, P0),
    kind="analysis",
    text=(
        f"Joint least-squares fit of the full Mueller model to all 864 single-shot "
        f"points (Scan A + Scan B), gauge o_h:=0.\n"
        f"IDEAL waveplates: o_q={o_q:.3f}±{e[0]:.3f}° (mod 90 = {oq_mod:.2f}°), "
        f"o_p={o_p:.3f}±{e[1]:.3f}° => C=2·o_h−o_p={C:.3f}° (mod 180), "
        f"P0={P0:.4f}±{e[2]:.4f} W, C_bg={Cbg:.5f}±{e[3]:.5f} W (consistent with 0). "
        f"Reduced χ²={rc:.2f}, residual RMS={rms:.5f} W ≈ detector σ ({SIGMA} W): the ideal "
        f"model explains the data to the noise floor.\n"
        f"EXTENDED fit (free retardances): ret_q={xe[4]:.2f}±{ee[4]:.2f}° (ideal 90), "
        f"ret_h={xe[5]:.2f}±{ee[5]:.2f}° (ideal 180) — the waveplates are ideal within error, "
        f"so no need to correct for imperfect retardance.\n"
        f"These global values agree with the per-scan results (o_q 78.04°, C 73.05°) and reproduce "
        f"the independent qwp=hwp=0 sweep peak (85.0°). Reminder: o_q is mod 90° (fast/slow "
        f"unresolved by intensity), and only C is physical — o_h, o_p are individually gauge."),
    data={"o_q_deg": o_q, "o_q_err": e[0], "o_q_mod90": oq_mod,
          "o_p_gauge_oh0_deg": o_p, "o_p_err": e[1], "C_deg": C,
          "P0": P0, "P0_err": e[2], "Cbg": Cbg, "Cbg_err": e[3],
          "redchi2": rc, "resid_rms": rms,
          "ret_q": xe[4], "ret_q_err": ee[4], "ret_h": xe[5], "ret_h_err": ee[5],
          "residuals": res},
    references=["batch_20260702_110714_779653", "batch_20260702_110949_886033"],
    script=open(__file__).read(),
    figures=["global_fit_validation.png"],
)
np.savez(f"{SCR}/global_fit.npz", o_q=o_q, o_q_err=e[0], o_p=o_p, C=C, P0=P0, Cbg=Cbg,
         redchi2=rc, rms=rms, ret_q=xe[4], ret_h=xe[5])
print("DONE")

"""Master calibration summary figure (2x3) + final ELN summary entry.

Panels: (1) Scan A heatmap, (2) QWP amplitude fit -> o_q, (3) Scan B heatmap,
(4) degeneracy collapse -> C, (5) global-fit validation, (6) results text.
"""
import sys
sys.path.insert(0, "/agent/shared/scripts"); sys.path.insert(0, "/agent/workspace")
import numpy as np, matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from physics_model import predict_power
from auto_log_client import log_analysis, AUTO_LOG_DIR

SCR = "/tmp/claude-1001/-agent-workspace/5b062f49-458b-430d-8d6c-7c8ec035f7d0/scratchpad"
A = np.load(f"{SCR}/scanA.npz"); B = np.load(f"{SCR}/scanB.npz")
fa = np.load(f"{SCR}/scanA_fit.npz"); fb = np.load(f"{SCR}/scanB_fit.npz")
g  = np.load(f"{SCR}/global_fit.npz")
o_q = float(g["o_q"]); o_q_err = float(fa["o_q_err"]); C = float(g["C"]); C_err = float(fb["C_err"])
P0 = float(g["P0"]); Cbg = float(g["Cbg"]); redchi2 = float(g["redchi2"]); rms = float(g["rms"])
ret_q = float(g["ret_q"]); ret_h = float(g["ret_h"])

plt.rcParams.update({"font.size": 10})
fig, axes = plt.subplots(2, 3, figsize=(15.5, 9.2))

# (1) Scan A heatmap
ax = axes[0,0]
im = ax.pcolormesh(A["pol_ax"], A["qwp_ax"], A["power"], shading="auto", cmap="viridis")
for qn in fa["neutrals"]: ax.axhline(qn, color="r", ls="--", lw=1)
for qc in fa["circulars"]: ax.axhline(qc, color="w", ls=":", lw=1)
ax.set_xlabel("polarizer set-angle (°)"); ax.set_ylabel("QWP set-angle (°)")
ax.set_title("① Scan A — power(QWP, pol), HWP=0\nred=linear(neutral), white=circular")
fig.colorbar(im, ax=ax, label="P (W)")

# (2) QWP amplitude fit
ax = axes[0,1]
amp, qax = fa["amp"], fa["qwp"]
ax.plot(qax, amp, "o", ms=4, color="#1f77b4", label="sweep amplitude")
xs = np.linspace(0,180,500)
ax.plot(xs, (P0/2)*np.abs(np.cos(np.deg2rad(2*(xs+o_q)))), "r-",
        label=r"$(P_0/2)\,|\cos 2(q+o_q)|$")
for qn in fa["neutrals"]: ax.axvline(qn, color="r", ls="--", lw=1)
ax.set_xlabel("QWP set-angle (°)"); ax.set_ylabel("amplitude (W)")
ax.set_title(f"② QWP calibration → $o_q$={o_q:.2f}±{o_q_err:.2f}° (mod 90)")
ax.legend(fontsize=8, loc="upper right")

# (3) Scan B heatmap
ax = axes[0,2]
im = ax.pcolormesh(B["pol_ax"], B["hwp_ax"], B["power"], shading="auto", cmap="magma")
hh = np.linspace(0,180,400)
for k in (-2,-1,0,1,2): ax.plot(2*hh + C + 180*k, hh, "c--", lw=1)
ax.set_xlim(B["pol_ax"].min(), B["pol_ax"].max()); ax.set_ylim(B["hwp_ax"].min(), B["hwp_ax"].max())
ax.set_xlabel("polarizer set-angle (°)"); ax.set_ylabel("HWP set-angle (°)")
ax.set_title("③ Scan B — power(HWP, pol), QWP neutral\ncyan: peak line slope +2")
fig.colorbar(im, ax=ax, label="P (W)")

# (4) degeneracy collapse
ax = axes[1,0]
H, Pp = np.meshgrid(B["hwp_ax"], B["pol_ax"], indexing="ij")
u = ((Pp - 2*H) % 180.0).ravel(); pw = B["power"].ravel()
def col(u_, Aa, Bb, Cc): return Aa + Bb*np.cos(2*np.deg2rad(u_ - Cc))
popt,_ = curve_fit(col, u, pw, p0=[0.5,0.5,C])
ax.plot(u, pw, ".", ms=5, alpha=0.5, color="#6a3d9a", label="all 432 points")
us = np.linspace(0,180,400); ax.plot(us, col(us,*popt), "r-", lw=2, label="single sinusoid")
ax.set_xlabel(r"$u=\mathrm{pol}_{set}-2\,\mathrm{HWP}_{set}$ (mod 180°)"); ax.set_ylabel("P (W)")
ax.set_title(f"④ Degeneracy: P depends only on u\n→ C=2$o_h$−$o_p$={C:.2f}±{C_err:.2f}° (RMS≈σ)")
ax.legend(fontsize=8)

# (5) global fit validation
ax = axes[1,1]
qA = A["qwp_ax"]; pA = A["pol_ax"]; QA,PPA = np.meshgrid(qA,pA,indexing="ij")
q_all = np.concatenate([QA.ravel(), np.full(H.size, float(B["qwp_neutral"]))])
h_all = np.concatenate([np.zeros(QA.size), H.ravel()])
p_all = np.concatenate([PPA.ravel(), Pp.ravel()])
P_all = np.concatenate([A["power"].ravel(), B["power"].ravel()])
model = predict_power(q_all, h_all, p_all, o_q=o_q, o_h=0.0, o_p=-C, P0=P0, Cbg=Cbg)
ax.plot(P_all, model, ".", ms=4, alpha=0.4, color="#2ca02c")
ax.plot([0,1.05],[0,1.05], "r-", lw=1)
ax.set_xlabel("measured P (W)"); ax.set_ylabel("model P (W)")
ax.set_title(f"⑤ Global fit validation (864 pts)\nred. χ²={redchi2:.2f}, RMS={rms:.4f} W ≈ σ")

# (6) results text
ax = axes[1,2]; ax.axis("off")
txt = (
    "CALIBRATION RESULTS\n"
    "(true angle = set angle + offset;\n 0° = horizontal, from fixed input polarizer)\n"
    "──────────────────────────────\n"
    f"λ/4  QWP:   $o_q$ = {o_q:.2f} ± {o_q_err:.2f}°  (mod 90°)\n"
    f"   fast axis ∥ H at set ≈ {(-o_q)%90:.2f}°\n"
    f"   (fast/slow unresolved → 78° or 168°)\n\n"
    f"λ/2 HWP  &  polarizer:  DEGENERATE\n"
    f"   only  C = 2·$o_h$ − $o_p$ = {C:.2f} ± {C_err:.2f}°\n"
    f"   (mod 180°) is observable.\n"
    f"   e.g.  $o_h$=0 → $o_p$=−{C:.1f}° ;\n"
    f"         $o_p$=0 → $o_h$=+{C/2:.1f}°\n\n"
    "SYSTEM\n"
    f"   $P_0$ = {P0:.3f} W,  background ≈ 0\n"
    f"   noise: additive σ ≈ 0.005 W\n"
    f"   retardances: QWP {ret_q:.1f}°, HWP {ret_h:.1f}°\n"
    f"   (ideal within ±0.1°)\n"
    f"   reduced χ² = {redchi2:.2f}"
)
ax.text(0.0, 0.98, txt, va="top", ha="left", family="monospace", fontsize=10.3,
        bbox=dict(boxstyle="round", fc="#f4f4f4", ec="gray"))

fig.suptitle("Absolute-angle calibration of QWP, HWP and polarizer  "
             "(H-pol → QWP → HWP → polarizer → detector)", fontsize=13, y=0.995)
fig.tight_layout(rect=[0,0,1,0.97])
fig.savefig(f"{AUTO_LOG_DIR}/CALIBRATION_SUMMARY.png", dpi=120); plt.close(fig)
print("saved CALIBRATION_SUMMARY.png")

log_analysis(
    title="CALIBRATION SUMMARY — QWP o_q=78.0° (mod90); HWP/pol combination C=73.0° (mod180)",
    kind="analysis",
    text=(
        "FINAL absolute-angle calibration (convention: true = set + offset, 0° = horizontal "
        "defined by the fixed input polarizer).\n\n"
        f"• λ/4 QWP: o_q = {o_q:.2f} ± {o_q_err:.2f}° (mod 90°). Fast axis is horizontal when the "
        f"QWP is set to ≈{(-o_q)%90:.2f}° (or +90°). Only mod 90° is knowable from intensity "
        f"(handedness/fast-vs-slow is invisible to a power meter), so the fast-axis offset is "
        f"either 78° or 168° (mod 180°).\n\n"
        f"• λ/2 HWP and rotatable polarizer are FUNDAMENTALLY DEGENERATE with a power meter: they "
        f"form a single effective analyzer at 2·hwp_true − pol_true, so only the combination "
        f"C = 2·o_h − o_p = {C:.2f} ± {C_err:.2f}° (mod 180°) is measurable. Demonstrated by the "
        f"collapse of all 432 (HWP,pol) points onto u = pol−2·HWP at noise level. To assign "
        f"individual offsets one extra assumption is needed (e.g. o_h=0 ⇒ o_p=−{C:.1f}°, or "
        f"o_p=0 ⇒ o_h=+{C/2:.1f}°); breaking it physically requires a polarization-resolving "
        f"detector or a known reference polarization.\n\n"
        f"• System: P0={P0:.3f} W, background≈0, additive noise σ≈0.005 W, QWP/HWP retardances "
        f"{ret_q:.1f}°/{ret_h:.1f}° (ideal within ±0.1°). Global Mueller fit: reduced χ²={redchi2:.2f}, "
        f"residual RMS={rms:.4f} W = σ. All independent methods (per-scan + global + the "
        f"qwp=hwp=0 cross-check at 85.0°) agree."),
    data={"o_q_deg": o_q, "o_q_err_deg": o_q_err, "o_q_mod": 90,
          "C_2oh_minus_op_deg": C, "C_err_deg": C_err, "C_mod": 180,
          "P0_W": P0, "Cbg_W": Cbg, "sigma_noise_W": 0.005,
          "ret_q_deg": ret_q, "ret_h_deg": ret_h, "reduced_chi2": redchi2, "resid_rms_W": rms,
          "qwp_fast_axis_horizontal_set_deg": (-o_q) % 90},
    references=["batch_20260702_110714_779653", "batch_20260702_110949_886033"],
    script=open(__file__).read(),
    figures=["CALIBRATION_SUMMARY.png"],
)
print("DONE")

cd /agent/shared/auto_log && ls -1 analysis_*.json 2>/dev/null | tail -3; echo "--- inspect latest summary analysis record ---"; f=$(ls -1t analysis_*.json 2>/dev/null | head -1); echo "file: $f"; python3 -c "
import json
d=json.load(open('$f'))
print('title:', d.get('title'))
print('figures field:', d.get('figures'))
print('keys:', list(d.keys()))
"; echo "--- any copied figure dirs? ---"; ls -1 | grep -iE 'fig|CALIBRATION' | head

log_analysis(
    title="CALIBRATION SUMMARY — QWP o_q=78.0° (mod90); HWP/pol combination C=73.0° (mod180)",

if "--nolog" in sys.argv:
    print("skip logging (figure regenerated in place; existing ELN record now points to it)")
    sys.exit(0)

log_analysis(
    title="CALIBRATION SUMMARY — QWP o_q=78.0° (mod90); HWP/pol combination C=73.0° (mod180)",

"""Blind hold-out validation of the calibration.

Pick random (qwp, hwp, pol) triples over the FULL 0-360 deg range (the scans
only used 0-180 and were on coarse grids, so these are genuinely held out).
Predict each power from the calibrated Mueller model and compare to a fresh
averaged measurement. If predictions match at the noise level, the calibration
(o_q, C, P0, Cbg + ideal waveplates) is validated end-to-end.

Gauge used: o_h=0, o_p=-C  (any gauge with the same C predicts identical power).
"""
import sys
sys.path.insert(0, "/agent/shared/scripts"); sys.path.insert(0, "/agent/workspace")
import numpy as np, matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from tools_client import set_angle, measure_power
from auto_log_client import start_batch, stop_batch, log_analysis, AUTO_LOG_DIR
from physics_model import predict_power

SCR = "/tmp/claude-1001/-agent-workspace/5b062f49-458b-430d-8d6c-7c8ec035f7d0/scratchpad"
SIGMA = 0.005
g = np.load(f"{SCR}/global_fit.npz")
o_q, C, P0, Cbg = float(g["o_q"]), float(g["C"]), float(g["P0"]), float(g["Cbg"])

def P():
    d = measure_power()["power"]
    return float(d["value"]) if isinstance(d, dict) else float(d)

rng = np.random.default_rng(20260702)
N = 16
Q = rng.uniform(0, 360, N); H = rng.uniform(0, 360, N); Pp = rng.uniform(0, 360, N)
# add 3 diagnostic configs: QWP neutral + HWP that rotates H->V + crossed/aligned analyzer
# with o_h=0 gauge, hwp_true=hwp_set; light after HWP linear at 2*hwp_set; peak at pol_true=2*hwp_set
diag = np.array([
    [11.96,  0.0,   -C % 360],          # neutral QWP, HWP=0 -> light H -> analyzer aligned (max)
    [11.96,  0.0,   (90 - C) % 360],     # same but analyzer crossed (min ~ 0)
    [11.96, 22.5,   (45 - C) % 360],     # HWP rotates 45 deg -> analyzer aligned (max)
])
Q = np.concatenate([Q, diag[:,0]]); H = np.concatenate([H, diag[:,1]]); Pp = np.concatenate([Pp, diag[:,2]])

pred = predict_power(Q, H, Pp, o_q=o_q, o_h=0.0, o_p=-C, P0=P0, Cbg=Cbg)

start_batch("Hold-out validation: 19 random/diagnostic (qwp,hwp,pol) triples (avg 5)")
meas = np.empty(Q.size)
for i in range(Q.size):
    set_angle(float(Q[i]),  "lambda_quarter")
    set_angle(float(H[i]),  "lambda_half")
    set_angle(float(Pp[i]), "polarizer")
    meas[i] = np.mean([P() for _ in range(5)])
stop_batch()

resid = meas - pred
rms = float(np.sqrt(np.mean(resid**2))); maxerr = float(np.max(np.abs(resid)))
sigma_mean = SIGMA/np.sqrt(5)
print("hold-out points:", Q.size)
for i in range(Q.size):
    tag = " <-diag" if i >= N else ""
    print("  q=%6.1f h=%6.1f p=%6.1f | pred=%.4f meas=%.4f d=%+.4f%s"
          % (Q[i], H[i], Pp[i], pred[i], meas[i], resid[i], tag))
print("RMS=%.5f W  max|err|=%.5f W  (single-shot σ=%.4f, avg5 σ=%.4f)"
      % (rms, maxerr, SIGMA, sigma_mean))

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11,4.3))
ax1.errorbar(pred, meas, yerr=sigma_mean, fmt="o", ms=6, capsize=3, color="#d62728", label="held-out")
ax1.plot([0,1.05],[0,1.05], "k-", lw=1, label="ideal y=x")
ax1.set_xlabel("PREDICTED power from calibration (W)"); ax1.set_ylabel("MEASURED power (W)")
ax1.set_title(f"Blind hold-out validation (19 configs, full 0–360°)\nRMS={rms:.4f} W ≈ σ/√5={sigma_mean:.4f}")
ax1.legend(fontsize=8)
ax2.axhline(0, color="k", lw=0.8)
ax2.errorbar(np.arange(Q.size), resid, yerr=sigma_mean, fmt="o", ms=5, capsize=2, color="#1f77b4")
ax2.axhspan(-2*sigma_mean, 2*sigma_mean, color="gray", alpha=0.2, label="±2σ band")
ax2.set_xlabel("configuration #"); ax2.set_ylabel("measured − predicted (W)")
ax2.set_title("Residuals within noise"); ax2.legend(fontsize=8)
fig.tight_layout(); fig.savefig(f"{AUTO_LOG_DIR}/holdout_validation.png", dpi=120); plt.close(fig)

log_analysis(
    title="Hold-out validation PASSED — calibration predicts blind configs to σ level (RMS=%.4f W)" % rms,
    kind="analysis",
    text=(
        f"Blind test on 19 (qwp,hwp,pol) triples over the full 0–360° range (16 random + 3 "
        f"diagnostic), none on the calibration grids. Predicting power from the calibrated model "
        f"(o_q={o_q:.2f}°, gauge o_h=0 ⇒ o_p=−C=−{C:.2f}°, P0={P0:.3f} W, Cbg={Cbg:.4f} W, ideal "
        f"waveplates) vs fresh measurements (avg 5): RMS={rms:.5f} W, max|err|={maxerr:.5f} W — "
        f"both at the averaged-noise level (σ/√5={sigma_mean:.4f} W). The diagnostic configs behaved "
        f"as designed (aligned analyzer → ~P0, crossed → ~0). This confirms the full forward model "
        f"and the calibrated offsets over the entire operating range; the o_q/C parametrization is "
        f"sufficient to predict any measurement (individual o_h,o_p need not — and cannot — be known)."),
    data={"qwp": Q, "hwp": H, "pol": Pp, "predicted": pred, "measured": meas,
          "residuals": resid, "rms_W": rms, "max_abs_err_W": maxerr, "sigma_mean_W": sigma_mean},
    references=["batch_20260702_110714_779653", "batch_20260702_110949_886033"],
    script=open(__file__).read(),
    figures=["holdout_validation.png"],
)
np.savez(f"{SCR}/holdout.npz", Q=Q, H=H, Pp=Pp, pred=pred, meas=meas, rms=rms, maxerr=maxerr)
print("DONE")

# Optics setup — calibration reference

Beam path: **Laser → fixed horizontal polarizer → λ/4 (QWP) → λ/2 (HWP) → rotatable polarizer → detector.**
Tools: `set_angle(angle, component)` for `lambda_quarter` / `lambda_half` / `polarizer`; `measure_power()` → `{'power': {'value': W, 'unit': 'W'}}`.

## Angle convention
`true_angle = set_angle + offset`. **0° = horizontal**, the axis defined by the fixed input polarizer.
For waveplates, true 0° means the fast axis is horizontal. To command a true angle: `set = true − offset`.

## System constants
- **P₀ ≈ 1.00 W** — full power reaching the detector (fully transmitted).
- **Background ≈ 0** — extinction goes to ~0 W with linear light; there is no dark offset.
- **Detector noise: additive Gaussian, σ ≈ 0.005 W, independent of power level.** (So relative noise is ~0.5 % near full scale but ~10 % near extinction. Average N shots → σ/√N. Use uniform-weight least squares.)
- **Both waveplates are ideal:** retardances 90.1° (λ/4) and 179.9° (λ/2).

## Calibrated offsets
### λ/4 QWP — fully calibrated (mod 90°)
- **offset o_q = 78.0° (mod 90°)**, ±0.03°.
- Fast/slow axis is **horizontal at set-angle ≈ 12.0°** (and 102.0°).
- Makes **circular light at set-angle ≈ 57.0°** (and 147.0°).
- ⚠️ Only known **mod 90°**: a power meter cannot sense handedness, so fast-vs-slow is unresolved (the fast-axis offset is 78° *or* 168°). This does not affect intensity predictions.

### λ/2 HWP + rotatable polarizer — degenerate (only the combination is measurable)
- ⚠️ **The HWP and analyzing polarizer offsets cannot be separated with a power meter.** They act as one effective analyzer at `2·hwp_true − pol_true`, so only the combination is observable:
  - **C = 2·o_h − o_p = 73.0° (mod 180°)**, ±0.03°.
- Power depends on the HWP and polarizer settings **only through `u = pol_set − 2·hwp_set`**.
- To assign individual numbers you must pick a gauge, e.g. `o_h = 0 → o_p = −73.0°`, or `o_p = 0 → o_h = +36.5°`. Separating them physically needs a polarization-resolving detector or a known reference polarization.

## Ready-to-use operating points (gauge-free, in set-angles)
- **Preserve horizontal linear light** (QWP acts neutral): `lambda_quarter = 12.0°` (or 102.0°).
- **Circular light after QWP:** `lambda_quarter = 57.0°` (or 147.0°).
- **With the QWP neutral**, light after the HWP is linear at true angle `2·hwp_true`, analyzed by the polarizer. The detected power follows Malus's law in `u = pol_set − 2·hwp_set`:
  - **P ≈ P₀ · cos²(u − 73.0°)**, i.e. **maximum at `u = 73.0°`**, **extinction at `u = 163.0°`** (both mod 180°).
- **Rotate the output linear polarization** by δ (true): increase `lambda_half` by δ/2 (HWP rotates polarization at 2× its own rotation).

## Confidence
Global Mueller-model fit over 864 points: reduced χ² = 1.01, residual RMS = σ. Blind hold-out test over the full 0–360° range: predicted vs measured RMS = 0.002 W (noise-limited). Full raw data, analysis, figures and scripts are in `/agent/shared/auto_log/` (see `CALIBRATION_SUMMARY.png`).