# SCRIPT: SINE_GUARD_CHECK
# At a point of a projector family with Killing direction x, take a unit path direction V (in g).
# Claim (derived): Gram determinant of (V, d_x) = g_par - C, with C = g_par sin^2(theta) the Clairaut
# constant, so the curvature ceiling 2 sqrt(det) = 2 sqrt(g_par) cos(theta) = 2 Y/omega, the wave admittance.
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 P(x,y):  # IQGT 8.6 family, gap narrows on y=0
    w,v=np.linalg.eigh(np.cos(3*x)*sx+np.sin(3*x)*sy+2.5*y*sz); return np.outer(v[:,0],v[:,0].conj())
h=1e-5; rng=np.random.default_rng(1); worst=0; worstO=0
for _ in range(200):
    x,y=rng.uniform(-1,1,2)
    dx=(P(x+h,y)-P(x-h,y))/(2*h); dy=(P(x,y+h)-P(x,y-h))/(2*h); Pi=P(x,y)
    g=lambda A,B: 0.5*np.trace(A@B).real
    Om=lambda A,B: (1j*np.trace(Pi@(A@B-B@A))).real
    gxx,gyy,gxy=g(dx,dx),g(dy,dy),g(dx,dy)
    th=rng.uniform(-1.5,1.5)                       # angle from the interface normal, in g
    V=(np.sin(th)/np.sqrt(gxx))*dx+(np.cos(th)/np.sqrt(gyy))*dy   # unit in g (gxy=0 here)
    C=gxx*np.sin(th)**2
    det=g(V,V)*gxx-g(V,dx)**2
    worst=max(worst,abs(det-(gxx-C)), abs(gxy))
    worstO=max(worstO, abs(Om(V,dx))-2*np.sqrt(max(det,0)))
print(f"max |det - (g_par - C)| = {worst:.1e};  max (|Omega(V,dx)| - 2 sqrt(det)) = {worstO:.1e} (must be <= 0)")