program composite_wall_convection
    implicit none

    ! Parameters
    real, parameter :: sigma = 5.67037e-8  ! Stefan-Boltzmann constant, W/(m^2 K^4)

    ! Inputs
    real :: Th, Tc          ! Inside and outside fluid temperatures, C
    real :: hi, ho          ! Inside and outside convection coefficients, W/(m^2 K)
    integer :: use_rad      ! 0 = No radiation, 1 = Yes radiation
    real :: eps_i, eps_o    ! Inside and outside surface emissivities
    real :: A               ! Wall area, m^2
    integer :: N            ! Number of layers
    real, allocatable :: L(:), k(:), R(:) ! Thickness (mm), conductivity (W/m-K), resistance (m^2-K/W)
    real :: RH_in           ! Inside room relative humidity, %

    ! Solved variables
    real :: L_m             ! Layer thickness in meters
    real :: R_wall          ! Total conduction resistance of the wall
    real :: R_si, R_se      ! Effective surface resistances
    real :: R_total         ! Overall thermal resistance, m^2 K/W
    real :: U               ! Overall heat transfer coefficient, W/m^2 K
    real :: q               ! Heat flux, W/m^2
    real :: Q_total         ! Total heat transfer rate, W
    real :: T_si, T_so      ! Inner and outer surface temperatures, C
    real, allocatable :: T_node(:) ! Node temperatures at interfaces
    real :: Tdew            ! Dew point temperature, C
    
    ! Iterative variables
    real :: Tsi_K, Tso_K, Th_K, Tc_K
    real :: hr_i, hr_o
    real :: T_si_old, T_so_old
    integer :: iter
    logical :: converged
    
    ! UI display arrays
    real, allocatable :: R_pct(:) ! Percentage resistance contribution
    
    ! Loop variables
    integer :: i
    real :: alpha_val

    ! 1. Read Inputs from stdin
    read(*,*) Th
    read(*,*) Tc
    read(*,*) hi
    read(*,*) ho
    read(*,*) use_rad
    if (use_rad == 1) then
        read(*,*) eps_i
        read(*,*) eps_o
    else
        eps_i = 0.0
        eps_o = 0.0
    end if
    read(*,*) A
    read(*,*) N

    allocate(L(N))
    allocate(k(N))
    allocate(R(N))
    allocate(R_pct(N))
    allocate(T_node(N+1))

    do i = 1, N
        read(*,*) L(i)
        read(*,*) k(i)
    end do
    read(*,*) RH_in

    ! 2. Compute Wall Conduction Resistances
    R_wall = 0.0
    do i = 1, N
        L_m = L(i) / 1000.0  ! convert mm to meters
        R(i) = L_m / k(i)
        R_wall = R_wall + R(i)
    end do

    ! 3. Iterative solver for Surface Temperatures (to resolve radiation non-linearity)
    T_si = Th - 0.1 * (Th - Tc)
    T_so = Tc + 0.1 * (Th - Tc)
    
    iter = 0
    converged = .false.
    do while (.not. converged .and. iter < 100)
        iter = iter + 1
        T_si_old = T_si
        T_so_old = T_so
        
        if (use_rad == 1) then
            ! Radiation requires temperatures in Kelvin
            Th_K = Th + 273.15
            Tc_K = Tc + 273.15
            Tsi_K = T_si + 273.15
            Tso_K = T_so + 273.15
            
            hr_i = eps_i * sigma * (Th_K + Tsi_K) * (Th_K**2 + Tsi_K**2)
            hr_o = eps_o * sigma * (Tso_K + Tc_K) * (Tso_K**2 + Tc_K**2)
        else
            hr_i = 0.0
            hr_o = 0.0
        end if
        
        R_si = 1.0 / (hi + hr_i)
        R_se = 1.0 / (ho + hr_o)
        
        R_total = R_si + R_wall + R_se
        q = (Th - Tc) / R_total
        
        T_si = Th - q * R_si
        T_so = Tc + q * R_se
        
        if (abs(T_si - T_si_old) < 1e-5 .and. abs(T_so - T_so_old) < 1e-5) then
            converged = .true.
        end if
    end do

    U = 1.0 / R_total
    Q_total = q * A

    ! 4. Compute Interface Temperatures
    T_node(1) = T_si
    do i = 1, N
        T_node(i+1) = T_node(i) - q * (L(i) / 1000.0) / k(i)
    end do

    ! Calculate percentage resistance of each component
    do i = 1, N
        R_pct(i) = (R(i) / R_total) * 100.0
    end do

    ! 5. Dew Point Calculation (Magnus-Tetens)
    alpha_val = (17.625 * Th) / (243.04 + Th) + log(RH_in / 100.0)
    Tdew = (243.04 * alpha_val) / (17.625 - alpha_val)

    ! 6. Output Engineering Report
    write(*,*) "=========================================================================="
    write(*,*) "         COMPOSITE WALL THERMAL ANALYSIS REPORT (CONVECTION & RAD)        "
    write(*,*) "=========================================================================="
    write(*, '(A, F7.2, A, F7.2, A)') " Boundary Temperatures:  Inside = ", Th, " C, Outside = ", Tc, " C"
    write(*, '(A, F7.2, A, F7.2, A)') " Convection Coeffs:      Inside = ", hi, " W/m2-K, Outside = ", ho, " W/m2-K"
    if (use_rad == 1) then
        write(*, '(A, F5.3, A, F5.3)') " Surface Emissivities:   Inside = ", eps_i, ", Outside = ", eps_o
        write(*, '(A, F7.3, A, F7.3, A)') " Radiative Coeffs (hr):  Inside = ", hr_i, ", Outside = ", hr_o, " W/m2-K"
    else
        write(*,*) " Surface Radiation:      Disabled"
    end if
    write(*, '(A, F8.2, A)') " Wall Surface Area:     ", A, " m2"
    write(*, '(A, I3)') " Number of Wall Layers:  ", N

    write(*,*) "--------------------------------------------------------------------------"
    write(*,*) " LAYER RESISTANCE & PROPERTIES"
    write(*,*) "--------------------------------------------------------------------------"
    write(*,*) " Layer  Thick(mm)  Conductivity(W/m-K)  Resistance(m2-K/W)  Resistance %"
    do i = 1, N
        write(*, '(I4.2, F12.1, F21.3, F20.4, F13.1, A)') &
            i, L(i), k(i), R(i), R_pct(i), " %"
    end do
    write(*, '(A, F20.4, F13.1, A)') " Inside Surface Resistance (R_si):   ", R_si, (R_si / R_total * 100.0), " %"
    write(*, '(A, F20.4, F13.1, A)') " Outside Surface Resistance (R_se):  ", R_se, (R_se / R_total * 100.0), " %"

    write(*,*) "--------------------------------------------------------------------------"
    write(*,*) " THERMAL PERFORMANCE METRICS"
    write(*,*) "--------------------------------------------------------------------------"
    write(*, '(A, F10.5, A)') " Total Thermal Resistance (R_total):   ", R_total, " m2-K/W"
    write(*, '(A, F10.5, A)') " Overall Heat Transfer Coeff (U):       ", U, " W/m2-K"
    write(*, '(A, F10.2, A)') " Heat Flux (q):                         ", q, " W/m2"
    write(*, '(A, F10.2, A)') " Total Heat Loss/Gain (Q):              ", Q_total, " W"

    write(*,*) "--------------------------------------------------------------------------"
    write(*,*) " INTERFACE TEMPERATURE PROFILE"
    write(*,*) "--------------------------------------------------------------------------"
    write(*, '(A, F7.2, A)') " Inner Surface Temp (T_si):           ", T_si, " C"
    do i = 1, N
        write(*, '(A, I2, A, I2, A, F7.2, A)') " Interface ", i, " Temp (T_", i, "):             ", T_node(i+1), " C"
    end do
    write(*, '(A, F7.2, A)') " Outer Surface Temp (T_so):           ", T_so, " C"

    write(*,*) "--------------------------------------------------------------------------"
    write(*,*) " CONDENSATION RISK ANALYSIS"
    write(*,*) "--------------------------------------------------------------------------"
    write(*, '(A, F5.1, A, F7.2, A)') " Inside Air RH: ", RH_in, " %, Dew Point (Tdew): ", Tdew, " C"
    write(*, '(A, F7.2, A)') " Inner Surface Temperature (T_si):    ", T_si, " C"
    if (T_si <= Tdew) then
        write(*,*) " WARNING: Surface temperature is below dew point! Condensation risk."
    else
        write(*,*) " Status: Surface temperature is safe. No condensation risk."
    end if
    write(*,*) "=========================================================================="

    deallocate(L)
    deallocate(k)
    deallocate(R)
    deallocate(R_pct)
    deallocate(T_node)
end program composite_wall_convection