"""Independent continuous one-zone calculation. No App/Engine import or VICE execution."""
import math
import numpy as np
from scipy.integrate import quad,solve_ivp
from scipy.optimize import brentq

def histories(t,p):
 tau=p['tau'];mass=p['formedMass'];burst=p['burstFraction']
 return mass*(1-burst)*math.exp(-t/tau)/(tau*(-math.expm1(-8/tau)))+mass*burst*math.exp(-.5*((t-3.7)/.2)**2)/(.2*math.sqrt(2*math.pi)),p['inflowPeak']*math.exp(-.5*((t-3)/.3)**2)
def formed(t,p):return p['formedMass']*(1-p['burstFraction'])*(-math.expm1(-t/p['tau']))/(-math.expm1(-8/p['tau']))+p['formedMass']*p['burstFraction']*.5*(math.erf((t-3.7)/(.2*math.sqrt(2)))-math.erf(-3.7/(.2*math.sqrt(2))))
def incoming(t,p):return p['inflowPeak']*.3*math.sqrt(math.pi/2)*(math.erf((t-3)/(.3*math.sqrt(2)))-math.erf(-3/(.3*math.sqrt(2))))
def delay_cdf(t):
 if t<=.15:return 0.
 if t>=15:return 1.
 return (.15**-.1-t**-.1)/(.15**-.1-15**-.1)
def response(t,p,tight=False):
 if t<=.15:return 0.
 A=.1/(.15**-.1-15**-.1);points=[u for u in [t-4.1,t-3.7,t-3.3]if .15<u<t]
 return quad(lambda delay:histories(t-delay,p)[0]*A*delay**-1.1,.15,min(t,15),points=points,epsabs=2e-13 if tight else 2e-11,epsrel=2e-11 if tight else 2e-9)[0]
def simulate(p,fixed,save=.02,tight=False):
 keys=['formedMass','tau','eta','inflowPeak','burstFraction'];bounds=[(.2,.55),(.45,4.5),(0,1.5),(0,.7),(0,.35)]
 for key,(lower,upper)in zip(keys,bounds):
  if not math.isfinite(p[key])or not lower<=p[key]<=upper:raise ValueError('DOMAIN: '+key)
 R=fixed['returnFraction'];consumption=1-R+p['eta']
 def gas(t):return 1+incoming(t,p)-consumption*formed(t,p)
 def slope(t):a,b=histories(t,p);return b-consumption*a
 grid=np.linspace(0,8,1601);critical=[0,8]
 for a,b in zip(grid[:-1],grid[1:]):
  if slope(a)*slope(b)<0:critical.append(brentq(slope,a,b,xtol=1e-13))
 minimum=min(gas(t)for t in critical)
 if minimum<=0:raise ValueError('GAS_EXHAUSTED')
 if minimum<.05:raise ValueError('GAS_SUPPORT')
 def rhs(t,y):
  psi,flow=histories(t,p);delayed=response(t,p,tight);locked=(1-R)*psi;wind=p['eta']*psi;zmg,zfe=y[4]/y[0],y[5]/y[0]
  pm,pf=fixed['netCC']['mg']*psi,fixed['netCC']['fe']*psi;dm,df=fixed['netIa']['mg']*delayed,fixed['netIa']['fe']*delayed
  return [flow-locked-wind,locked,wind,flow,pm+dm-(locked+wind)*zmg,pf+df-(locked+wind)*zfe,locked*zmg,locked*zfe,wind*zmg,wind*zfe,pm,pf,dm,df]
 initial=[1,0,0,0,fixed['initialFeFraction']*fixed['netCC']['mg']/fixed['netCC']['fe'],fixed['initialFeFraction']]+[0]*8
 count=round(8/save);times=np.array([i*8/count for i in range(count+1)])
 # Integrate exactly to the known source-onset before evaluating later dense states.
 options=dict(method='DOP853',rtol=4e-12 if tight else 2e-10,atol=4e-15 if tight else 2e-13,max_step=.025 if tight else .05,dense_output=True)
 early=solve_ivp(rhs,(0,.15),initial,**options);late=solve_ivp(rhs,(.15,8),early.y[:,-1],**options)
 if not early.success or not late.success:raise ValueError('INTEGRATION')
 states=np.array([early.sol(t)if t<=.15 else late.sol(t)for t in times]);fractions=states[:,4:6]/states[:,0,None]
 brackets=np.column_stack((np.log10(fractions[:,1]/fixed['solarMassFractions']['fe']),np.log10((fractions[:,0]/fractions[:,1])/(fixed['solarMassFractions']['mg']/fixed['solarMassFractions']['fe']))))
 locked=np.diff(states[:,1]);mg=np.diff(states[:,6])/locked;fe=np.diff(states[:,7])/locked
 cohorts=np.column_stack((times[:-1],times[1:],locked,mg,fe,np.log10(fe/fixed['solarMassFractions']['fe']),np.log10((mg/fe)/(fixed['solarMassFractions']['mg']/fixed['solarMassFractions']['fe']))))
 if not np.all(np.isfinite(states))or not np.all(np.isfinite(cohorts)):raise ValueError('NONFINITE')
 return {'schema':'chemical-reader-record-v1','name':p['name'],'parameters':p,'saveEveryGyr':save,'timeGyr':times.tolist(),'stateFields':['gas','locked','wind','inflow','mgGas','feGas','mgLocked','feLocked','mgWind','feWind','mgCC','feCC','mgIa','feIa'],'states':states.tolist(),'abundances':brackets.tolist(),'flowFields':['formed','inflow','locked','wind','delayedResponse'],'flows':[[*histories(t,p),(1-R)*histories(t,p)[0],p['eta']*histories(t,p)[0],response(t,p,tight)]for t in times],'cohortFields':['startGyr','endGyr','lockedMass','magnesiumBirthFraction','ironBirthFraction','ironHydrogenDex','magnesiumIronDex'],'cohorts':cohorts.tolist(),'minimumGas':minimum,'units':{'mass':'1e9 solar masses','time':'Gyr','rate':'1e9 solar masses/Gyr = solar masses/year','abundance':'dex, fixedtracehydrogenAsplundreference'},'interpretation':'Independentownprescribedhistory, not observedstars or calibratedgalaxy history'}