# SCRIPT: RS01B_LEAD_DIAGNOSIS
# Diagnosis of the RS-01 K2 kill. Hypothesis: the declared lead joins the curved metric with a
# kink (a' jumps at y_L), which reflects a small O(1/omega) echo that beats with the main
# reflection and makes the lateral shift oscillate. Test: replace the kink by a C-infinity blend
# of width ell. If the hypothesis is right the error becomes small and falls with omega;
# if the oscillation survives the blend, the hypothesis is wrong.
import json, numpy as np
from scipy.integrate import solve_ivp, quad
from scipy.optimize import brentq
K=3.0; YL=1.2; B=0.25
def S(t):
    t=np.clip(t,0,1)
    if t<=0: return 0.0
    if t>=1: return 1.0
    e1=np.exp(-1/t); e2=np.exp(-1/(1-t)); return e1/(e1+e2)
def dS(t,h=1e-7): return (S(t+h)-S(t-h))/(2*h)
def make(ell):
    YE=YL+ell
    aex=lambda y: 0.25*K*K*np.sin(y)**2
    daex=lambda y: 0.25*K*K*np.sin(2*y)
    aE=aex(YE) if ell>0 else aex(YL)
    def a(y):
        if y<=YL: return aex(y)
        if ell==0 or y>=YE: return aE
        t=(y-YL)/ell; s=S(t); return aex(y)*(1-s)+aE*s
    def da(y):
        if y<=YL: return daex(y)
        if ell==0 or y>=YE: return 0.0
        t=(y-YL)/ell; s=S(t); return daex(y)*(1-s)+(aE-aex(y))*dS(t)/ell
    return a,da,(YE if ell>0 else YL),aE
def run(ell, omegas, thetas):
    a,da,YE,aE=make(ell)
    Q=lambda y,w,p: B*(w*w-p*p/a(y))
    def rhs(y,s,w,p): return [s[1], -da(y)/(2*a(y))*s[1]-Q(y,w,p)*s[0]]
    def refl(w,p):
        yt=np.arcsin(np.sqrt(4*p*p/(w*w)/(K*K)))
        f=lambda y0: quad(lambda y: np.sqrt(max(-Q(y,w,p),0)),y0,yt,limit=200)[0]-30
        y0=1e-6 if f(1e-6)<0 else brentq(f,1e-6,yt*(1-1e-9))
        sol=solve_ivp(rhs,[y0,YE],[1.0,np.sqrt(-Q(y0,w,p))],args=(w,p),method='DOP853',rtol=1e-12,atol=1e-30,max_step=0.01)
        u,du=sol.y[0,-1],sol.y[1,-1]; kE=np.sqrt(B*(w*w-p*p/aE))
        return 0.5*(u+du/(1j*kE))/(0.5*(u-du/(1j*kE)))
    def geo(C):
        yt=np.arcsin(np.sqrt(4*C/(K*K)))
        g=lambda s: 2*s*np.sqrt(B*C)/(np.sqrt(a(yt+s*s))*np.sqrt(a(yt+s*s)-C))
        return 2*quad(g,0,np.sqrt(YE-yt),limit=600,epsabs=1e-13,epsrel=1e-12)[0]
    rows=[]
    for th in thetas:
        C=aE*np.sin(np.radians(th))**2; dx=geo(C)
        for w in omegas:
            p=w*np.sqrt(C); dp=p*1e-5
            D=-np.angle(refl(w,p+dp)/refl(w,p-dp))/(2*dp)
            rows.append((th,w,D,dx,(D-dx)/dx))
    return rows
res={}
omegas=[20,40,80,160,320,640]; thetas=[20,35,50,65]
for ell in [0.0,0.1,0.3]:
    rows=run(ell,omegas,thetas); res[str(ell)]=rows
    print(f"--- blend width ell = {ell}")
    for th in thetas:
        r=[x for x in rows if x[0]==th]
        errs=' '.join(f"{x[4]:+.2e}" for x in r)
        sl=np.polyfit(np.log([x[1] for x in r]),np.log([abs(x[4]) for x in r]),1)[0]
        print(f"theta {th:2d}: rel_err over omega {omegas}: {errs}  slope {sl:+.2f}")
json.dump(res,open('rs01b_results.json','w'),indent=1,default=float)