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
# All numbers are deliberately hypothetical teaching inputs, not a named material.
M=.02; CP=500.; LATENT=50000.; T0=20.; TM=50.; POWER=10.
AREA=.1; AVAILABLE=2.; WVTRS=[.2,1.,5.]
def energy(m,cp,latent,t0,tm,p):
    if min(m,cp,latent,p)<=0 or tm<t0:raise ValueError('Positive properties/power and Tm>=T0 required')
    sensible=m*cp*(tm-t0);fusion=m*latent
    return sensible,fusion,(sensible+fusion)/p

def moisture(wvtr,area,days,available):
    if min(wvtr,area)<=0 or min(days,available)<0:raise ValueError('Positive rate/area and nonnegative time/water required')
    rate=wvtr*area;exhaustion=available/rate;potential=rate*days;loss=min(available,potential)
    return dict(wvtr_g_m2_day=wvtr,area_m2=area,time_days=days,initial_available_water_g=available,potential_constant_rate_loss_g=potential,capped_water_loss_g=loss,remaining_water_g=available-loss,exhaustion_time_days=exhaustion,reservoir_exhausted=days>=exhaustion)
qs,ql,tt=energy(M,CP,LATENT,T0,TM,POWER)
heating=[]
for t in np.linspace(0,tt,131):
    q=POWER*t;temp=T0+q/(M*CP) if q<qs else TM;fraction=0 if q<qs else min(1,(q-qs)/ql)
    heating.append(dict(time_s=float(t),input_energy_J=float(q),temperature_C=float(temp),liquid_fraction=float(fraction),stage='warming' if q<qs else 'melting' if q<qs+ql else 'complete'))
sweep=[]
for mass in [.01,.02,.05]:
    for p in [5.,10.,20.]:
        a,b,t=energy(mass,CP,LATENT,T0,TM,p);sweep.append(dict(mass_kg=mass,ideal_coupled_power_W=p,sensible_energy_J=a,fusion_energy_J=b,total_energy_J=a+b,ideal_completion_time_s=t))
water=[moisture(w,AREA,float(t),AVAILABLE) for w in WVTRS for t in np.linspace(0,30,121)]
def write(name,rows):
    with (OUT/name).open('w',newline='',encoding='utf-8') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
write('inclusion-heating.csv',heating);write('heating-budget.csv',sweep);write('intact-film-water.csv',water)
# Independent enthalpy reconstruction from the exported state, not temperature alone.
max_e=max(abs(M*CP*(r['temperature_C']-T0)+M*LATENT*r['liquid_fraction']-r['input_energy_J']) for r in heating)
assert max_e<1e-9
assert qs==300 and ql==1000 and tt==130
assert energy(2*M,CP,LATENT,T0,TM,POWER)[2]==2*tt
assert energy(M,CP,LATENT,T0,TM,2*POWER)[2]==tt/2
for r in water:
    assert 0<=r['capped_water_loss_g']<=AVAILABLE
    assert abs(r['remaining_water_g']+r['capped_water_loss_g']-AVAILABLE)<1e-12
    # Unit conversion independently from g/(m² day) to kg/(m² s), back to grams.
    si_loss=r['wvtr_g_m2_day']/1000/86400*r['area_m2']*(r['time_days']*86400)*1000
    assert abs(si_loss-r['potential_constant_rate_loss_g'])<1e-12
for w in WVTRS:
    t=AVAILABLE/(AREA*w);assert abs(moisture(w,AREA,t,AVAILABLE)['remaining_water_g'])<1e-12
meta={'title':'Changing phase takes energy; slowing moisture takes a barrier','kind':'Hypothetical deterministic teaching inputs, not material measurements','heating_parameters':{'mass_kg':M,'solid_specific_heat_J_kg_K':CP,'specific_latent_heat_J_kg':LATENT,'initial_C':T0,'melting_C':TM,'ideal_coupled_power_W':POWER},'moisture_parameters':{'effective_intact_area_m2':AREA,'initial_available_water_g':AVAILABLE,'illustrative_wvtr_g_m2_day':WVTRS,'max_days':30,'boundary_conditions':'Fixed unspecified temperature and humidity conditions are assumed by the hypothetical WVTR inputs; no empirical WVTR or actual condition rating is claimed.'},'formulas':{'sensible':'Qs=m*cp*(Tm-T0)','fusion':'Ql=m*Lf','completion_time':'t=(Qs+Ql)/P','water':'loss=min(m_available,A*WVTR*t_days)','exhaustion':'t_exhaust=m_available/(A*WVTR)'},'row_counts':{'heating_trace':len(heating),'heating_budgets':len(sweep),'water':len(water)},'validation':{'enthalpy_reconstruction_max_error_J':max_e,'baseline_energy_J':[qs,ql],'baseline_time_s':tt,'mass_and_power_scaling':True,'water_conservation_bounds':True,'SI_unit_conversion_parity':True,'exhaustion_transition':True},'limitations':['All parameters are freely chosen examples, not properties of a named alloy, composite or film.','Uniform inclusion temperature, constant solid heat capacity, sharp equilibrium melting temperature, constant latent heat, perfect power coupling, no losses or matrix heating.','The trace stops at complete melting; there is no post-melt heat capacity, cooling or resolidification model.','Ideal completion time is a lower-time reference under the stated supplied power with no competing heat sources; real heat losses or surrounding material heat capacity increase requirements.','Liquid fraction is an energy-accounting state, not a predicted stiffness, leakage or durability curve.','WVTR is an imposed steady mass-transmission rate per area at fixed conditions, not permeability; thickness cannot be varied without a separate model or measured rate.','Water cap is a reservoir accounting stop: after exhaustion the imposed constant driving conditions cease to be physically supportable. Flat capped curves do not predict a dry package equilibrium.','No seals, holes, adsorption, changing humidity, initial diffusion transient, oxygen transmission, microbial reduction or shelf life is computed.'],'sources':[{'title':'OpenStax Physics, Phase Change and Latent Heat','url':'https://openstax.org/books/physics/pages/11-3-phase-change-and-latent-heat'},{'title':'NASA Glenn, Heat Transfer','url':'https://www.grc.nasa.gov/WWW/BGH/heat.html'},{'title':'NIST/NBS, Permeability of membranes to water vapor with special reference to packaging materials','url':'https://www.nist.gov/publications/permeability-membranes-water-vapor-special-reference-packaging-materials'}]}
(OUT/'study.json').write_text(json.dumps(meta,indent=2)+'\n',encoding='utf-8')
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False})
fig,ax=plt.subplots(1,3,figsize=(13.6,4.7),constrained_layout=True)
ax[0].plot([r['time_s'] for r in heating],[r['temperature_C'] for r in heating],lw=2,color='#393f94');ax[0].axvline(qs/POWER,ls=':',color='#666');ax[0].text(55,40,'Energy goes into melting\nwhile temperature stays fixed',fontsize=9);ax[0].set(title='A  |  Warming, then melting',xlabel='Time under ideal 10 W input (s)',ylabel='Inclusion temperature (°C)')
ax[1].bar(['Warm by 30 K','Melt at 50°C'],[qs,ql],color=['#393f94','#b42e8d']);ax[1].set(title='B  |  Latent heat is part of the budget',ylabel='Required ideal energy (J)');ax[1].text(.5,1100,'Total 1,300 J → 130 s at 10 W',ha='center',fontsize=9);ax[1].set_ylim(0,1200)
for w,c in zip(WVTRS,['#393f94','#087f8c','#b42e8d']):
    rr=[r for r in water if r['wvtr_g_m2_day']==w];before=[r for r in rr if not r['reservoir_exhausted']];ex=AVAILABLE/(AREA*w)
    xx=[r['time_days'] for r in before]+([ex] if ex<=30 else []);yy=[r['capped_water_loss_g'] for r in before]+([AVAILABLE] if ex<=30 else [])
    ax[2].plot(xx,yy,lw=2,color=c,label=f'{w:g} g/(m² day)')
    if ex<=30:ax[2].plot([ex,30],[AVAILABLE,AVAILABLE],ls=':',color=c);ax[2].plot(ex,AVAILABLE,'o',color=c)
ax[2].set(title='C  |  Never lose more water than exists',xlabel='Time (days)',ylabel='Cumulative water loss (g)',ylim=(0,2.3));ax[2].legend(fontsize=8);ax[2].text(8,.7,'Dotted: accounting cap after exhaustion;\nconstant-rate premise has ended.',fontsize=8,bbox=dict(facecolor='white',edgecolor='none',alpha=.95))
for a in ax:a.grid(alpha=.2,axis='y')
fig.suptitle('Illustrative material budgets — no measured material properties',fontsize=15,weight='bold')
fig.supxlabel('Heating: 20 g hypothetical inclusion. Moisture: 0.1 m² intact film, 2 g available water; hypothetical rates at fixed conditions.',fontsize=9)
for ext in ['png','svg','pdf']:fig.savefig(OUT/f'thermal-barrier.{ext}',dpi=180)
print(json.dumps(meta['validation']))