#!/usr/bin/env python3
"""Andrade Safe | Engenharia 003 — perda térmica em cilindro com isolante.

Modelo didático, regime permanente. Resistências em série por convecção interna,
condução radial aço e isolante; fronteira externa: convecção + radiação (não-linear).
Python 3.9+ com biblioteca padrão. Para gráficos, requer matplotlib (opcional).

Uso: python Laboratorio_Engenharia_003_Modelo_Termico.py
     python Laboratorio_Engenharia_003_Modelo_Termico.py --graficos
"""
from __future__ import annotations
import argparse
import csv
import math
from dataclasses import dataclass, replace
from pathlib import Path

SIGMA = 5.670374419e-8  # W/m2/K4, Stefan-Boltzmann

@dataclass(frozen=True)
class Parametros:
    comprimento_m: float = 10.0
    raio_interno_m: float = 0.0125
    raio_tubo_m: float = 0.015
    espessura_m: float = 0.02
    temp_fluido_C: float = 90.0
    temp_ambiente_C: float = 25.0
    h_interno: float = 500.0
    h_externo: float = 10.0
    k_aco: float = 45.0
    k_isolante: float = 0.04
    emissividade: float = 0.80
    radiacao: bool = True
    horas_ano: float = 4000
    tarifa_R_kWh: float = 0.75
    eficiencia_geracao: float = 0.85

@dataclass(frozen=True)
class Resultado:
    perda_W: float
    perda_W_m: float
    T_superficie_C: float
    T_interna_tubo_C: float
    T_externa_tubo_C: float
    R_in: float
    R_aco: float
    R_isolante: float
    q_conveccao_W: float
    q_radiacao_W: float
    perda_energia_termica_kWh_ano: float
    custo_energia_reposicao_R_ano: float


def calcular(p: Parametros) -> Resultado:
    if not (p.comprimento_m > 0 and p.raio_interno_m > 0 and p.raio_tubo_m > p.raio_interno_m
            and p.espessura_m >= 0 and p.h_interno > 0 and p.h_externo > 0
            and p.k_aco > 0 and p.k_isolante > 0 and 0 <= p.emissividade <= 1
            and 0 < p.eficiencia_geracao <= 1 and p.horas_ano >= 0 and p.tarifa_R_kWh >= 0
            and p.temp_fluido_C > p.temp_ambiente_C and p.temp_ambiente_C > -273.15):
        raise ValueError('Parametros fora do dominio do modelo. Requer fluido mais quente que ambiente.')
    r1, r2, r3, L = p.raio_interno_m, p.raio_tubo_m, p.raio_tubo_m+p.espessura_m, p.comprimento_m
    R_in = 1 / (p.h_interno * 2 * math.pi * r1 * L)
    R_aco = math.log(r2/r1) / (2*math.pi*p.k_aco*L)
    R_iso = math.log(r3/r2) / (2*math.pi*p.k_isolante*L)
    R_int = R_in + R_aco + R_iso
    area = 2*math.pi*r3*L
    ta = p.temp_ambiente_C
    delta = p.temp_fluido_C-ta

    def fronteira(Ts_C: float):
        qconv = p.h_externo*area*(Ts_C-ta)
        qrad = (p.emissividade*SIGMA*area*((Ts_C+273.15)**4-(ta+273.15)**4)) if p.radiacao else 0.0
        return qconv, qrad
    # (Tf - Ts)/R_int = Q_conv(Ts)+Q_rad(Ts).
    lo, hi = ta, p.temp_fluido_C
    for _ in range(90):
        ts = (lo+hi)/2
        qc, qr=fronteira(ts)
        balance = (p.temp_fluido_C-ts)/R_int-(qc+qr)
        if balance > 0:
            lo=ts
        else:
            hi=ts
    ts=(lo+hi)/2
    qc,qr=fronteira(ts)
    q=(p.temp_fluido_C-ts)/R_int
    tin=p.temp_fluido_C-q*R_in
    tout=tin-q*R_aco
    energy=q/1000*p.horas_ano
    return Resultado(q, q/L, ts, tin, tout, R_in, R_aco, R_iso, qc, qr,
                     energy, energy/p.eficiencia_geracao*p.tarifa_R_kWh)


def perfil_radial(p: Parametros, n: int=180):
    s=calcular(p)
    r1,r2,r3=p.raio_interno_m,p.raio_tubo_m,p.raio_tubo_m+p.espessura_m
    rlist=[r1+(r3-r1)*i/(n-1) for i in range(n)]
    tlist=[]
    for r in rlist:
        if r<=r2:
            t=s.T_interna_tubo_C-s.perda_W*math.log(r/r1)/(2*math.pi*p.k_aco*p.comprimento_m)
        else:
            t=s.T_externa_tubo_C-s.perda_W*math.log(r/r2)/(2*math.pi*p.k_isolante*p.comprimento_m)
        tlist.append(t)
    return rlist,tlist


def produzir_csv(dest: Path, p: Parametros):
    with dest.open('w',newline='',encoding='utf-8-sig') as f:
        w=csv.writer(f,delimiter=';')
        w.writerow(['espessura_mm','perda_W','perda_W_m','temperatura_superficie_C',
                    'conveccao_W','radiacao_W','custo_reposicao_R_ano', 'economia_relativa_percentual'])
        base=calcular(replace(p,espessura_m=0))
        for mm in range(0,51,2):
            s=calcular(replace(p,espessura_m=mm/1000))
            w.writerow([mm,round(s.perda_W,3),round(s.perda_W_m,3),round(s.T_superficie_C,3),
                        round(s.q_conveccao_W,3),round(s.q_radiacao_W,3),round(s.custo_energia_reposicao_R_ano,2),
                        round(100*(1-s.perda_W/base.perda_W),2)])


def gerar_graficos(dest: Path,p: Parametros):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    from matplotlib.ticker import AutoMinorLocator
    plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,
                         'axes.spines.right':False,'axes.grid':True,'grid.alpha':.16,
                         'savefig.facecolor':'white','figure.facecolor':'white'})
    azul,laranja,verde,cinza='#1467a0','#e36a31','#208772','#364152'
    x=list(range(0,51))
    loss=[calcular(replace(p,espessura_m=i/1000)).perda_W for i in x]
    base=loss[0]
    # gráfico 01
    fig, ax=plt.subplots(figsize=(10.6,5.7))
    ax.plot(x,loss,lw=3,color=azul)
    ax.scatter([0,10,20,40],[loss[i] for i in [0,10,20,40]],s=60,color=laranja,zorder=4)
    for i in [0,10,20,40]:
        ax.annotate(f'{loss[i]:.0f} W',(i,loss[i]), xytext=(6,9),textcoords='offset points', fontsize=9)
    ax.set(xlabel='Espessura do isolamento (mm)',ylabel='Perda térmica nos 10 m (W)',
           title='Tubulação a 90 °C: quanto isolamento reduz a perda de calor?')
    ax.set_xlim(0,50); ax.set_ylim(0,max(loss)*1.12)
    fig.text(.12,.015,'Modelo radial estacionário • h externo constante • radiação resolvida na superfície • Andrade Safe',color=cinza,fontsize=9)
    fig.tight_layout(rect=[0,.04,1,1]); fig.savefig(dest/'grafico_01_perda_espessura.png',dpi=170);plt.close(fig)
    # gráfico 02
    fig,ax=plt.subplots(figsize=(10.6,5.7))
    for e,color in [(0,laranja),(10,'#e4a44c'),(20,azul),(40,verde)]:
        pp=replace(p,espessura_m=e/1000)
        radii,ts=perfil_radial(pp)
        ax.plot([r*1000 for r in radii],ts,label=f'{e} mm de isolamento',lw=2.6,color=color)
    ax.axvline(p.raio_tubo_m*1000,ls='--',color=cinza,alpha=.5)
    ax.axhline(p.temp_ambiente_C,ls=':',color=cinza,alpha=.8,label='Ar ambiente = 25 °C')
    ax.set(xlabel='Distância radial ao eixo (mm)',ylabel='Temperatura (°C)',
           title='Queda de temperatura através do isolamento')
    ax.legend(loc='upper right',fontsize=9);ax.set_xlim(p.raio_interno_m*1000,60)
    fig.text(.12,.015,'Perfis radiais no aço e isolante; cada curva termina na sua superfície externa • Andrade Safe',color=cinza,fontsize=9)
    fig.tight_layout(rect=[0,.04,1,1]);fig.savefig(dest/'grafico_02_perfil_temperatura.png',dpi=170);plt.close(fig)
    # grafico 03
    fig,ax=plt.subplots(figsize=(10.6,5.7))
    for e,color in [(0,laranja),(10,'#e4a44c'),(20,azul),(40,verde)]:
        s=calcular(replace(p,espessura_m=e/1000))
        ax.bar(str(e),s.custo_energia_reposicao_R_ano,color=color,width=.58)
        ax.text(str(e),s.custo_energia_reposicao_R_ano+40,f'R$ {s.custo_energia_reposicao_R_ano:,.0f}',ha='center')
    ax.set(xlabel='Isolamento (mm)',ylabel='Custo anual estimado de reposição de energia (R$/ano)',
           title='Cenário ilustrativo: 4.000 h/ano, R$ 0,75/kWh e eficiência térmica de 85%')
    ax.set_ylim(0,max(calcular(replace(p,espessura_m=e/1000)).custo_energia_reposicao_R_ano for e in [0,10,20,40])*1.18)
    fig.text(.12,.015,'Não inclui custo de isolamento, manutenção, mudanças de processo nem preço real de combustível • Andrade Safe',color=cinza,fontsize=9)
    fig.tight_layout(rect=[0,.04,1,1]);fig.savefig(dest/'grafico_03_custo_operacional.png',dpi=170);plt.close(fig)
    # grafico 04 — radius critical special case
    r0=.001;ka=.04;he=10;Tc=90;Ta=25
    rs=[r0+i*.0002 for i in range(0,300)]
    qs=[(Tc-Ta)/(math.log(r/r0)/(2*math.pi*ka)+1/(he*2*math.pi*r)) for r in rs]
    fig,ax=plt.subplots(figsize=(10.6,5.7))
    ax.plot([r*1000 for r in rs],qs,color=azul,lw=3)
    ax.axvline(4,color=laranja,ls='--',label='Raio crítico = k/h = 4 mm')
    ax.set(xlabel='Raio externo resultante (mm)',ylabel='Perda térmica por metro (W/m)',
           title='Raio crítico: quando acrescentar isolante pode inicialmente elevar a perda?')
    ax.legend();fig.text(.12,.015,'EXEMPLO SEPARADO: fio de 1 mm de raio, sem radiação; h = 10 W/(m²·K), k = 0,04 W/(m·K)',color=cinza,fontsize=9)
    fig.tight_layout(rect=[0,.04,1,1]);fig.savefig(dest/'grafico_04_raio_critico.png',dpi=170);plt.close(fig)


def main():
    ap=argparse.ArgumentParser()
    ap.add_argument('--graficos',action='store_true',help='Gerar 4 graficos com matplotlib')
    ap.add_argument('--saida',type=Path,default=Path(__file__).parent)
    args=ap.parse_args();args.saida.mkdir(parents=True,exist_ok=True)
    p=Parametros(); base=calcular(replace(p,espessura_m=0))
    print('ENGENHARIA 003 | Modelo radial com conveccao + radiacao')
    for e in [0,10,20,40,50]:
        s=calcular(replace(p,espessura_m=e/1000))
        print(f'{e:2d} mm | Q={s.perda_W:8.2f} W | {s.perda_W_m:7.2f} W/m | ' 
              f'Tsup={s.T_superficie_C:5.2f} C | queda={(1-s.perda_W/base.perda_W)*100:5.2f}% | ' 
              f'custo=R$ {s.custo_energia_reposicao_R_ano:,.2f}/ano')
    print('Raio critico convectivo teorico:',p.k_isolante/p.h_externo*1000,'mm; raio tubo =',p.raio_tubo_m*1000,'mm')
    produzir_csv(args.saida/'tabela_003_cenarios.csv',p)
    if args.graficos: gerar_graficos(args.saida,p)

if __name__=='__main__':main()