"""4.3-cinder. Nonlinear porous-flow damping added to unforced two-mass dynamics.
Tanaka f=175/Re+1.6, hydraulic diameter eps*wire/(1-eps). Flow approximated
Q=Ad*vd; screen experimental calibration and full chamber mass-flow coupling
are absent. A bounded oscillation in this screen is not engine validation.
"""
import analyze as a
from scipy.integrate import solve_ivp,trapezoid
import numpy as np,math,json,platform
from pathlib import Path
D=a.D;O=a.OUT;area=math.pi*(D['regenerator_outer_mm']**2-D['regenerator_inner_mm']**2)*1e-6
por=D['regenerator_porosity'];wire=50e-6;dh=por*wire/(1-por);rho=D['mean_pressure_Pa']/(a.R*a.Tr);mu=3.2e-5;length=D['regenerator_length_mm']/1000
mass=D['mean_pressure_Pa']*sum(a.Vbase/a.T)/a.R
vbp=D['piston_bounce_cm3']*1e-6;vbd=D['displacer_bounce_cm3']*1e-6

def make(kd,rl):
 def f(t,y):
  xp,vp,xd,vd,i=y;V=a.volumes(xp,xd)
  if min(V)<=0 or vbp-a.Ap*xp<=0 or vbd-a.Ar*xd<=0:return np.zeros(5)
  p=mass*a.R/sum(V/a.T);pb=D['mean_pressure_Pa']*(vbp/(vbp-a.Ap*xp))**a.gamma;pd=D['mean_pressure_Pa']*(vbd/(vbd-a.Ar*xd))**a.gamma
  u=a.Ad*vd/(por*area);dp=.5*rho*length/dh*(175*mu/(rho*dh)*u+1.6*u*abs(u))
  fp=a.Ap*(p-pb)-D['piston_mechanical_spring_N_m']*xp-D['piston_friction_N_s_m']*vp-D['alternator_K_N_A']*i
  fd=a.Ar*(p-pd)-kd*xd-D['displacer_friction_N_s_m']*vd-a.Ad*dp
  return np.array([vp,fp/D['piston_moving_mass_kg'],vd,fd/D['displacer_moving_mass_kg'],(D['alternator_K_N_A']*vp-(rl+D['alternator_R_ohm'])*i)/D['alternator_L_H']])
 return f
if __name__=='__main__':
 rows=[];candidates=[]
 for kd in range(14000,40001,1000):
  for rl in [2,3,4,5,6,8,10,15,20,30]:
   f=make(kd,rl);eps=1e-7;J=np.column_stack([(f(0,np.eye(5)[k]*eps)-f(0,-np.eye(5)[k]*eps))/(2*eps) for k in range(5)]);eig=np.linalg.eigvals(J);growth=max(x.real for x in eig)
   row={'kd_N_m':kd,'load_ohm':rl,'max_growth_per_s':growth};rows.append(row)
   if 2<growth<9:candidates.append(row)
 # Preserve breadth of loading; no claim that a favorable candidate is guaranteed.
 selected=[]
 for rl in [2,3,4,5,6,8,10,15,20,30]:
  options=[r for r in candidates if r['load_ohm']==rl]
  if options:selected.append(min(options,key=lambda r:abs(r['max_growth_per_s']-4)))
 trials=[]
 for c in selected:
  f=make(c['kd_N_m'],c['load_ohm'])
  def hit(t,y):return min(.0065-abs(y[0]),.0045-abs(y[2]))
  hit.terminal=True;hit.direction=-1
  sol=solve_ivp(f,[0,2],[5e-5,0,5e-5,0,0],max_step=1/3500,rtol=1e-6,atol=1e-9,events=hit,dense_output=True)
  tt=np.linspace(0,sol.t[-1],12001);y=sol.sol(tt);tail=tt>max(0,tt[-1]-.25);before=(tt>max(0,tt[-1]-.5))&(tt<=max(0,tt[-1]-.25))
  ap=max(abs(y[0,tail]));ad=max(abs(y[2,tail]));old=max(abs(y[0,before])) if before.any() else 0
  u=a.Ad*y[3,tail]/(por*area);dp=.5*rho*length/dh*(175*mu/(rho*dh)*u+1.6*u*abs(u));pump=float(np.mean(dp*a.Ad*y[3,tail]));load=float(np.mean(y[4,tail]**2)*c['load_ohm']);cu=float(np.mean(y[4,tail]**2)*D['alternator_R_ohm'])
  result={**c,'end_s':float(tt[-1]),'overstroke':bool(len(sol.t_events[0])),'piston_peak_mm':float(ap*1000),'displacer_peak_mm':float(ad*1000),'amplitude_change_fraction_last_windows':float((ap-old)/max(ap,1e-10)),'load_W':load,'copper_W':cu,'regenerator_pressure_loss_W':pump,'rectified_minus_15W_aux_W':load*.96-15}
  result['stable_window_screen']=not result['overstroke'] and abs(result['amplitude_change_fraction_last_windows'])<.02
  trials.append(result);print(result,flush=True)
  np.savetxt(O/('dynamic-k%d-r%s.csv'%(c['kd_N_m'],c['load_ohm'])),np.c_[tt,y.T],delimiter=',',header='time_s,xp_m,vp_m_s,xd_m,vd_m_s,current_A',comments='')
 result={'version':'4.3-cinder','host':platform.node(),'linear_cases':len(rows),'linear_screen':rows,'nonlinear_trials':trials,'scope':'Engineering stability screening, uncalibrated isothermal model and approximate Tanaka pressure loss. No combustion, finite heat transfer, leakage, field solution or physical prototype acceptance.'}
 (O/'dynamic-screen.json').write_text(json.dumps(result,indent=2))