"""n36 — hierarchical (Dyson) record ladder: gates 1-3 pre-registered.
Per-scale channels: rate gamma_l = g0/l^(1+2s); KK forces Omega_l = w0/l^(1+2s);
KMS weight eta_l = eta*l^(s-1) (declared choice: local thermal partner).
Per Fourier mode q: G(q)=sum gamma_l eta_l s_l ; W(q)=sum 2 Omega_l s_l ;
R(q)=sum 2 Omega_l eta_l^2 s_l ;  s_l(q)=4 sin^2(q l/2).
Effective 2x2 per mode:  tau' = -G tau + R eps + lam h ; eps' = -G eps - W tau.
GATE1 p in [1.95,2.05] (chi_eps = W/(G^2+RW) ~ q^-2) for all s.
GATE2 z -> [0.95,1.05] as s -> 0 (omega = sqrt(WR)), underdamped G/omega < 1.
GATE3 drag: moving source, F(v): s=1 must show linear (Demery-Dean) drag;
       s=0.05 must show suppressed drag below c_eff (Cherenkov window)."""
import numpy as np

N=256; G0=0.2; W0=1.0; ETA=0.5; LAM=0.15
ls=np.arange(1,N//2)

def kernels(s):
    gl=G0/ls**(1+2*s); Ol=W0/ls**(1+2*s); el=ETA*ls**(s-1.0)
    qs=2*np.pi*np.arange(N)/N
    sl=4*np.sin(np.outer(qs,ls)/2)**2
    G=sl@(gl*el); W=sl@(2*Ol); R=sl@(2*Ol*el**2)
    return qs,G,W,R

print("GATES 1-2: p and z per s (fit over modes 3..14)")
print(f"{'s':>5} {'p':>8} {'z':>8} {'G/omega':>9} {'c_eff(q~2pi/12)':>16}")
for s in (1.0,0.5,0.25,0.1,0.05):
    qs,G,W,R=kernels(s)
    sel=np.arange(3,15)
    chi=W[sel]/(G[sel]**2+R[sel]*W[sel])
    om=np.sqrt(W[sel]*R[sel])
    p=-np.polyfit(np.log(qs[sel]),np.log(chi),1)[0]
    z= np.polyfit(np.log(qs[sel]),np.log(om),1)[0]
    kf=N//12
    ce=np.sqrt(W[kf]*R[kf])/qs[kf]
    print(f"{s:>5} {p:>8.4f} {z:>8.4f} {(G[sel]/om).mean():>9.3f} {ce:>16.3f}")

print()
print("GATE 3: drag on uniformly moving source (time-averaged self-force)")
def drag(s,v,T=600.0,dt=0.004,W_src=2.0):
    qs,G,W,R=kernels(s)
    # real-space circulant matrices via FFT of kernels
    def mat(K):
        return np.real(np.fft.ifft(np.fft.fft(np.eye(N),axis=0)*K[:,None],axis=0))
    Gm,Wm,Rm=mat(G),mat(W),mat(R)
    tau=np.zeros(N); eps=np.zeros(N)
    x=np.arange(N); F=[]; X=0.0
    steps=int(T/dt)
    for st in range(steps):
        d=(x-X+N/2)%N-N/2
        pr=np.exp(-0.5*(d/W_src)**2); h=pr/pr.sum(); gh=(d/W_src**2)*pr/pr.sum()
        src=LAM*(h-h.mean())
        for _ in range(1):
            dtau=-Gm@tau+Rm@eps+src
            deps=-Gm@eps-Wm@tau
            tau=tau+dt*dtau; eps=eps+dt*deps
        X=(X+v*dt)%N
        if st>steps*0.6 and st%25==0:
            F.append(-LAM*float(gh@(eps-eps.mean())))
    return float(np.mean(F))
for s in (1.0,0.05):
    qs,G,W,R=kernels(s); kf=N//12; ce=np.sqrt(W[kf]*R[kf])/qs[kf]
    print(f"  s={s} (c_eff~{ce:.2f}):")
    for v in (0.05,0.15,0.4,1.2):
        Fd=drag(s,v)
        print(f"    v={v:>5}: <F_self> = {Fd:+.3e}   (v/c_eff = {v/ce:.2f})")

print()
print("wavefront test: pulse at t=0, front position vs t (linear=wave, sqrt=diffusive)")
def front(s,T=120.0,dt=0.004):
    qs,G,W,R=kernels(s)
    def mat(K): return np.real(np.fft.ifft(np.fft.fft(np.eye(N),axis=0)*K[:,None],axis=0))
    Gm,Wm,Rm=mat(G),mat(W),mat(R)
    tau=np.zeros(N); eps=np.zeros(N)
    x=np.arange(N); d=(x-0+N/2)%N-N/2
    pr=np.exp(-0.5*(d/2.0)**2); tau+= pr/pr.sum()   # tau pulse
    out=[]
    for st in range(int(T/dt)):
        dtau=-Gm@tau+Rm@eps; deps=-Gm@eps-Wm@tau
        tau=tau+dt*dtau; eps=eps+dt*deps
        if st%2500==0 and st>0:
            prof=np.abs(eps-eps.mean())
            thr=0.05*prof.max()
            idx=np.where(prof>thr)[0]
            dd=np.minimum(idx,N-idx)
            out.append((st*dt,dd.max()))
    return out
for s in (1.0,0.05):
    o=front(s)
    ts=np.array([a for a,_ in o]); rs=np.array([b for _,b in o])
    al=np.polyfit(np.log(ts),np.log(rs+1e-9),1)[0]
    print(f"  s={s}: front ~ t^{al:.2f}  ({'WAVE' if al>0.8 else 'diffusive' if al<0.65 else 'mixed'})")
print("done.")