"""Restricted classical stability and homogeneous energy accounting; dimensionless."""
from pathlib import Path
import json, hashlib
import numpy as np
import sympy as s
from scipy.integrate import solve_ivp
from scipy.optimize import minimize_scalar
out=Path(__file__).resolve().parent
x,r,q=s.symbols('x r q', real=True)
U=x*x*(1-x)**2
Vc=U+r*(1-x)**2
Vq=U+r*s.sqrt(q*q+(1-x)**4)
checks={
 'Uprime':s.simplify(s.diff(U,x)-2*x*(1-x)*(1-2*x))==0,
 'Usecond':s.simplify(s.diff(U,x,2)-2*(1-6*x+6*x*x))==0,
 'cold_derivative':s.simplify(s.diff(Vc,x)+2*(1-x)*(x*(2*x-1)+r))==0,
 'finite_p_derivative':s.simplify(s.diff(Vq,x)-2*(1-x)*(x*(1-2*x)-r*(1-x)**2/s.sqrt(q*q+(1-x)**4)))==0,
}
a=(1-s.sqrt(1-8*r))/4; b=(1+s.sqrt(1-8*r))/4
checks['cold_barrier']=s.simplify(Vc.subs(x,b)-Vc.subs(x,a)-(1-8*r)**s.Rational(3,2)/16)==0
assert all(checks.values())
rows=[]
for rho in [0,.05,.1,.125,.2]:
 roots=[1.]
 if rho<=.125:
  d=np.sqrt(1-8*rho); roots=[(1-d)/4,(1+d)/4,1.]
 for a in sorted(set(roots)):
  rows.append({'r':rho,'chi':a,'mass_ratio':(1-a)**2,'curvature':2*(1-6*a+6*a*a)+2*rho,'potential':a*a*(1-a)**2+rho*(1-a)**2})
fold=[]
for p in [0,.3,1,3]:
 def required_density(a):return a*(1-2*a)*np.sqrt(p*p+(1-a)**4)/(1-a)**2
 z=minimize_scalar(lambda a:-required_density(a),bounds=(0,.5),method='bounded',options={'xatol':1e-13})
 fold.append({'p_over_m0':p,'chi_fold':z.x,'r_fold':-z.fun})
# Independent integrations of homogeneous, unforced vacuum scalar trajectories.
trajectories=[]
for x0 in [.01,.9, .99]:
 def rhs(t,y):return [y[1],-2*y[0]*(1-y[0])*(1-2*y[0])]
 ts=np.linspace(0,30,3001)
 sol=solve_ivp(rhs,(0,30),[x0,0],method='DOP853',t_eval=ts,rtol=1e-11,atol=1e-13,max_step=.04)
 assert sol.success
 chi,v=sol.y; energy=.5*v*v+chi*chi*(1-chi)**2
 err=float(np.max(abs(energy-energy[0])))
 assert err<1e-9
 trajectories.append({'chi_initial':x0,'chi_min':float(chi.min()),'chi_max':float(chi.max()),'energy_error':err})
 np.savetxt(out/'data'/f'vacuum_{x0}.csv',np.c_[ts,chi,v,energy],delimiter=',',header='t,chi,chi_dot,energy',comments='')
# Prescribed smooth conversion cycle with explicitly integrated generalized work.
# Field plus particle energy density: E=.5 chi_dot^2+U+n sqrt(p^2+m^2).
p=.3;n=.1

def path(t):
 if t<10: a=t/10; sign=1.; base=0.
 elif t<=15:return 1.,0.,0.
 else:a=(t-15)/10;sign=-1.;base=1.
 return base+sign*(6*a**5-15*a**4+10*a**3),sign*(30*a*a*(a-1)**2)/10,sign*(120*a**3-180*a*a+60*a)/100

def energy(t):
 a,v,_=path(t);return .5*v*v+a*a*(1-a)**2+n*np.sqrt(p*p+(1-a)**4)

def work_rhs(t,y):
 a,v,acc=path(t);H=np.sqrt(p*p+(1-a)**4)
 particle_source=-2*n*(1-a)**3/H
 Jext=acc+2*a*(1-a)*(1-2*a)+particle_source
 return [Jext*v,particle_source*v]
for step in [.1,.05]:
 sums=np.zeros(2);maxerr=0.;samples=[]
 for ta,tb in [(0,10),(10,15),(15,25)]:
  ts=np.linspace(ta,tb,1001)
  sol=solve_ivp(work_rhs,(ta,tb),sums,method='DOP853',t_eval=ts,max_step=step,rtol=1e-11,atol=1e-13)
  assert sol.success
  sums=sol.y[:,-1]
  e=np.array([energy(t) for t in ts]);maxerr=max(maxerr,float(np.max(abs(e-energy(0)-sol.y[0]))))
  samples.extend(zip(ts,e,sol.y[0],sol.y[1]))
 if step==.05:
  np.savetxt(out/'data'/'controlled_cycle.csv',samples,delimiter=',',header='t,field_plus_particle_energy,external_work,particle_energy_change',comments='')
  cycle={'p':p,'n':n,'energy_initial':energy(0),'energy_platform':energy(12),'external_work_entry':energy(12)-energy(0),'net_external_work':float(sums[0]),'net_particle_work':float(sums[1]),'energy_balance_max_error':maxerr,'field_kinetic_peak':.5*.1875**2,'unloaded_potential_barrier':1/16}
 assert maxerr<1e-9
cold_r=.05;chi=(1-np.sqrt(1-8*cold_r))/4;K=2*(1-6*chi+6*chi*chi)+2*cold_r;k=.1;C=4*cold_r;A=k*k+K
wm=(A-np.sqrt(A*A+4*C*k*k))/2
mass=(1-chi)**2; mp=-2*(1-chi)
linear_matrix=np.array([[0,1,0,0],[-k*k-K,0,-mp,0],[0,0,0,1],[-cold_r*mp/mass*k*k,0,0,0]])
growth=np.linalg.eigvals(linear_matrix).real.max()
assert abs(growth-np.sqrt(-wm))<1e-12
result={'assumptions':'c=hbar=f=Lambda=m0=1; fixed density results are restricted curvature tests, not full coupled stability','symbolic_checks':checks,'cold_equilibria':rows,'finite_momentum_folds':fold,'vacuum_dynamics':trajectories,'controlled_cycle':cycle,'dust_instability_example':{'r':cold_r,'chi':chi,'k':k,'K':K,'omega_minus_squared':wm,'growth_rate':float(np.sqrt(-wm))},'source_sha256':hashlib.sha256(Path(__file__).read_bytes()).hexdigest()}
(out/'data'/'results.json').write_text(json.dumps(result,ensure_ascii=False,indent=2))
print(json.dumps(result,ensure_ascii=False,indent=2))