"""n34 — dynamic matter against the FULL n33 model (linear => full model IS the judge).
Sections: B chi check | C two falling bodies + N3(t) | G lag & drag | H self-force
controls Delta=0, sign flip | E EP dynamic | F force noise | I entanglement (BMV)."""
import numpy as np
from scipy.linalg import solve_continuous_lyapunov as lyap

ETA=0.5
def build_field(N,g,Delta,kappa):
    B=np.zeros((N,N))
    for e in range(N): B[e,e]=1; B[e,(e+1)%N]=-1
    E=N; n=2*N+2*E; A=np.zeros((n,n)); sq=np.sqrt(2)*g
    A[0:N,2*N+E:]+= sq*ETA*B.T
    A[N:2*N,2*N:2*N+E]-= sq*B.T
    A[2*N:2*N+E,2*N+E:]+=Delta*np.eye(E); A[2*N:2*N+E,N:2*N]+=sq*ETA*B
    A[2*N:2*N+E,2*N:2*N+E]-=kappa/2*np.eye(E)
    A[2*N+E:,2*N:2*N+E]-=Delta*np.eye(E); A[2*N+E:,0:N]-=sq*B
    A[2*N+E:,2*N+E:]-=kappa/2*np.eye(E)
    return A,B

N=48; KAP=4.0; G=0.6; LAM=0.15; E0=1.0; W=1.5; M=500.0
gam=G*G*KAP/(4+KAP**2/4); Om0=G*G*2/(4+KAP**2/4)   # |Delta|=2
GamD=gam*ETA
print(f"params: gamma_eff={gam:.4f} Omega(|D|=2)={Om0:.4f} Gamma={GamD:.4f}")

def prof(X):
    x=np.arange(N); d=(x-X+N/2)%N-N/2
    p=np.exp(-0.5*(d/W)**2); return E0*p/p.sum(), E0*(d/W**2)*p/p.sum()

def run_two_body(Delta,T=3000,dt=0.05,XA0=18.0,XB0=30.0,PA0=0.0,drag_test=False):
    A,B=build_field(N,G,Delta,KAP); n=A.shape[0]
    R=np.zeros(n); XA,PA,XB,PB=XA0,PA0,XB0,0.0
    out=[]
    for s in range(int(T/dt)):
        hA,dA=prof(XA); hB,dB=prof(XB)
        b=np.zeros(n); src=LAM*(hA+hB); b[0:N]=src-src.mean()
        eps=R[N:2*N]-R[N:2*N].mean()
        FA=-LAM*float(dA@eps)*(-1); FB=-LAM*float(dB@eps)*(-1)
        # F = -d<H_int>/dX = -lam * sum_i eps_i dh/dX ; dh/dX = -dprof/di => sign
        FA=LAM*float(dA@eps); FB=LAM*float(dB@eps)
        k1=A@R+b; k2=A@(R+dt/2*k1)+b; k3=A@(R+dt/2*k2)+b; k4=A@(R+dt*k3)+b
        R=R+dt/6*(k1+2*k2+2*k3+k4)
        XA+=dt*PA/M; PA+=dt*FA
        if not drag_test: XB+=dt*PB/M; PB+=dt*FB
        if s%200==0: out.append((s*dt,XA,XB,PA,PB,FA,FB))
    return np.array(out)

print("\n== C: two bodies released at rest, Delta=-2 (predicted attractive) ==")
o=run_two_body(-2.0)
t=o[:,0]; XA=o[:,1]; XB=o[:,2]; FA=o[:,5]; FB=o[:,6]
print(f"  separation: {XB[0]-XA[0]:.2f} -> {XB[-1]-XA[-1]:.2f}  ({'FALLING' if XB[-1]-XA[-1]<XB[0]-XA[0] else 'repelling'})")
mask=t>300
n3=np.abs(FA+FB)/(0.5*(np.abs(FA)+np.abs(FB))+1e-30)
print(f"  N3(t): max(after buildup)={n3[mask].max():.2%}  mean={n3[mask].mean():.2%}  (transient max={n3.max():.2%})")
print(f"  P_A+P_B drift={o[-1,3]+o[-1,4]:+.3e} vs |P_A|={abs(o[-1,3]):.3e} -> momentum into field during transient")

print("\n== controls: Delta=0 and Delta=+2, same setup ==")
o0=run_two_body(0.0,T=1500); op=run_two_body(+2.0,T=1500)
F0=np.abs(o0[:,5]).max(); Fm=np.abs(o[:,5]).max()
print(f"  |F|max Delta=0: {F0:.3e} vs Delta=-2: {Fm:.3e}   ratio={F0/Fm:.2e} (gate <1e-2)")
dp=op[-1,2]-op[-1,1]; dm=o[int(len(op))-1,2]-o[int(len(op))-1,1]
print(f"  Delta=+2: separation {op[0,2]-op[0,1]:.2f} -> {dp:.2f}  ({'REPELLING' if dp>op[0,2]-op[0,1] else 'attract'}) — dynamic sign flip")

print("\n== G: drag on a single moving body (self-force of the wake) ==")
for D,tag in ((-2.0,"Delta=-2"),(0.0,"Delta=0")):
    od=run_two_body(D,T=1200,XA0=10.0,PA0=M*0.02,XB0=1e9%N,drag_test=True)  # B parked far & frozen
    # B frozen at some site still sources; instead set E0 of B by using only A: quick hack: measure dP_A over window
    dP=od[-1,3]-od[int(len(od)*0.3),3]; dT=od[-1,0]-od[int(len(od)*0.3),0]
    print(f"  {tag}: <dP_A/dt> = {dP/dT:+.3e} at v=0.02  (drag coefficient ~ {dP/dT/0.02:+.3e})")

print("\n== E: dynamic EP — two probes, same E same M, widths 1.2 vs 2.4, same external field ==")
Af,_=build_field(N,G,-2.0,KAP)
hS,_=prof(6.0)
b=np.zeros(Af.shape[0]); b[0:N]=LAM*3.0*(hS-hS.mean())   # heavy source E=3
Rss=np.linalg.lstsq(Af,-b,rcond=None)[0]; eps_ext=Rss[N:2*N]-Rss[N:2*N].mean()
def prof_w(X,wd):
    x=np.arange(N); d=(x-X+N/2)%N-N/2
    p=np.exp(-0.5*(d/wd)**2); return E0*(d/wd**2)*p/p.sum()
FF=[LAM*float(prof_w(24.0,wd)@eps_ext) for wd in (1.2,2.4)]
print(f"  F(w=1.2)={FF[0]:+.5e}  F(w=2.4)={FF[1]:+.5e}  Delta_EP={abs(FF[0]-FF[1])/(0.5*(abs(FF[0])+abs(FF[1]))):.2%}")
print("  NOTE (analytic): force per unit E universal; ACCELERATION-EP additionally requires E prop M_inertial,")
print("  a relation absent from nonrelativistic toy by construction (the 1/c^2 gap).")

print("\n== F: force noise (Lyapunov NESS covariance, N=24) ==")
N2=24
A2,_=build_field(N2,G,-2.0,KAP)
D2=np.zeros_like(A2); D2[2*N2:,2*N2:]=KAP/2*np.eye(2*N2)   # vacuum noise on aux quadratures
S=lyap(A2,-D2)
See=S[N2:2*N2,N2:2*N2]
x=np.arange(N2); d=(x-12.0+N2/2)%N2-N2/2; p=np.exp(-0.5*(d/W)**2); dprof=E0*(d/W**2)*p/p.sum()
varF=LAM**2*float(dprof@See@dprof)
# mean pair force at d=6 from unit source (static formula from NESS solve)
hB2=np.exp(-0.5*(((x-6.0+N2/2)%N2-N2/2)/W)**2); hB2/=hB2.sum()
b2=np.zeros(A2.shape[0]); b2[0:N2]=LAM*(hB2-hB2.mean())
e2=np.linalg.lstsq(A2,-b2,rcond=None)[0][N2:2*N2]; e2-=e2.mean()
Fmean=LAM*float(dprof@e2)
print(f"  sqrt(Var F)={np.sqrt(varF):.3e}   |<F>|(d=6, E=1)={abs(Fmean):.3e}   ratio={np.sqrt(varF)/abs(Fmean):.1f}")
print(f"  ratio scales as 1/(lam*E_source): classical limit needs lam*E >> vacuum scale")

print("\n== I: entanglement through the mediator (Gaussian BMV), transient E_N(t) ==")
def bmv(Delta,T=400,dt=0.02):
    Nn=16; A,B=build_field(Nn,G,Delta,KAP); nf=A.shape[0]
    m=1.0; om=0.3; c=0.25   # two matter oscillators, dipole-coupled to eps at sites 4 and 12
    n=nf+4; AA=np.zeros((n,n)); AA[:nf,:nf]=A
    dvec={}
    for k,site in enumerate((4,12)):
        d=(np.arange(Nn)-site+Nn/2)%Nn-Nn/2; pr=np.exp(-0.5*(d/1.5)**2); pr/=pr.sum()
        dvec[k]=c*(d/1.5**2)*pr
    for k in (0,1):
        ix,ip=nf+2*k,nf+2*k+1
        AA[ix,ip]=1/m; AA[ip,ix]=-m*om*om
        AA[ip,Nn:2*Nn]-=dvec[k]                      # p_X reads eps
        AA[0:Nn,ix]+=dvec[k]                          # x_X writes into tau (dH/deps struct? via eps-coupling)
    DD=np.zeros((n,n)); DD[2*Nn:2*Nn+2*Nn,2*Nn:2*Nn+2*Nn]=KAP/2*np.eye(2*Nn)
    S=np.zeros((n,n))
    # initial: matter vacuum, clock/aux small vacuum-ish
    for k in (0,1):
        S[nf+2*k,nf+2*k]=1/(2*m*om); S[nf+2*k+1,nf+2*k+1]=m*om/2
    S[2*Nn:2*Nn+2*Nn,2*Nn:2*Nn+2*Nn]+=0.5*np.eye(2*Nn)
    S[0:Nn,0:Nn]+=0.5*np.eye(Nn); S[Nn:2*Nn,Nn:2*Nn]+=0.5*np.eye(Nn)
    ENs=[]
    for s in range(int(T/dt)):
        dS=AA@S+S@AA.T+DD
        S=S+dt*dS
        if s%400==0:
            sm=np.zeros((4,4))
            idx=[nf,nf+1,nf+2,nf+3]
            for a in range(4):
                for bq in range(4): sm[a,bq]=S[idx[a],idx[bq]]
            Aa=sm[:2,:2]; Bb=sm[2:,2:]; Cc=sm[:2,2:]
            Dt=np.linalg.det(Aa)+np.linalg.det(Bb)-2*np.linalg.det(Cc)
            disc=Dt*Dt-4*np.linalg.det(sm)
            if disc<0: ENs.append(0.0); continue
            nu2=(Dt-np.sqrt(disc))/2
            ENs.append(max(0.0,-np.log(2*np.sqrt(max(nu2,1e-30)))))
    return max(ENs)
for D in (-0.25,-1.0,-4.0):
    print(f"  Delta={D:+.2f} (|D|/kappa={abs(D)/KAP:.2f}): max E_N = {bmv(D):.4f}")
print("done.")