# SCRIPT: BH-MIXER-DATA-01
# Measured tables for the black-hole mixer: the Sun's dynamo spectrum on the real axis and its eigenvalue paths around
# loops enclosing the first n birth points (n = 1, 2, 3), from the operator of BH-SUNLADDER-16b; the 8_19 constants
# from BH-KNOT-14; the Ramanujan radial-limit table from BH-SHADOW-07 (as printed in the paper).
import os as _os
_HERE=_os.path.dirname(_os.path.abspath(__file__))
_FIG=_os.path.join(_HERE,'..','paper','fig') if _os.path.isdir(_os.path.join(_HERE,'..','paper','fig')) else _os.path.join(_HERE,'..','fig')
import numpy as np, json, scipy.sparse as sp, scipy.sparse.linalg as sla
N=120; h=np.pi/(N+1); x=np.arange(1,N+1)*h
D2=sp.diags([np.ones(N-1),-2*np.ones(N),np.ones(N-1)],[-1,0,1])/h**2
D1=sp.diags([-np.ones(N-1),np.ones(N-1)],[-1,1])/(2*h); Cc=sp.diags(np.cos(x)); K=12
def spec(a):
    A=sp.bmat([[D2,a*Cc],[D1,D2]]).tocsc().astype(complex)
    w=sla.eigs(A,k=K+4,sigma=0,which='LM',return_eigenvectors=False); return w[np.argsort(-w.real)][:K]
births=[3.8717,5.2802,7.1576,9.0850,11.0266]
al=np.linspace(2.5,12.0,240); W=np.array([spec(a) for a in al])[:,:11]
def loop(n,ht=0.2):
    base=3.0; a1=(births[n-1]+births[n])/2
    L1=a1-base; Lc=np.pi*ht; T=2*ht+2*L1+Lc
    def P(s):
        if s<=ht: return base-1j*s
        s-=ht
        if s<=L1: return base+s-1j*ht
        s-=L1
        if s<=Lc: return a1+ht*np.exp(1j*(-np.pi/2+s/ht))
        s-=Lc
        if s<=L1: return a1-s+1j*ht
        s-=L1
        return base+1j*(ht-s)
    s=0; cur=spec(P(0)); start=cur.copy(); tr=[cur.copy()]; ps=[P(0)]; ds=0.01
    while s<T-1e-12:
        step=min(ds,T-s)
        while True:
            new=spec(P(s+step)); Dm=np.abs(cur[:,None]-new[None,:]); asg=np.argmin(Dm,axis=1)
            good=len(set(asg))==K and all(np.sort(Dm[j])[0]<=0.25*np.sort(Dm[j])[1] for j in range(K))
            if good or step<1e-6: break
            step/=2
        cur=new[asg]; s+=step; tr.append(cur.copy()); ps.append(P(s)); ds=min(0.01,step*1.5)
    tr=np.array(tr); ps=np.array(ps)
    idx=np.linspace(0,len(tr)-1,420).astype(int)
    perm=[int(np.argmin(np.abs(start-c))) for c in cur]
    return {'alpha':[[round(z.real,4),round(z.imag,4)] for z in ps[idx]],
            'lam':[[[round(z.real,4),round(z.imag,4)] for z in tr[i,:9]] for i in idx],'perm':perm[:9]}
data={'sun':{'alpha':[round(a,4) for a in al],'re':np.round(W.real,4).tolist(),'im':np.round(W.imag,4).tolist(),'births':births,
       'loops':{str(n):loop(n) for n in (1,2,3)}},
      'k819':json.load(open(_os.path.join(_HERE,'bh_knot_14.json')))['c640'],
      'shadow':[['-1','1e-3',180.26,'3.996012','4'],['-1','3e-4',597.24,'3.998801','4'],['i','3e-4',150.67,'0.002390 + 3.994025 i','4i'],['e^{iπ/3}','3e-4',67.91,'−1.988709 − 3.450721 i','−2 − 2√3 i']]}
for n in (1,2,3): print(n,data['sun']['loops'][str(n)]['perm'])
json.dump(data,open('mixer_data.json','w'),separators=(',',':'))