"""Constant-wall liquid heating. Python 3 + Matplotlib; writes beside this file.
All values are illustrative calculations, not measurements.
"""
import csv
import json
import math
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

OUT = Path(__file__).resolve().parent
CP, UA, INLET, WALL = 4200., 10., 20., 80.

def temperature(flow, position=1):
    return WALL+(INLET-WALL)*math.exp(-UA*position/(flow*CP))

def numerical(flow, n):
    # Independent integration of dT/ds = UA/(mdot cp)*(Tw-T), s in [0,1].
    T, ds = INLET, 1/n
    f = lambda v: UA/(flow*CP)*(WALL-v)
    for _ in range(n):
        a=f(T); b=f(T+ds*a/2); c=f(T+ds*b/2); d=f(T+ds*c)
        T += ds*(a+2*b+2*c+d)/6
    return T

def main():
    rows=[]
    for i in range(1,51):
        flow=i/1000
        for j in range(101):
            s=j/100; T=temperature(flow,s)
            rows.append(dict(mass_flow_kg_s=flow,fraction_of_passage=s,temperature_C=T,heat_received_so_far_W=flow*CP*(T-INLET)))
    with (OUT/'liquid.csv').open('w',newline='',encoding='utf-8') as f:
        w=csv.DictWriter(f,fieldnames=rows[0]);w.writeheader();w.writerows(rows)
    error=max(abs(numerical(i/1000,100)-temperature(i/1000)) for i in range(1,51))
    assert error<1e-6
    fig,axes=plt.subplots(1,2,figsize=(10,4.5),layout='constrained')
    for flow,col in [(0.001,'#aa267f'),(.01,'#12658a'),(.05,'#826b00')]:
        r=[r for r in rows if r['mass_flow_kg_s']==flow]
        axes[0].plot([r['fraction_of_passage'] for r in r],[r['temperature_C'] for r in r],color=col,label=f'{flow:g} kg/s')
    axes[0].set(xlabel='Fraction of heated passage (0–1)',ylabel='Liquid temperature (°C)',title='a  Warming along the passage',ylim=(20,85))
    axes[0].axhline(WALL,ls='--',color='#555',label='Wall: 80 °C');axes[0].legend(fontsize=8)
    r=[r for r in rows if r['fraction_of_passage']==1]
    axes[1].plot([r['mass_flow_kg_s'] for r in r],[r['heat_received_so_far_W'] for r in r],color='#12658a')
    axes[1].set(xlabel='Mass flow rate (kg/s)',ylabel='Total heat-transfer rate (W)',title='b  More mass can carry more energy')
    for ax in axes: ax.grid(alpha=.2);ax.spines[['top','right']].set_visible(False)
    fig.suptitle('Constant-wall teaching model · UA = 10 W/K · cp = 4,200 J/(kg K)',fontsize=11)
    plt.rcParams['svg.fonttype']='none'
    for ext in ('svg','png','pdf'):fig.savefig(OUT/f'liquid-study.{ext}',dpi=180)
    plt.close(fig)
    report=dict(rows=len(rows),source='https://web.mit.edu/16.unified/www/SPRING/propulsion/notes/node131.html',constants=dict(cp_J_kgK=CP,UA_W_K=UA,inlet_C=INLET,wall_C=WALL),model='T(s)=Tw+(Tin-Tw)exp[-UA*s/(mdot*cp)]',verification=dict(RK4_100_segments_max_outlet_error_K=error),limits=['Constant wall temperature and conductance, independent of flow rate in this comparison.','No phase change, axial conduction, pressure drop, entrance effects or heater control.','A colder outlet at larger mass flow does not necessarily mean a lower total heat-transfer rate.'])
    (OUT/'liquid-study.json').write_text(json.dumps(report,indent=2)+'\n',encoding='utf-8')
    print(json.dumps(report['verification']))

if __name__=='__main__':main()