"""n48 — motion-gated decoherence, per n48_preregistration.md (+ Amendment 1).
Overdamped drainage chasing a modulated mass superposition; decoherence
rate = time-averaged squared violation difference between branches."""
import numpy as np

N = 200; MU2 = 1e-6; F = 0.3; EPS = 0.5

def K_of(dk):
    k = 1.0 + dk
    K = MU2*np.eye(N)
    for e in range(N):
        i, j = e, (e+1) % N
        K[i,i] += k[e]; K[j,j] += k[e]; K[i,j] -= k[e]; K[j,i] -= k[e]
    return K

# info-metric drainage map dk*(c) (n46 machinery)
w2, V = np.linalg.eigh(K_of(np.zeros(N)))
wv = np.sqrt(np.maximum(w2, 1e-30))
U = np.zeros((N, N))
for e in range(N): U[e] = V[e] - V[(e+1) % N]
W1, W2 = np.meshgrid(wv, wv, indexing='ij')
kern = 1.0/((W1+W2)**2 * W1*W2)
B = np.einsum('en,em->enm', U, U).reshape(N, N*N)
G = 0.125*(B*kern.reshape(1, N*N)) @ B.T
A = np.zeros((N, N))
for i in range(N): A[i,(i-1)%N] = 1; A[i,i] += 1
Ginv_At = np.linalg.solve(G, A.T)
Wm = A @ Ginv_At
Wp = np.linalg.pinv(Wm, rcond=1e-10)
Dmap = -Ginv_At @ Wp            # dk* = Dmap @ c

def pattern(sites):
    c = np.zeros(N)
    for s in sites: c[s] = 1.0
    return c

cL, cR = pattern([60, 161]), pattern([62, 163])
dkL0, dkR0 = Dmap @ cL, Dmap @ cR
# violation of the lagging field: v = A dk + c*f(t); with A dk* = -c (neutral),
# v(t) = A (dk - f(t) dk*0). Integrate branch lag x = dk - f(t) dk*0 exactly:
# d(dk)/dt = -gamma (dk - f(t) dk*0)  (linear, scalar per mode along dk*0)
def Gamma_dec(nu, gamma):
    """time-averaged |Delta v|^2 in steady state (numerical integration)."""
    if nu == 0.0:
        # relax to static target from vacuum, measure residual after 5/gamma... integrate
        T = 5.0/gamma; dt = 0.002/gamma
        aL = 0.0; aR = 0.0   # amplitude of dk along dk*0 (starts undrained)
        acc = 0.0; steps = int(T/dt)
        for i in range(steps):
            aL += dt*(-gamma*(aL - F)); aR += dt*(-gamma*(aR - F))
        dv = A @ (aL*dkL0 - aR*dkR0) + F*(cL - cR) - (A @ (F*dkL0) + F*cL - A@(F*dkR0) - F*cR)
        # measure over one further unit of time
        vL = A@(aL*dkL0) + F*cL; vR = A@(aR*dkR0) + F*cR
        return np.mean((vL - vR)**2) * N   # site-summed
    Tper = 1.0/nu
    dt = min(0.002/gamma, Tper/2000)
    # transient 5/gamma + 2 periods, then average 3 periods
    t = 0.0; aL = F; aR = F      # start on the static solution
    acc = 0.0; navg = 0
    T1 = 5.0/gamma + 2*Tper; T2 = T1 + 3*Tper
    while t < T2:
        f = F*(1 + EPS*np.sin(2*np.pi*nu*t))
        aL += dt*(-gamma*(aL - f)); aR += dt*(-gamma*(aR - f))
        t += dt
        if t >= T1:
            vL = A@(aL*dkL0) + f*cL
            vR = A@(aR*dkR0) + f*cR
            acc += np.sum((vL - vR)**2); navg += 1
    return acc/navg

gamma = 1.0
nus = np.array([0.0] + list(10.0**np.linspace(-2, 2, 17)))
Gs = np.array([Gamma_dec(nu, gamma) for nu in nus])
Gref = Gs[np.argmin(np.abs(nus - gamma))]
print("n48: Gamma_dec(nu), gamma = 1")
for nu, g in zip(nus, Gs):
    print(f"  nu/gamma = {nu:8.3f}: Gamma = {g:.6e}")
p1 = Gs[0] < 1e-12*Gref
print(f"P1 no static floor: Gamma(0)/Gamma(gamma) = {Gs[0]/Gref:.2e} -> {'PASS' if p1 else 'FAIL'}")
lo = (nus >= gamma/100) & (nus <= gamma/10)
slope_lo = np.polyfit(np.log(nus[lo]), np.log(Gs[lo]), 1)[0]
print(f"P2 quadratic onset: slope = {slope_lo:.4f} -> {'PASS' if abs(slope_lo-2) < 0.05 else 'FAIL'}")
hi = (nus >= 10*gamma)
slope_hi = np.polyfit(np.log(nus[hi]), np.log(Gs[hi]), 1)[0]
print(f"P3 saturation: high-nu slope = {slope_hi:.4f} -> {'PASS' if abs(slope_hi) < 0.2 else 'FAIL'}")
# knee: max derivative of log Gamma vs log nu transition midpoint
dlog = np.diff(np.log(Gs[1:]))/np.diff(np.log(nus[1:]))
mid = nus[1:][np.argmin(np.abs(dlog - 1.0))]
print(f"P3 knee location ~ nu = {mid:.3f} gamma -> {'PASS' if gamma/3 <= mid <= 3*gamma else 'FAIL'}")
# P4: gamma scaling at nu = 0.01
g4 = [Gamma_dec(0.01, g) for g in (1.0, 2.0, 4.0)]
expo = np.polyfit(np.log([1.0, 2.0, 4.0]), np.log(g4), 1)[0]
print(f"P4 gamma exponent at nu << gamma: {expo:.3f} (page said -1; registered: measure truth) "
      f"-> {'within gate' if abs(expo + 2) < 0.1 else ('CORRECTS THE PAGE' if abs(expo+2) < 0.3 else 'unexpected')}")