# SCRIPT: RS02C_REALIZATION
# Are the steps of RS-02 realised by actual projector families (the IQGT object)?
# (1) d = (sin b cos kx, sin b sin kx, cos b), b = beta(y): a = k^2 sin^2 beta / 4, b = beta'^2/4.
# (2) Sigma = cos3x sx + sin3x sy + m sz: a * gap^2 = 9, so n = sqrt(a) = 3/gap, and the
#     normal-incidence Gamma between two gaps is tanh(0.5 ln(gap2/gap1)), the peel form.
import numpy as np
sx=np.array([[0,1],[1,0]],complex); sy=np.array([[0,-1j],[1j,0]]); sz=np.diag([1.,-1.]).astype(complex)
def lowP(H):
    w,v=np.linalg.eigh(H); return np.outer(v[:,0],v[:,0].conj())
def metric(Hf,x,y,h=1e-5):
    dx=(lowP(Hf(x+h,y))-lowP(Hf(x-h,y)))/(2*h); dy=(lowP(Hf(x,y+h))-lowP(Hf(x,y-h)))/(2*h)
    g=lambda A,B: 0.5*np.trace(A@B).real
    return g(dx,dx),g(dy,dy),g(dx,dy)
rng=np.random.default_rng(20261020); worst1=0
for _ in range(50):
    k=rng.uniform(1,4); b0=rng.uniform(0.2,2.9); b1=rng.uniform(-3,3)
    beta=lambda y: b0+b1*y
    Hf=lambda x,y: np.sin(beta(y))*np.cos(k*x)*sx+np.sin(beta(y))*np.sin(k*x)*sy+np.cos(beta(y))*sz
    gxx,gyy,gxy=metric(Hf,0.3,0.0)
    worst1=max(worst1,abs(gxx-k*k*np.sin(b0)**2/4),abs(gyy-b1*b1/4),abs(gxy))
print(f"(1) two-level family: max error in a=k^2 sin^2(beta)/4, b=beta'^2/4, g_xy=0 over 50 draws: {worst1:.1e}")
worst2=0; worst3=0
for m1,m2 in [(0.,1.),(0.5,2.),(-1.,3.),(0.2,0.25)]:
    Hf=lambda m: (lambda x,y: np.cos(3*x)*sx+np.sin(3*x)*sy+m*sz)
    a1=metric(Hf(m1),0.4,0.0)[0]; a2=metric(Hf(m2),0.4,0.0)[0]
    g1=2*np.sqrt(1+m1*m1); g2=2*np.sqrt(1+m2*m2)
    worst2=max(worst2,abs(a1*g1*g1-9),abs(a2*g2*g2-9))
    G=(np.sqrt(a1)-np.sqrt(a2))/(np.sqrt(a1)+np.sqrt(a2)); peel=np.tanh(0.5*np.log(g2/g1))
    worst3=max(worst3,abs(G-peel))
    print(f"   m {m1:+.2f} -> {m2:+.2f}: gaps {g1:.4f} -> {g2:.4f}; Gamma(n=sqrt g_par) = {G:+.8f}; tanh(0.5 ln(gap2/gap1)) = {peel:+.8f}")
print(f"(2) max |a*gap^2 - 9| = {worst2:.1e}; max |Gamma - peel| = {worst3:.1e}")
print("Realization:", "PASS" if max(worst1,worst2,worst3)<1e-6 else "KILL")