from pathlib import Path
import math,csv,json,hashlib,sys
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
P=Path(__file__).resolve().parent

def tracer(a,t):
 if not all(math.isfinite(x) for x in (a,t)) or a<0 or t<0:raise ValueError('Nonnegative finite rate and time required')
 return math.exp(-a*t)
def thermal(t,tau):
 if not all(math.isfinite(x) for x in (t,tau)) or t<0 or tau<=0:raise ValueError('Finite nonnegative time and positive time constant required')
 return -math.expm1(-t/tau)
tr=[dict(exchange_per_hour=j/10,time_hours=i/20,normalized_concentration=tracer(j/10,i/20))for j in range(31)for i in range(121)]
th=[dict(time_constant_units=tau,time_units=i/10,normalized_response=thermal(i/10,tau))for tau in (1,3,10)for i in range(301)]
for name,rows in [('tracer.csv',tr),('thermal.csv',th)]:
 with (P/name).open('w',newline='',encoding='utf8')as f:
  w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
# Independently integrate the defining ODE by explicit midpoint; demonstrate second order convergence.
checks=[]
for kind,rate,end in [('tracer',.5,6),('tracer',3,6),('thermal',1,3),('thermal',.1,30)]:
 errors=[]
 for steps in (120,240,480):
  dt=end/steps;y=1. if kind=='tracer'else 0.
  def rhs(v):return -rate*v if kind=='tracer'else rate*(1-v)
  for _ in range(steps):y+=dt*rhs(y+dt*rhs(y)/2)
  reference=tracer(rate,end)if kind=='tracer'else thermal(end,1/rate)
  errors.append(abs(y-reference))
 assert errors[0]>errors[1]>errors[2] and 3.5<errors[1]/errors[2]<4.5
 checks.append(dict(model=kind,rate=rate,end=end,steps=[120,240,480],absolute_errors=errors))
assert tracer(0,6)==1 and tracer(3,0)==1 and thermal(0,1)==0
assert math.isclose(thermal(1,1),1-math.exp(-1))
assert tracer(3,1000)==0 and thermal(1000,1)==1
for j in range(31):
 vals=[r['normalized_concentration']for r in tr[j*121:(j+1)*121]]
 assert all(0<=v<=1 for v in vals)and all(x>=y for x,y in zip(vals,vals[1:]))
for j in range(3):
 vals=[r['normalized_response']for r in th[j*301:(j+1)*301]]
 assert all(0<=v<=1 for v in vals)and all(x<=y for x,y in zip(vals,vals[1:]))
guards=0
for f,args in [(tracer,(-1,1)),(tracer,(1,-1)),(tracer,(math.nan,1)),(tracer,(1,math.inf)),(thermal,(-1,1)),(thermal,(1,0)),(thermal,(1,-1)),(thermal,(1,math.inf))]:
 try:f(*args)
 except ValueError:guards+=1
 else:raise AssertionError(args)
plt.rcParams.update({'font.size':11,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':180})
fig,ax=plt.subplots(figsize=(8,5),layout='constrained')
for rate,style in zip((0,.5,1,2),(':','--','-','-.')):
 t=np.linspace(0,6,121);ax.plot(t,[tracer(rate,x)for x in t],ls=style,label=f'Assigned exchange = {rate:g} h⁻¹')
ax.set(xlabel='Time (hours)',ylabel='Normalized tracer concentration C/C₀',ylim=(-.02,1.05),title='An inert tracer leaves a perfectly mixed room\nNo ongoing source; incoming air has zero tracer');ax.legend();ax.grid(alpha=.2)
for ext in ('png','svg','pdf'):fig.savefig(P/f'tracer-balance.{ext}')
plt.close(fig)
fig,ax=plt.subplots(figsize=(8,5),layout='constrained')
for tau,style in zip((1,3,10),('-','--','-.')):
 t=np.linspace(0,30,301);ax.plot(t,[thermal(x,tau)for x in t],ls=style,label=f'Assigned time constant = {tau:g} units')
ax.set(xlabel='Time (arbitrary time units)',ylabel='Normalized thermal response θ',ylim=(-.02,1.05),title='Different assigned response times\nIndependent lumped examples; no spatial temperature prediction');ax.legend();ax.grid(alpha=.2)
for ext in ('png','svg','pdf'):fig.savefig(P/f'thermal-response.{ext}')
plt.close(fig)
(P/'validation.json').write_text(json.dumps(dict(all_pass=True,tracer_rows=len(tr),thermal_rows=len(th),ode_convergence=checks,guards=guards,limits_monotonicity_bounds=True),indent=2)+'\n')
(P/'study.json').write_text(json.dumps(dict(runtime=dict(python=sys.version,numpy=np.__version__,matplotlib=matplotlib.__version__),generator_sha256=hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),synthetic=True,tracer=dict(rate_per_hour=[0,3,.1],time_hours=[0,6,.05],formula='C/C0=exp(-a*t)',assumptions=['constant volume and exchange','perfect mixing','no source','zero incoming tracer','nonreactive tracer']),thermal=dict(time_constants=[1,3,10],time_units=[0,30,.1],formula='theta=1-exp(-t/tau)',assumptions=['independent first-order lumps','constant boundary condition','arbitrary consistent time units','no spatial field']),not_a_health_or_treatment_model=True),indent=2)+'\n')
print('PASS: 3751 tracer + 903 thermal = 4654 rows')