from pathlib import Path
import csv,json,math
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
OUT=Path(__file__).resolve().parent
ALPHA=1e-7; L=.01; TI=20.; TS=5.

def theta(x,fo,terms=64):
    x=np.asarray(x,dtype=float)
    n=np.arange(terms); lam=(2*n+1)*np.pi/2
    return np.sum((4/np.pi*(-1.)**n/(2*n+1)*np.exp(-lam*lam*fo))[:,None]*np.cos(lam[:,None]*x.reshape(1,-1)),axis=0)

def mean_theta(fo,terms=64):
    n=np.arange(terms); k=2*n+1
    return float(np.sum(8/(np.pi*np.pi*k*k)*np.exp(-k*k*np.pi*np.pi*fo/4)))

def fd_check(intervals):
    x=np.linspace(-1,1,intervals+1); dx=2/intervals
    u=np.ones(intervals+1);u[[0,-1]]=0
    fo=0.;results=[]
    for target in [.01,.05,.2,1.]:
        steps=math.ceil((target-fo)/(.4*dx*dx));dt=(target-fo)/steps;r=dt/dx**2
        for _ in range(steps):u[1:-1]+=r*(u[:-2]-2*u[1:-1]+u[2:])
        exact=theta(x,target)
        results.append({'fourier_number':target,'max_temperature_error_C':float(np.max(np.abs(u-exact))*(TI-TS)),'centre_temperature_error_C':float(abs(u[intervals//2]-exact[intervals//2])*(TI-TS)),'time_step_Fo':dt,'stability_ratio':r})
        fo=target
    return {'intervals':intervals,'results':results}

def write_csv(name,fields,rows):
    with (OUT/name).open('w',newline='',encoding='utf-8') as f:
        w=csv.DictWriter(f,fieldnames=fields);w.writeheader();w.writerows(rows)

x=np.linspace(-1,1,201);fos=[.01,.05,.2,1.]
profiles=[]
for fo in fos:
    vals=TS+(TI-TS)*theta(x,fo)
    for xx,v in zip(x,vals):profiles.append({'fourier_number':fo,'time_s':fo*L*L/ALPHA,'position_mm':xx*L*1000,'temperature_C':float(v)})
write_csv('slab-profiles.csv',list(profiles[0]),profiles)
time=[]
for fo in np.geomspace(.001,2.,241):
    time.append({'fourier_number':float(fo),'time_s':float(fo*L*L/ALPHA),'centre_C':float(TS+(TI-TS)*theta([0],fo)[0]),'surface_C':TS,'mean_C':TS+(TI-TS)*mean_theta(fo)})
write_csv('slab-time.csv',list(time[0]),time)
convergence=max(float(np.max(np.abs(theta(x,fo,64)-theta(x,fo,128)))*(TI-TS)) for fo in np.geomspace(.001,2,100))
checks=[fd_check(n) for n in [100,200,400]]
assert convergence<1e-10
assert all(r['max_temperature_error_C']<.015 for r in checks[-1]['results'])
for i in range(4):assert checks[-1]['results'][i]['max_temperature_error_C']<checks[0]['results'][i]['max_temperature_error_C']/10
assert min(row['temperature_C'] for row in profiles)>=TS-1e-10
assert max(row['temperature_C'] for row in profiles)<=TI+1e-10
assert all(time[i+1]['centre_C']<=time[i]['centre_C']+1e-10 for i in range(len(time)-1))
meta={'model':'One-dimensional constant-property inert slab; two fixed-temperature faces; analytical Fourier series','parameters':{'thermal_diffusivity_m2_s':ALPHA,'half_thickness_m':L,'initial_temperature_C':TI,'fixed_face_temperature_C':TS,'characteristic_time_s':L*L/ALPHA},'series_terms':64,'minimum_Fo':.001,'formula':'theta=sum_n 4*(-1)^n/((2*n+1)*pi)*cos((2*n+1)*pi*x/(2*L))*exp(-(2*n+1)^2*pi^2*Fo/4)','validation':{'64_vs_128_terms_max_difference_C':convergence,'explicit_finite_difference':checks,'bounds_and_monotonic_centre_checks':'passed'},'limitations':['Synthetic calculation, not experimental measurements.','No phase change, fluid motion, heat production, finite convection resistance or temperature-dependent properties.','The surface jumps instantaneously at t=0; plotted domain is t>0.','An infinitely wide plane wall removes lateral edge losses.','No biological, food safety or cryopreservation claims.']}
meta['source']='https://ocw.mit.edu/courses/3-044-materials-processing-spring-2013/6b75f358c8da03e1a3ea70a6e2098dd5_MIT3_044S13_TranHeaCondSol.pdf'
(OUT/'slab-study.json').write_text(json.dumps(meta,indent=2),encoding='utf-8')
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False,'axes.labelcolor':'#15253c','text.color':'#15253c','axes.titleweight':'bold','savefig.facecolor':'white'})
fig,ax=plt.subplots(1,2,figsize=(11,4.5),constrained_layout=True)
colors=['#353c92','#087f8c','#b02b8d','#db7410']
for fo,c in zip(fos,colors):ax[0].plot(x*L*1000,TS+(TI-TS)*theta(x,fo),label=f'{fo*L*L/ALPHA:g} s (Fo={fo:g})',color=c,lw=2)
ax[0].set(xlabel='Position across slab (mm)',ylabel='Temperature (°C)',title='A  |  The faces cool before the centre',ylim=(4.5,20.5));ax[0].legend(fontsize=8);ax[0].grid(alpha=.18)
for key,label,c,style in [('centre_C','Centre','#353c92','-'),('mean_C','Slab average','#087f8c','--'),('surface_C','Both faces held at 5°C','#b02b8d',':')]:ax[1].semilogx([r['time_s'] for r in time],[r[key] for r in time],label=label,color=c,ls=style,lw=2)
ax[1].set(xlabel='Time after face-temperature change (s)',ylabel='Temperature (°C)',title='B  |  A surface reading can miss stored heat',ylim=(4.5,20.5));ax[1].legend(fontsize=8);ax[1].grid(alpha=.18)
fig.suptitle('Cooling an inert slab: spatial heat diffusion',fontsize=15,weight='bold')
for ext in ['png','svg','pdf']:fig.savefig(OUT/f'slab-cooling.{ext}',dpi=180)
print(json.dumps(meta['validation'],indent=2))