1import json, numpy as np
 2from vdp_cell import simulate, radius
 3
 4def rate(z,mu,R):
 5    n=len(z)//2; x=z[:n]; y=z[n:]
 6    return mu*(1-np.dot(z,z)/(R*R))*np.dot(y,y)
 7
 8def fit(x,y):
 9    return float(np.dot(x-x.mean(),y-y.mean())/np.dot(x-x.mean(),x-x.mean()))
10R=1.; omega=2.; mu=1.
11signs=[]
12for r in [.5,.8,1.,1.2,1.5]:
13    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]]
14    signs.append({'r':r,'mean_dE':float(np.mean(vals)),'predicted_sign':int(np.sign(1-r*r))})
15z=np.array([.6,.4]); mus=np.array([.1,.25,.5,.75,1.])
16ys=np.array([rate(z,m,R) for m in mus]); slope=fit(mus,ys)
17expected=(1-np.dot(z,z))*z[1]**2
18reg=[]
19for r0 in [.1,.4,.8,1.2,2.,4.]:
20    tr=simulate(np.array([r0,0.]),30000,.01,omega,1.,R); tail=radius(tr)[-5000:]
21    reg.append({'initial_r':r0,'tail_mean':float(tail.mean()),'tail_std':float(tail.std()),'max_r':float(radius(tr).max())})
22# RK4 accuracy/stability sweep: compare h*omega to predicted linear stability limit ~2.785.
23def amp(h):
24    z=np.array([1.,0.]); tr=simulate(z,1000,h,omega,0.,R); return float(radius(tr)[-1])
25steps=[]
26for h in [.5,.9,1.2,1.4,1.5]: steps.append({'h_omega':h*omega,'final_radius':amp(h)})
27result={'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}
28open('results.json','w').write(json.dumps(result,indent=2))
29print(json.dumps(result,indent=2))