import json, math
import numpy as np

# Homogeneous cubic field: F_x=alpha*x^3, F_y=(alpha+beta)*x^2*y.
# At u*=e1: a=alpha, G=0, and the tangent Jacobian eigenvalue is beta.
def field(z, alpha, beta):
    x, y = z
    return np.array([alpha*x**3, (alpha+beta)*x*x*y], dtype=float)

def leading_on_unit(u, alpha, beta):
    return field(u, alpha, beta)

def controller(z, alpha, beta, rc=2.0, delta=.15, k=8.0, qtol=.20):
    r = np.linalg.norm(z)
    if r == 0: return field(z, alpha, beta), False
    u = z/r
    v = leading_on_unit(u, alpha, beta)
    a = float(u @ v)
    q = np.linalg.norm(v-a*u)
    # smooth radial damping gated by radius and the near-equilibrium test
    kap = k * np.logaddexp(0., (r-rc)/delta)
    active = (r > rc and q < qtol and a > 0)
    return field(z, alpha, beta) - (kap*u if active else 0.), active

def rk4_step(z, dt, fun):
    k1=fun(z); k2=fun(z+dt*k1/2); k3=fun(z+dt*k2/2); k4=fun(z+dt*k3)
    return z + dt*(k1+2*k2+2*k3+k4)/6

def integrate(z0, T, dt, fun, rcap=1e5):
    z=np.array(z0,dtype=float); maxr=np.linalg.norm(z); interventions=0; first=None
    n=int(math.ceil(T/dt)); actual=0.
    for i in range(n):
        h=min(dt,T-actual)
        if h<=0: break
        znew=rk4_step(z,h,fun)
        actual += h; z=znew
        r=np.linalg.norm(z); maxr=max(maxr,r)
        if getattr(fun,'last_active',False): interventions += 1; first=first if first is not None else actual
        if not np.all(np.isfinite(z)) or r>rcap:
            return {'z':z,'r':r,'maxr':maxr,'t':actual,'failed':True,'interventions':interventions,'first':first}
    return {'z':z,'r':np.linalg.norm(z),'maxr':maxr,'t':actual,'failed':False,'interventions':interventions,'first':first}

def make_control_fun(alpha,beta):
    state={'active':False}
    def f(z):
        v, act=controller(z,alpha,beta)
        state['active']=act
        f.last_active=act
        return v
    f.last_active=False
    return f

def math_check():
    # Exact infinity-atlas quantities at e1, checked by finite differences.
    alpha,beta=1.7,.63; u=np.array([1.,0.]); v=leading_on_unit(u,alpha,beta)
    a=float(u@v); eps=1e-6
    def G(th):
        uu=np.array([math.cos(th),math.sin(th)])
        vv=leading_on_unit(uu,alpha,beta)
        return vv-(uu@vv)*uu
    num=(G(eps)[1]-G(-eps)[1])/(2*eps)
    return {'a_at_u_star':a,'predicted_alpha':alpha,'tangent_eigenvalue_fd':num,'predicted_beta':beta,
            'relative_errors':[abs(a-alpha)/alpha,abs(num-beta)/beta]}

def radial_scaling():
    # Along exact e1, blow-up time is t*=1/(2 alpha r0^2). Measure when r reaches R.
    r0=.8; R=20.; dt=5e-5
    rows=[]
    for alpha in [.5, .8, 1.2, 1.8]:
        target=(1/r0**2-1/R**2)/(2*alpha)
        T=target + .05
        out=integrate([r0,0],T,dt,lambda z:field(z,alpha,0),rcap=R)
        rows.append({'alpha':alpha,'predicted_time_to_R':target,'observed_time':out['t'],'relative_error':abs(out['t']-target)/target})
    return rows

def angular_scaling():
    # For small angle, d log(theta)/d tau = beta. tau=int r^2 dt, and use a short interval.
    r0=.35; theta0=1e-4; T=.25; dt=1e-5
    rows=[]
    for beta in [-.8,.2,.5,1.0,1.5]:
        alpha=0.4
        z0=[r0*math.cos(theta0),r0*math.sin(theta0)]
        out=integrate(z0,T,dt,lambda z:field(z,alpha,beta))
        th0=math.atan2(z0[1],z0[0]); th=math.atan2(out['z'][1],out['z'][0])
        tau=(r0**2)*T # leading approximation, adequate at small T
        measured=math.log(abs(th/th0))/tau
        rows.append({'beta':beta,'predicted_slope':beta,'observed_slope':measured,'relative_error':abs(measured-beta)/max(abs(beta),1e-12)})
    return rows

def controller_sweep():
    # Perturbed near-equilibrium direction: vanilla bursts, controller damps radial growth.
    rows=[]; r0=.8; T=.8; dt=2e-5; beta=.5; theta=1e-3
    for alpha in [.8,1.2,1.6,2.0]:
        z0=[r0*math.cos(theta),r0*math.sin(theta)]
        vanilla=integrate(z0,T,dt,lambda z:field(z,alpha,beta),rcap=1e4)
        cf=make_control_fun(alpha,beta)
        controlled=integrate(z0,T,dt,cf,rcap=1e4)
        rows.append({'alpha':alpha,'vanilla_max_r':vanilla['maxr'],'controlled_max_r':controlled['maxr'],
                     'vanilla_failed':vanilla['failed'],'controlled_failed':controlled['failed'],
                     'controller_interventions':controlled['interventions'],'first_intervention':controlled['first']})
    return rows

def onset_boundary():
    # Along e1, reaching rc by T requires alpha >= (1/r0^2-1/rc^2)/(2T).
    r0=.8; rc=2.; T=.8; dt=5e-5
    predicted=(1/r0**2-1/rc**2)/(2*T)
    rows=[]
    for alpha in [.6,.8,.9,1.0,1.2]:
        out=integrate([r0,0],T,dt,lambda z:field(z,alpha,.5),rcap=1e4)
        rows.append({'alpha':alpha,'reached_rc':bool(out['maxr']>rc),'max_r':out['maxr']})
    return {'predicted_alpha_boundary':predicted,'sweep':rows}

def main():
    result={'math_check':math_check(),'radial_scaling':radial_scaling(),
            'angular_scaling':angular_scaling(),'onset_boundary':onset_boundary(),'controller_sweep':controller_sweep()}
    print(json.dumps(result,indent=2))

if __name__=='__main__': main()