"""
n15 — Is the vacuum's record resolution R0 derivable? Replace pointwise
monitoring by smeared Gaussian probes of width R; for each R bisect the
classicality threshold kappa*(R) and measure the masses' heating there.
"Lazy vacuum" hypothesis: R0 = argmin of heating under fixed classicality.

MEASURED (July 2026): NO optimum. While kappa*(R) varies x6, the heating
at threshold is flat to ~3% across R = 0.5..16 in three configurations
(and the residual slope flips sign between them). Reason: an information-
disturbance equality -- for linear monitoring the same coefficient Lambda
gives the record rate AND the momentum diffusion on the mass, and the
threshold pins Lambda* = K/2 (n14). So the source-heating floor is
apparatus-independent, and R0 is a STRUCTURALLY free parameter of the
classical-channel class: only the experimental corridor
(Donadi X-rays from below, Eot-Wash short-range gravity from above)
constrains it. Depends on n14_classical_channel.py. Runtime: ~10 minutes.
"""

import numpy as np
from scipy.linalg import solve_continuous_lyapunov
import n14_classical_channel as m

G_FID = 0.00960
R_SQUEEZE = 2.0
TMAX, NT_STEPS = 600.0, 50

def probes(win, R):
    lo, hi = win
    vs = []
    idx = np.arange(m.N)
    for c0 in range(lo, hi, 2):
        v = np.exp(-0.5 * ((idx - c0) / R) ** 2)
        v[np.abs(idx - c0) > max(4 * R, 2)] = 0.0
        n = np.linalg.norm(v)
        if n > 0:
            vs.append(v / n)
    return vs

def drift_diff_probes(K, kappa, win, R):
    A, D = m.drift_diff(K, 0.0, win)          # bains de bord seuls
    if kappa > 0:
        for v in probes(win, R):
            D[m.NT:m.NT + m.N, m.NT:m.NT + m.N] += kappa * np.outer(v, v)
    return A, D

def run_probes(d, g, kappa, R):
    a, b = m.CENTER - d // 2, m.CENTER - d // 2 + d
    K = m.K_total(a, b, g)
    win = (a - 10, b + 10)
    A, D = drift_diff_probes(K, kappa, win, R)
    siginf = solve_continuous_lyapunov(A, -D)
    sig0 = m.sigma0_product(a, b, g, R_SQUEEZE)
    lam, V, Vinv = m.evolve_setup(A)
    dsig0 = sig0 - siginf
    ts = np.linspace(0, TMAX, NT_STEPS + 1)[1:]
    ENmax, hA_fin = 0.0, 0.0
    for t in ts:
        s = m.sigma_t(t, lam, V, Vinv, dsig0, siginf)
        ENmax = max(ENmax, m.log_negativity(s))
        if t == ts[-1]:
            hA_fin = (s[m.NT + m.N, m.NT + m.N] - 0.5 * m.OMEGA_M
                      + m.OMEGA_M ** 2 * (s[m.N, m.N] - 0.5 / m.OMEGA_M))
    return ENmax, hA_fin

def kstar(d, g, R, eps=1e-3):
    lo, hi = 0.0, 4e-3
    while run_probes(d, g, hi, R)[0] > eps:
        hi *= 2
        if hi > 200: return None
    for _ in range(12):
        mid = 0.5 * (lo + hi)
        if run_probes(d, g, mid, R)[0] > eps: lo = mid
        else: hi = mid
    return 0.5 * (lo + hi)

def lam_record(d, g, kappa, R):
    a, b = m.CENTER - d // 2, m.CENTER - d // 2 + d
    Kf = m.K_total(a, b, 0.0)[:m.N, :m.N]
    imprint = g * np.linalg.solve(Kf, m.cgrad(a))
    return kappa * sum((v @ imprint) ** 2 for v in probes((a - 10, b + 10), R))

def scan(d, g, label):
    a, b = m.CENTER - d // 2, m.CENTER - d // 2 + d
    Keff = abs(m.K_eff(a, b, g))
    _, h0 = run_probes(d, g, 0.0, 1.0)
    print(f"--- {label} : d={d}, g={g:.5f}, |K_eff|={Keff:.3e}, "
          f"chauffe(0)={h0:.4f} ---")
    print(f"{'R':>6} {'kappa*':>10} {'Lambda*/K':>10} {'DH(R)':>12}")
    for R in (0.5, 1.0, 2.0, 4.0, 8.0, 16.0):
        ks = kstar(d, g, R)
        if ks is None:
            print(f"{R:>6} {'> 200':>10}")
            continue
        _, hA = run_probes(d, g, ks, R)
        lam = lam_record(d, g, ks, R)
        print(f"{R:>6} {ks:>10.3e} {lam/Keff:>10.3f} {hA - h0:>12.4e}")
    print()

if __name__ == "__main__":
    scan(12, G_FID, "fiduciel")
    # quelle echelle R0 suit-il ? variation de la portee du champ et de d
    m.MU2 = 1e-2          # mu 0.05 -> 0.1 : portee 20 -> 10
    scan(12, G_FID, "mu x2 (portee /2)")
    m.MU2 = 2.5e-3
    scan(20, G_FID, "d=20")