program steady_heat_conduction_plane_walls
    implicit none
    
    integer :: n_layers, i
    real(8) :: T_inner, T_outer, Q, R_total, area, heat_flux
    real(8), allocatable :: thickness(:), k(:), R(:), T(:)
    
    ! Read number of layers
    read *, n_layers
    
    ! Allocate arrays
    allocate(thickness(n_layers))
    allocate(k(n_layers))
    allocate(R(n_layers))
    allocate(T(0:n_layers))
    
    ! Read layer properties (thickness, conductivity)
    do i = 1, n_layers
        read *, thickness(i)
        read *, k(i)
    end do
    
    ! Read boundary temperatures and area
    read *, T_inner
    read *, T_outer
    read *, area
    
    ! Calculate thermal resistances in K/W
    R_total = 0.0d0
    do i = 1, n_layers
        ! R_i = L_i / (k_i * A)
        R(i) = thickness(i) / (k(i) * area)
        R_total = R_total + R(i)
    end do
    
    ! Calculate heat transfer rate Q (W)
    Q = (T_inner - T_outer) / R_total
    heat_flux = Q / area
    
    ! Calculate interface temperatures
    T(0) = T_inner
    T(n_layers) = T_outer
    do i = 1, n_layers - 1
        T(i) = T_inner - (T_inner - T_outer) * sum(R(1:i)) / R_total
    end do
    
    ! Display results
    print *, '========================================='
    print *, 'STEADY STATE THERMAL CONDUCTION'
    print *, 'GEOMETRY: PLANE WALLS'
    print *, '========================================='
    print *, ''
    print *, 'INPUT PARAMETERS:'
    print *, '----------------------------------------'
    print '(A,F10.4,A)', ' Wall Area:                 ', area, ' m2'
    print '(A,F10.2,A)', ' Inner Temperature:         ', T_inner, ' deg-C'
    print '(A,F10.2,A)', ' Outer Temperature:         ', T_outer, ' deg-C'
    print '(A,F10.2,A)', ' Temperature Difference:    ', (T_inner - T_outer), ' deg-C'
    print '(A,I3)', ' Number of Layers:          ', n_layers
    print *, ''
    
    print *, '========================================='
    print *, 'MAIN RESULTS'
    print *, '========================================='
    print *, ''
    print '(A,F12.6,A)', ' Total Thermal Resistance:    ', R_total, ' K/W'
    print '(A,F12.2,A)', ' Total Heat Flux (Q):         ', Q, ' W'
    print '(A,F12.2,A)', ' Heat Flux Per Area (q):      ', heat_flux, ' W/m2'
    print *, ''
    
    print *, '========================================='
    print *, 'LAYER ANALYSIS'
    print *, '========================================='
    print *, ''
    
    do i = 1, n_layers
        print '(A,I2,A)', ' === Layer ', i, ' ==='
        print '(A,F10.4,A)', '   Thickness (L):           ', thickness(i), ' m'
        print '(A,F10.4,A)', '                            ', thickness(i)*1000.0d0, ' mm'
        print '(A,F10.4,A)', '   Conductivity (k):        ', k(i), ' W/m.K'
        print '(A,F10.6,A)', '   Resistance (R):          ', R(i), ' K/W'
        print '(A,F8.2,A)', '   % of Total Resistance:   ', (R(i)/R_total)*100.0d0, ' %'
        print '(A,F10.2,A)', '   Inner Temperature:       ', T(i-1), ' deg-C'
        print '(A,F10.2,A)', '   Outer Temperature:       ', T(i), ' deg-C'
        print '(A,F10.2,A)', '   Temperature Drop:        ', (T(i-1) - T(i)), ' deg-C'
        print *, ''
    end do
    
    print *, '========================================='
    print *, 'INTERFACE TEMPERATURES'
    print *, '========================================='
    print *, ''
    print '(A,F10.2,A)', ' Inner Surface (T0):        ', T(0), ' deg-C'
    do i = 1, n_layers - 1
       print '(A,I2,A,I2,A,F10.2,A)', ' Interface layer ', i, '-', i+1, ': ', T(i), ' deg-C'
    end do
    print '(A,F10.2,A)', ' Outer Surface (Tn):        ', T(n_layers), ' deg-C'
    print *, ''
    
    print *, '========================================='
    print *, 'PERFORMANCE METRICS'
    print *, '========================================='
    print *, ''
    
    ! Energy loss calculations
    print *, 'ENERGY LOSSES (estimates):'
    print '(A,F12.2,A)', ' Per hour:    ', Q * 3600.0d0 / 1000.0d0, ' kJ/h'
    print '(A,F12.2,A)', ' Per day:     ', Q * 86400.0d0 / 1000.0d0, ' kJ/day'
    print '(A,F12.2,A)', '              ', Q * 24.0d0 / 1000.0d0, ' kWh/day'
    print '(A,F12.2,A)', ' Per year:    ', Q * 8760.0d0 / 1000.0d0, ' kWh/year'
    print *, ''
    
    print *, '========================================='
    print *, 'END OF CALCULATION'
    print *, '========================================='
    
    deallocate(thickness, k, R, T)
end program steady_heat_conduction_plane_walls