program steady_heat_conduction_cylinders_spheres
    implicit none
    
    ! Variable declarations
    integer :: n_layers, i, n_nodes, j, geometry_type
    real(8) :: T_inner, T_outer, Q, R_total, T_interface
    real(8), allocatable :: r_inner(:), r_outer(:), k(:), R(:), T(:)
    real(8) :: r_val, dr_val, length_cyl, heat_flux, r_mid
    real(8) :: surface_inner, surface_outer
    character(len=50) :: filename
    character(len=20) :: geometry_name
    
    ! Read geometry type: 1 = Cylinder, 2 = Sphere
    read *, geometry_type
    
    if (geometry_type == 1) then
        geometry_name = "CYLINDER"
        read *, length_cyl  ! Length for cylinder
    else
        geometry_name = "SPHERE"
        length_cyl = 1.0d0  ! Not used for sphere
    end if
    
    ! Read number of layers
    read *, n_layers
    
    ! Allocate arrays
    allocate(r_inner(n_layers))
    allocate(r_outer(n_layers))
    allocate(k(n_layers))
    allocate(R(n_layers))
    allocate(T(0:n_layers))
    
    ! Read layer properties (inner radius, outer radius, conductivity)
    do i = 1, n_layers
        read *, r_inner(i)
        read *, r_outer(i)
        read *, k(i)
    end do
    
    ! Read boundary temperatures
    read *, T_inner
    read *, T_outer
    
    ! Calculate thermal resistances
    R_total = 0.0d0
    do i = 1, n_layers
        if (geometry_type == 1) then
            ! Cylindrical resistance: R = ln(r_o/r_i) / (2*pi*k*L)
            R(i) = log(r_outer(i) / r_inner(i)) / (2.0d0 * 3.14159265358979d0 * k(i) * length_cyl)
        else
            ! Spherical resistance: R = (1/r_i - 1/r_o) / (4*pi*k)
            R(i) = (1.0d0/r_inner(i) - 1.0d0/r_outer(i)) / (4.0d0 * 3.14159265358979d0 * k(i))
        end if
        R_total = R_total + R(i)
    end do
    
    ! Calculate heat transfer rate
    Q = (T_inner - T_outer) / R_total
    
    ! 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
    
    ! Calculate surface areas
    if (geometry_type == 1) then
        surface_inner = 2.0d0 * 3.14159265358979d0 * r_inner(1) * length_cyl
        surface_outer = 2.0d0 * 3.14159265358979d0 * r_outer(n_layers) * length_cyl
    else
        surface_inner = 4.0d0 * 3.14159265358979d0 * r_inner(1)**2
        surface_outer = 4.0d0 * 3.14159265358979d0 * r_outer(n_layers)**2
    end if
    
    ! Display results
    print *, '========================================='
    print *, 'STEADY STATE THERMAL CONDUCTION'
    print *, 'GEOMETRY: ', trim(geometry_name)
    print *, '========================================='
    print *, ''
    print *, 'INPUT PARAMETERS:'
    print *, '----------------------------------------'
    
    if (geometry_type == 1) then
        print '(A,F10.4,A)', ' Cylinder Length:           ', length_cyl, ' m'
        print '(A,F10.4,A)', ' Inner Radius:              ', r_inner(1), ' m'
        print '(A,F10.4,A)', ' Outer Radius:              ', r_outer(n_layers), ' m'
    else
        print '(A,F10.4,A)', ' Inner Radius:              ', r_inner(1), ' m'
        print '(A,F10.4,A)', ' Outer Radius:              ', r_outer(n_layers), ' m'
    end if
    
    print '(A,F10.2,A)', ' Inner Temperature:         ', T_inner, ' °C'
    print '(A,F10.2,A)', ' Outer Temperature:         ', T_outer, ' °C'
    print '(A,F10.2,A)', ' Temperature Difference:    ', (T_inner - T_outer), ' °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.4,A)', ' Inner Surface Area:          ', surface_inner, ' m²'
    print '(A,F12.4,A)', ' Outer Surface Area:          ', surface_outer, ' m²'
    print '(A,F12.2,A)', ' Inner Surface Heat Flux:     ', Q/surface_inner, ' W/m²'
    print '(A,F12.2,A)', ' Outer Surface Heat Flux:     ', Q/surface_outer, ' W/m²'
    print *, ''
    
    print *, '========================================='
    print *, 'LAYER ANALYSIS'
    print *, '========================================='
    print *, ''
    
    do i = 1, n_layers
        r_mid = (r_inner(i) + r_outer(i)) / 2.0d0
        
        print '(A,I2,A)', ' === Layer ', i, ' ==='
        print '(A,F10.4,A)', '   Inner Radius (r_i):      ', r_inner(i), ' m'
        print '(A,F10.4,A)', '                            ', r_inner(i)*1000, ' mm'
        print '(A,F10.4,A)', '   Outer Radius (r_o):      ', r_outer(i), ' m'
        print '(A,F10.4,A)', '                            ', r_outer(i)*1000, ' mm'
        print '(A,F10.4,A)', '   Thickness (r_o - r_i):   ', (r_outer(i)-r_inner(i)), ' m'
        print '(A,F10.4,A)', '                            ', (r_outer(i)-r_inner(i))*1000, ' 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, ' %'
        print '(A,F10.2,A)', '   Inner Temperature:       ', T(i-1), ' °C'
        print '(A,F10.2,A)', '   Outer Temperature:       ', T(i), ' °C'
        print '(A,F10.2,A)', '   Temperature Drop:        ', (T(i-1) - T(i)), ' °C'
        print *, ''
    end do
    
    print *, '========================================='
    print *, 'INTERFACE TEMPERATURES'
    print *, '========================================='
    print *, ''
    print '(A,F10.2,A)', ' Inner Surface (T0):        ', T(0), ' °C'
    do i = 1, n_layers - 1
       print '(A,I2,A,I2,A,F10.2,A)', ' Interface layer ', i, '-', i+1, ': ', T(i), ' °C'
    end do
    print '(A,F10.2,A)', ' Outer Surface (Tn):        ', T(n_layers), ' °C'
    print *, ''
    
    ! Generate temperature profile data file
    filename = 'temperature_profile.dat'
    open(unit=10, file=filename, status='replace')
    
    write(10, '(A)') '# Radius(m)  Temperature(°C)  Layer'
    
    do i = 1, n_layers
        n_nodes = 50
        dr_val = (r_outer(i) - r_inner(i)) / real(n_nodes, 8)
        
        do j = 0, n_nodes
            r_val = r_inner(i) + j * dr_val
            
            if (geometry_type == 1) then
                ! Cylindrical temperature distribution
                if (abs(r_outer(i) - r_inner(i)) > 1.0d-10) then
                    T_interface = T(i-1) - (T(i-1) - T(i)) * &
                                 log(r_val/r_inner(i)) / log(r_outer(i)/r_inner(i))
                else
                    T_interface = T(i-1)
                end if
            else
                ! Spherical temperature distribution
                if (abs(r_outer(i) - r_inner(i)) > 1.0d-10) then
                    T_interface = T(i-1) - (T(i-1) - T(i)) * &
                                 (1.0d0/r_inner(i) - 1.0d0/r_val) / &
                                 (1.0d0/r_inner(i) - 1.0d0/r_outer(i))
                else
                    T_interface = T(i-1)
                end if
            end if
            
            write(10, '(F12.6,2X,F12.4,2X,I3)') r_val, T_interface, i
        end do
    end do
    
    close(10)
    
    print *, '========================================='
    print *, 'DATA FILE'
    print *, '========================================='
    print *, ''
    print *, ' Temperature profile saved in:'
    print *, ' ', trim(filename)
    print *, ''
    print *, ' Format: Radius(m) Temperature(°C) Layer'
    print *, ''
    
    ! Additional performance metrics
    print *, '========================================='
    print *, 'PERFORMANCE METRICS'
    print *, '========================================='
    print *, ''
    
    ! Find the layer with maximum resistance (bottleneck)
    i = maxloc(R, dim=1)
    print '(A,I2)', ' Limiting Layer (max resistance): Layer ', i
    print '(A,F10.6,A)', ' Resistance of this layer:        ', R(i), ' K/W'
    print *, ''
    
    ! Calculate overall heat transfer coefficient
    print '(A,F10.4,A)', ' U-Coefficient (inner basis):     ', 1.0d0/(R_total*surface_inner), ' W/m²·K'
    print '(A,F10.4,A)', ' U-Coefficient (outer basis):     ', 1.0d0/(R_total*surface_outer), ' W/m²·K'
    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 *, ''
    
    ! Volume calculations
    if (geometry_type == 1) then
        print *, 'VOLUMETRIC CHARACTERISTICS:'
        print '(A,F12.6,A)', ' Inner Volume:      ', &
            3.14159265358979d0 * r_inner(1)**2 * length_cyl, ' m³'
        print '(A,F12.6,A)', ' Outer Volume:      ', &
            3.14159265358979d0 * r_outer(n_layers)**2 * length_cyl, ' m³'
        print '(A,F12.6,A)', ' Wall Volume:       ', &
            3.14159265358979d0 * (r_outer(n_layers)**2 - r_inner(1)**2) * length_cyl, ' m³'
    else
        print *, 'VOLUMETRIC CHARACTERISTICS:'
        print '(A,F12.6,A)', ' Inner Volume:      ', &
            (4.0d0/3.0d0) * 3.14159265358979d0 * r_inner(1)**3, ' m³'
        print '(A,F12.6,A)', ' Outer Volume:      ', &
            (4.0d0/3.0d0) * 3.14159265358979d0 * r_outer(n_layers)**3, ' m³'
        print '(A,F12.6,A)', ' Wall Volume:       ', &
            (4.0d0/3.0d0) * 3.14159265358979d0 * &
            (r_outer(n_layers)**3 - r_inner(1)**3), ' m³'
    end if
    print *, ''
    
    print *, '========================================='
    print *, 'END OF CALCULATION'
    print *, '========================================='
    
    ! Cleanup
    deallocate(r_inner, r_outer, k, R, T)
    
end program steady_heat_conduction_cylinders_spheres