# -*- coding: utf-8 -*-
"""
Converted from IPYNB to PY
"""

# %% [code] Cell 1
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

# ========================= INPUT PARAMETERS ==============================
# Grouping parameters in a dictionary to easily pass them into functions
params = {
    # Geometry
    'Di': 0.0471,             # Inner diameter (m)
    'Do': 0.0571,             # Outer diameter (m)
    
    # Flow
    'mdot': 0.053,            # kg/s
    
    # Exhaust gas properties
    'R': 287,                 # Gas constant (J/kg-K)
    'mu': 4.5e-5,             # Dynamic viscosity (kg/m-s)
    'k': 0.055,               # Thermal conductivity (W/m-K)
    'Pr': 0.70,               # Prandtl number
    
    # Environment
    'Tamb': 25 + 273.15,      # Ambient temperature (K)
    'VehicleSpeed': 20,       # Vehicle speed (m/s)
    
    # Radiation
    'eps': 0.80,
    'sigma': 5.670374419e-8
}

L = 4.812                     # Pipe length (m)
Tin_C = 800                   # Inlet Temperature (C)
Tin = Tin_C + 273.15          # Inlet Temperature (K)

# ========================= ODE FUNCTION ==================================
def exhaustODE(x, Tgas_arr, p):
    Tgas = Tgas_arr[0]
    
    # Gas Properties
    cp = 1000 + 0.1 * (Tgas - 300)
    rho = 130000 / (p['R'] * Tgas)
    Area = np.pi * p['Di']**2 / 4
    Velocity = p['mdot'] / (rho * Area)
    Re = rho * Velocity * p['Di'] / p['mu']
    Nu = 2 * 0.023 * Re**0.8 * p['Pr']**0.4
    hin = Nu * p['k'] / p['Di']
    
    # External Convection
    hout = 10 + 4 * p['VehicleSpeed']
    
    # Wall Temperature calculation
    def wallEq(Tw):
        return (hin * p['Di'] * (Tgas - Tw) 
                - hout * p['Do'] * (Tw - p['Tamb']) 
                - p['eps'] * p['sigma'] * p['Do'] * (Tw**4 - p['Tamb']**4))
    
    # fsolve returns an array, extract the first element
    Twall = fsolve(wallEq, x0=Tgas - 30)[0]
    
    # Gas Temperature Gradient
    dTdx = -(hin * np.pi * p['Di'] * (Tgas - Twall)) / (p['mdot'] * cp)
    return [dTdx]

# ========================= ODE SOLUTION ==================================
xspan = (0, L)

# Using dense_output or specifying a max_step ensures we get enough points for smooth plots
sol = solve_ivp(fun=lambda x, y: exhaustODE(x, y, params),
                t_span=xspan,
                y0=[Tin],
                method='RK45',
                rtol=1e-6,
                atol=1e-7,
                max_step=0.05) # Forces finer resolution for plotting

x = sol.t
Tgas = sol.y[0]

# =================== POST PROCESSING =====================================
N = len(x)
Twall = np.zeros(N)
hin   = np.zeros(N)
hout  = np.zeros(N)
Re    = np.zeros(N)
Nu    = np.zeros(N)
rho   = np.zeros(N)
Vel   = np.zeros(N)
cp    = np.zeros(N)
qConv = np.zeros(N)
qRad  = np.zeros(N)
qTotal = np.zeros(N)
Area = np.pi * params['Di']**2 / 4

for i in range(N):
    # Gas properties
    cp[i] = 1000 + 0.1 * (Tgas[i] - 300)
    rho[i] = 130000 / (params['R'] * Tgas[i])
    Vel[i] = params['mdot'] / (rho[i] * Area)
    Re[i] = rho[i] * Vel[i] * params['Di'] / params['mu']
    Nu[i] = 2 * 0.023 * Re[i]**0.8 * params['Pr']**0.4
    hin[i] = Nu[i] * params['k'] / params['Di']
    
    # External convection
    hout[i] = 10 + 4 * params['VehicleSpeed']
    
    # Solve wall temperature
    def wallEq_post(Tw):
        return (hin[i] * params['Di'] * (Tgas[i] - Tw) 
                - hout[i] * params['Do'] * (Tw - params['Tamb']) 
                - params['eps'] * params['sigma'] * params['Do'] * (Tw**4 - params['Tamb']**4))
        
    Twall[i] = fsolve(wallEq_post, x0=Tgas[i] - 30)[0]
    
    # Heat loss per unit length
    qConv[i] = hout[i] * np.pi * params['Do'] * (Twall[i] - params['Tamb'])
    qRad[i] = params['eps'] * params['sigma'] * np.pi * params['Do'] * (Twall[i]**4 - params['Tamb']**4)
    qTotal[i] = qConv[i] + qRad[i]

# ====================== SUMMARY ==========================================
print("\n================== RESULTS ==================")
print(f"Outlet Gas Temperature      : {Tgas[-1] - 273.15:.2f} C")
print(f"Outlet Wall Temperature     : {Twall[-1] - 273.15:.2f} C")
print(f"Average Reynolds Number     : {np.mean(Re):.0f}")
print(f"Average Nusselt Number      : {np.mean(Nu):.1f}")
print(f"Average Internal h          : {np.mean(hin):.1f} W/m2K")
print(f"Average External h          : {np.mean(hout):.1f} W/m2K")
print(f"Total Heat Loss             : {np.trapezoid(y=qTotal, x=x):.1f} W")
print(f"Convection Heat Loss        : {np.trapezoid(y=qConv, x=x):.1f} W")
print(f"Radiation Heat Loss         : {np.trapezoid(y=qRad, x=x):.1f} W")
print("=============================================\n")

# ========================== PLOTS ========================================
plt.figure(figsize=(8, 5))
plt.plot(x, Tgas - 273.15, 'b', linewidth=2, label='Gas')
plt.grid(True)
plt.xlabel('Distance Along Pipe (m)')
plt.ylabel('Temperature (°C)')
plt.title('Gas Temperature')
plt.legend()

plt.figure(figsize=(8, 5))
plt.plot(x, hin, linewidth=2)
plt.grid(True)
plt.xlabel('Distance (m)')
plt.ylabel('h_in (W/m^2K)')
plt.title('Internal Heat Transfer Coefficient')

plt.figure(figsize=(8, 5))
plt.plot(x, Re, linewidth=2)
plt.grid(True)
plt.xlabel('Distance (m)')
plt.ylabel('Reynolds Number')
plt.title('Reynolds Number')

plt.figure(figsize=(8, 5))
plt.plot(x, Nu, linewidth=2)
plt.grid(True)
plt.xlabel('Distance (m)')
plt.ylabel('Nusselt Number')
plt.title('Nusselt Number')

plt.figure(figsize=(8, 5))
plt.plot(x, Vel, linewidth=2)
plt.grid(True)
plt.xlabel('Distance (m)')
plt.ylabel('Velocity (m/s)')
plt.title('Gas Velocity')

plt.figure(figsize=(8, 5))
plt.plot(x, qConv, 'b', linewidth=2, label='Convection')
plt.plot(x, qRad, 'r', linewidth=2, label='Radiation')
plt.plot(x, qTotal, 'k', linewidth=2, label='Total')
plt.grid(True)
plt.xlabel('Distance (m)')
plt.ylabel('Heat Loss (W/m)')
plt.title('Heat Loss per Unit Length')
plt.legend()

plt.show()