# SCRIPT: BH-SUNLADDER-16
# One loop around the first n exceptional points of the Sun's dynamo operator: does it climb the ladder
# (as the Kerr cascade does) or is it a product of separate swaps?
# Operator (the solar paper's minimal Parker alpha-Omega operator): on x in [0, pi], Dirichlet,
#   dA/dt = alpha0 cos(x) B + A'',   dB/dt = G A' + B'',   G = eta = 1,
# second-order finite differences with N = 120 interior points: the same matrix as the solar paper's
# SUN-OPEN-02 (suns_v1_5_complete.zip), rebuilt here line for line.
# Why this decides the question: on a closed loop in the complex alpha0 plane the eigenvalues return as a
# permutation of themselves, and that permutation is fixed by which branch points the loop encloses. If
# neighbouring exceptional points share a branch (1-2 at the first, 2-3 at the second, ...), a loop around
# the first n points is a single cycle that carries the starting eigenvalue up the ladder; if each point
# belongs to its own pair, the loop is a product of disjoint swaps and every eigenvalue moves at most one
# place. Loops are thin stadiums hugging the real axis, at two heights, so the answer is checked not to
# depend on enclosing complex-alpha0 branch points by accident.
# Kill conditions, written before the run:
#  K1  control: the scan reproduces the solar paper's ladder 3.872, 5.280, 7.158, 9.085, 11.032 to 2e-3,
#      as points where the number of real eigenvalues among the top 14 drops by two. Kill otherwise.
#  K2  each single-point loop exchanges exactly one pair (one transposition), as in the solar paper. Kill otherwise.
#  K3  the permutation after a loop around the first n points, n = 1..5, is the same at both stadium heights.
#      Kill (no verdict) otherwise.
#  Verdict (not a kill): "ladder" if the permutation for n points contains an (n+1)-cycle through the
#  fundamental; "separate swaps" if it is a product of n disjoint transpositions; otherwise report it as found.
import numpy as np
N=120; h=np.pi/(N+1); x=np.arange(1,N+1)*h
D2=(np.diag(-2*np.ones(N))+np.diag(np.ones(N-1),1)+np.diag(np.ones(N-1),-1))/h**2
D1=(np.diag(np.ones(N-1),1)-np.diag(np.ones(N-1),-1))/(2*h)
C=np.diag(np.cos(x)); Z=np.zeros((N,N)); I=np.eye(N)
def L(a): return np.block([[D2,a*C],[D1,D2]])
K=14
def top(a):
    ev=np.linalg.eigvals(L(a)); return ev[np.argsort(-ev.real)][:K]
ok=True
def report(name,cond,val=''):
    global ok
    print(('PASS ' if cond else 'FAIL ')+name, val); ok=ok and cond
# K1: scan for pair births (real count drops by 2), refine by bisection
def nreal(a,tol=1e-7):
    ev=top(a); return int(np.sum(np.abs(ev.imag)<tol*np.maximum(1,np.abs(ev))))
grid=np.linspace(0.5,13,2501); counts=[nreal(a) for a in grid]
births=[]
for i in range(len(grid)-1):
    if counts[i+1]<counts[i]:
        lo,hi=grid[i],grid[i+1]
        for _ in range(40):
            mid=(lo+hi)/2
            if nreal(mid)==counts[i]: lo=mid
            else: hi=mid
        births.append((lo+hi)/2)
ref=[3.872,5.280,7.158,9.085,11.032]
print('   pair births found:', [round(b,4) for b in births])
report('K1 ladder reproduced', len(births)==5 and all(abs(b-r)<2e-3 for b,r in zip(births,ref)), '')
# loop tracking
def track(path):
    cur=top(path[0]); start=cur.copy()
    for a in path[1:]:
        new=top(a); used=set(); nxt=np.empty_like(cur)
        for j,c in enumerate(cur):
            d=np.abs(new-c); d[list(used)]=np.inf; k=int(np.argmin(d)); used.add(k); nxt[j]=new[k]
        cur=nxt
    # permutation: where did each starting eigenvalue land?
    perm=[int(np.argmin(np.abs(start-c))) for c in cur]
    err=max(np.min(np.abs(start-c)) for c in cur)
    return perm,err
def stadium(a0,a1,ht,ds=0.008):
    # counterclockwise: along the lower edge a0->a1, around the right cap, upper edge a1->a0, left cap;
    # steps of arc length ds (ds/ht <= 0.16, so the tracking step is small next to the distance to each branch point)
    M=int(np.ceil((a1-a0)/ds))+1; Mc=int(np.ceil(np.pi*ht/ds))+8
    s=np.linspace(0,1,M); pts=[]
    pts+=list(a0+(a1-a0)*s[:-1]-1j*ht)
    th=np.linspace(-np.pi/2,np.pi/2,Mc)[:-1]; pts+=list(a1+ht*np.exp(1j*th))
    pts+=list(a1+(a0-a1)*s[:-1]+1j*ht)
    th=np.linspace(np.pi/2,3*np.pi/2,Mc); pts+=list(a0+ht*np.exp(1j*th))
    return np.array(pts)
def cycles(perm):
    seen=set(); out=[]
    for i in range(len(perm)):
        if i in seen: continue
        c=[i]; seen.add(i); j=perm[i]
        while j!=i: c.append(j); seen.add(j); j=perm[j]
        if len(c)>1: out.append(c)
    return out
a_left=births[0]-0.4
res={}
for ht in (0.05,0.2):
    for n in range(1,6):
        a_right=births[n-1]+0.4 if n<5 else births[4]+0.4
        if n<5: a_right=min(a_right,(births[n-1]+births[n])/2)
        perm,err=track(stadium(a_left,a_right,ht))
        res[(ht,n)]=cycles(perm)
        print(f'   height {ht}: loop around first {n} point(s) [{a_left:.3f}, {a_right:.3f}]: cycles {cycles(perm)}  (return error {err:.1e})')
single=[]
for n in range(1,6):
    lo=births[n-1]-0.25 if n==1 else (births[n-2]+births[n-1])/2
    hi=births[n-1]+0.25 if n==5 else (births[n-1]+births[n])/2
    perm,err=track(stadium(lo,hi,0.05)); single.append(cycles(perm))
    print(f'   single loop around point {n}: cycles {cycles(perm)}')
report('K2 each single loop is one transposition', all(len(c)==1 and len(c[0])==2 for c in single), '')
report('K3 loop permutations agree at both heights', all(res[(0.05,n)]==res[(0.2,n)] for n in range(1,6)), '')
for n in range(1,6):
    cs=res[(0.05,n)]
    lad=any(len(c)==n+1 and 0 in c for c in cs)
    sw=len(cs)==n and all(len(c)==2 for c in cs)
    print(f'   verdict n={n}: '+('LADDER (one (n+1)-cycle through the fundamental)' if lad else 'SEPARATE SWAPS' if sw else 'OTHER: '+str(cs)))
print('ALL PASS' if ok else 'SOME FAILED')