import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# """
# Steel Heating Simulation according to EN 1993-1-2 (Eurocode 3: Fire Design)
#
# This script implements a simplified thermal model for steel elements
# exposed to fire, following the methodology of EN 1993-1-2.
#
# ---------------------------------------------------------------
# 1. Governing Equation (EN 1993-1-2, Eq. 4.3)
# ---------------------------------------------------------------
# The temperature rise of steel θ_a [°C] is determined from the heat balance:
#
#     dθ_a/dt = (k_sh * (A_m / V) * (q_net)) / (ρ_a * c_a)
#
# where:
# - A_m / V  : section factor [m⁻¹] (heated surface area / steel volume)
# - ρ_a      : density of steel [kg/m³] (≈ 7850)
# - c_a      : specific heat of steel [J/kg·K], temperature-dependent
# - k_sh     : correction factor for shadow effects (≈ 1.0 if ignored)
# - q_net    : net heat flux [W/m²] from convection + radiation
#
# ---------------------------------------------------------------
# 2. Heat Transfer Mechanisms
# ---------------------------------------------------------------
# Total heat flux is the sum of:
# - Convection (EN 1993-1-2 §4.2.5.1):
#       q_conv = h_c · (θ_g - θ_a)
#   where h_c ≈ 25 W/m²·K is the convective heat transfer coefficient.
#
# - Radiation (EN 1993-1-2 §4.2.5.2):
#       q_rad = φ · ε_m · σ · (Θ_g⁴ - Θ_a⁴)
#   where:
#     φ   = configuration factor (≈ 1.0 for uniform exposure),
#     ε_m = surface emissivity (≈ 0.7),
#     σ   = Stefan-Boltzmann constant (5.67×10⁻⁸ W/m²·K⁴),
#     Θ   = absolute temperature in Kelvin.
#
# ---------------------------------------------------------------
# 3. Gas Temperature (θ_g)
# ---------------------------------------------------------------
# The gas temperature-time curve is given by the *parametric fire curve*
# (EN 1991-1-2 Annex A). In this example, a simplified expression is used:
# External fire curve
#
#     θ_gas(t) = 660 · (1 - 0.687·e^(-0.32·t) - 0.313·e^(-3.8·t)) + 20
#
# where t is time in minutes.
#
# ---------------------------------------------------------------
# 4. Specific Heat of Steel (c_a)
# ---------------------------------------------------------------
# Temperature-dependent values are taken from EN 1993-1-2 Table 3.2:
#
# - θ < 600°C: polynomial fit
# - 600 ≤ θ < 735°C: 666 + 13002/(738 - θ)
# - 735 ≤ θ < 900°C: 545 + 17820/(θ - 731)
# - θ ≥ 900°C: constant 650 J/kg·K
#
# ---------------------------------------------------------------
# 5. Outputs
# ---------------------------------------------------------------
# - Animated visualization of steel temperature rise
# - Comparison with the applied gas temperature
# - Temperature distribution on a representative steel plate
#
# This model provides an illustrative calculation of steel heating
# under fire exposure in accordance with Eurocode 3 fire design rules.

# ************************************************************************************************************
# -------------------------
# Steel and fire parameters
# -------------------------
rho_a = 7850.0  # density of steel [kg/m3]

def c_steel(theta):
    if theta < 600:
        return 425 + 0.773 * theta - 1.69e-3 * theta ** 2 + 2.22e-6 * theta ** 3
    elif theta < 735:
        return 666 + 13002 / (738 - theta)
    elif theta < 900:
        return 545 + 17820 / (theta - 731)
    else:
        return 650.0

def theta_gas(t_min):
    return 660 * (1 - 0.687 * np.exp(-0.32 * t_min) - 0.313 * np.exp(-3.8 * t_min)) + 20

# -------------------------
# Input parameters
# -------------------------
Am_over_V = 146.36
h_c = 25.0
epsilon_m = 0.7
phi = 1.0
sigma = 5.67e-8
k_sh = 1.0

Delta_t = 5      # second
t_max = 900      # second = 15
# -------------------------
# Temperature calculation looping
# -------------------------
times = np.arange(0, t_max + Delta_t, Delta_t)
steel_temp = [20]

for i in range(1, len(times)):
    t_sec = times[i]
    t_min = t_sec / 60.0
    theta_g = theta_gas(t_min)
    theta_a = steel_temp[-1]

    c_a = c_steel(theta_a)
    Theta_g = theta_g + 273.15
    Theta_a = theta_a + 273.15

    q_conv = h_c * (theta_g - theta_a)
    q_rad = phi * epsilon_m * sigma * (Theta_g ** 4 - Theta_a ** 4)
    q_net = q_conv + q_rad

    dtheta_dt = (k_sh * Am_over_V * q_net) / (rho_a * c_a)
    theta_new = theta_a + dtheta_dt * Delta_t
    steel_temp.append(theta_new)

# -------------------------
# Animation setup
# -------------------------
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Left: plate
plate = ax1.imshow(np.ones((10, 15)) * steel_temp[0],
                   cmap="rainbow", vmin=20, vmax=max(steel_temp),
                   extent=[0, 262, 0, 146], origin="lower")
cbar = plt.colorbar(plate, ax=ax1)
cbar.set_label("Temperature [°C]")
ax1.set_title("Steel Plate Heating")
ax1.set_xlabel("Length [mm]")
ax1.set_ylabel("Width [mm]")
time_text = ax1.text(5, 170, "", color="black", fontsize=12,
                     bbox=dict(facecolor="white", alpha=0.7))


# Right: incremental temperature-time graph
ax2.set_xlim(0, t_max / 60)
ax2.set_ylim(min(steel_temp) - 10, max(steel_temp) + 100)
ax2.set_xlabel("Time [min]")
ax2.set_ylabel("Temperature [°C]")
ax2.set_title("Temperature History")
ax2.grid(True)

# plot steel curve (incremental result)
line, = ax2.plot([], [], color="red", lw=2, label="Steel temperature")

# plot gas temperature curve
gas_line, = ax2.plot(times / 60, [theta_gas(t / 60) for t in times],
                     color="blue", lw=1, linestyle="--", label="Gas temperature")

# add legend
ax2.legend(loc='lower right')


def update(frame):
    # Update plate color
    temp = steel_temp[frame]
    plate.set_data(np.ones((10, 15)) * temp)
    time_text.set_text(f"t = {times[frame] / 60:.1f} min\nT = {temp:.0f} °C")

    # Update incremental line
    line.set_data(times[:frame + 1] / 60, steel_temp[:frame + 1])

    return plate, time_text, line


ani = animation.FuncAnimation(fig, update, frames=len(times),
                              interval=20, blit=False, repeat=False)

# -------------------------
# Save as GIF
# -------------------------
ani.save("steel_heating.gif", writer="pillow", fps=10)

plt.tight_layout()
plt.show()

# exposed_sides = 2        # 1 = one side, 2 = both sides
# Am_over_V = exposed_sides / t

