import json, numpy as np
from vdp_cell import simulate, radius

def rate(z,mu,R):
    n=len(z)//2; x=z[:n]; y=z[n:]
    return mu*(1-np.dot(z,z)/(R*R))*np.dot(y,y)

def fit(x,y):
    return float(np.dot(x-x.mean(),y-y.mean())/np.dot(x-x.mean(),x-x.mean()))
R=1.; omega=2.; mu=1.
signs=[]
for r in [.5,.8,1.,1.2,1.5]:
    vals=[rate(np.array([r*np.cos(p),r*np.sin(p)]),mu,R) for p in np.linspace(0,2*np.pi,65)[:-1]]
    signs.append({'r':r,'mean_dE':float(np.mean(vals)),'predicted_sign':int(np.sign(1-r*r))})
z=np.array([.6,.4]); mus=np.array([.1,.25,.5,.75,1.])
ys=np.array([rate(z,m,R) for m in mus]); slope=fit(mus,ys)
expected=(1-np.dot(z,z))*z[1]**2
reg=[]
for r0 in [.1,.4,.8,1.2,2.,4.]:
    tr=simulate(np.array([r0,0.]),30000,.01,omega,1.,R); tail=radius(tr)[-5000:]
    reg.append({'initial_r':r0,'tail_mean':float(tail.mean()),'tail_std':float(tail.std()),'max_r':float(radius(tr).max())})
# RK4 accuracy/stability sweep: compare h*omega to predicted linear stability limit ~2.785.
def amp(h):
    z=np.array([1.,0.]); tr=simulate(z,1000,h,omega,0.,R); return float(radius(tr)[-1])
steps=[]
for h in [.5,.9,1.2,1.4,1.5]: steps.append({'h_omega':h*omega,'final_radius':amp(h)})
result={'energy_sign_sweep':signs,'mu_slope_observed':slope,'mu_slope_predicted':float(expected),'regulation':reg,'rk4_stability_sweep':steps,'rk4_linear_boundary_predicted_homega':2.785}
open('results.json','w').write(json.dumps(result,indent=2))
print(json.dumps(result,indent=2))