Critical Radius of Insulation

program critical_radius_insulation
    implicit none
    
    ! Variable declarations
    integer :: geometry_type, n_points, i
    real(8) :: r_inner, k_insulation, h_outer, T_inner, T_ambient
    real(8) :: r_critical, length_cyl, Q_bare, Q_critical, Q_insulated
    real(8) :: r_test, delta_r, Q_test, R_total
    real(8) :: surface_inner, surface_outer, pi
    character(len=50) :: filename
    character(len=20) :: geometry_name
    
    pi = 3.14159265358979d0
    
    ! 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 input parameters
    read *, r_inner           ! Inner radius [m]
    read *, k_insulation      ! Thermal conductivity of insulation [W/m·K]
    read *, h_outer           ! Convection coefficient [W/m²·K]
    read *, T_inner           ! Inner surface temperature [°C]
    read *, T_ambient         ! Ambient temperature [°C]
    
    ! Calculate critical radius
    if (geometry_type == 1) then
        ! Cylinder: r_cr = k/h
        r_critical = k_insulation / h_outer
    else
        ! Sphere: r_cr = 2k/h
        r_critical = 2.0d0 * k_insulation / h_outer
    end if
    
    ! Calculate heat transfer for bare pipe
    if (geometry_type == 1) then
        surface_inner = 2.0d0 * pi * r_inner * length_cyl
        Q_bare = h_outer * surface_inner * (T_inner - T_ambient)
    else
        surface_inner = 4.0d0 * pi * r_inner**2
        Q_bare = h_outer * surface_inner * (T_inner - T_ambient)
    end if
    
    ! Calculate heat transfer at critical radius
    if (geometry_type == 1) then
        R_total = log(r_critical/r_inner)/(2.0d0*pi*k_insulation*length_cyl) + &
                  1.0d0/(h_outer*2.0d0*pi*r_critical*length_cyl)
    else
        R_total = (1.0d0/r_inner - 1.0d0/r_critical)/(4.0d0*pi*k_insulation) + &
                  1.0d0/(h_outer*4.0d0*pi*r_critical**2)
    end if
    Q_critical = (T_inner - T_ambient) / R_total
    
    ! Display results
    print *, '========================================='
    print *, 'CRITICAL RADIUS OF INSULATION'
    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'
    end if
    
    print '(A,F10.4,A)', ' Inner Radius (r_i):          ', r_inner, ' m'
    print '(A,F10.4,A)', '                              ', r_inner*1000, ' mm'
    print '(A,F10.4,A)', ' Insulation Conductivity (k): ', k_insulation, ' W/m·K'
    print '(A,F10.2,A)', ' Convection Coefficient (h):  ', h_outer, ' W/m²·K'
    print '(A,F10.2,A)', ' Surface Temperature (T_s):   ', T_inner, ' °C'
    print '(A,F10.2,A)', ' Ambient Temperature (T_∞):   ', T_ambient, ' °C'
    print '(A,F10.2,A)', ' Temperature Difference:      ', (T_inner - T_ambient), ' °C'
    print *, ''
    
    print *, '========================================='
    print *, 'CRITICAL RADIUS'
    print *, '========================================='
    print *, ''
    
    if (geometry_type == 1) then
        print *, ' Formula (Cylinder): r_cr = k/h'
    else
        print *, ' Formula (Sphere):   r_cr = 2k/h'
    end if
    
    print *, ''
    print '(A,F10.6,A)', ' Critical Radius (r_cr):      ', r_critical, ' m'
    print '(A,F10.4,A)', '                              ', r_critical*1000, ' mm'
    print '(A,F10.4,A)', ' Critical Thickness:          ', (r_critical - r_inner), ' m'
    print '(A,F10.4,A)', '                              ', (r_critical - r_inner)*1000, ' mm'
    print *, ''
    
    ! Analysis
    print *, '========================================='
    print *, 'ANALYSIS'
    print *, '========================================='
    print *, ''
    
    if (r_inner < r_critical) then
        print *, ' WARNING: r_i < r_cr'
        print *, ' ----------------------------------------'
        print *, ' Adding insulation INCREASES heat losses'
        print *, ' up to the critical radius!'
        print *, ''
        print *, ' Recommendations:'
        print *, ' - If insulation needed: use'
        print *, '   r_o > r_cr to reduce losses'
        print *, ' - Or increase h (forced ventilation)'
        print *, ' - Or use better insulation'
        print *, '   (lower k)'
    else if (abs(r_inner - r_critical) < 0.001d0) then
        print *, ' LIMITING CASE: r_i ≈ r_cr'
        print *, ' ----------------------------------------'
        print *, ' The pipe is near the critical radius.'
        print *, ' Any insulation will be beneficial.'
    else
        print *, ' FAVORABLE: r_i > r_cr'
        print *, ' ----------------------------------------'
        print *, ' Any insulation thickness reduces'
        print *, ' heat losses.'
    end if
    print *, ''
    
    print *, '========================================='
    print *, 'HEAT TRANSFER'
    print *, '========================================='
    print *, ''
    print '(A,F12.2,A)', ' Q without insulation:        ', Q_bare, ' W'
    print '(A,F12.2,A)', ' Q at critical radius:        ', Q_critical, ' W'
    
    if (Q_critical > Q_bare) then
        print '(A,F10.2,A)', ' Increase at r_cr:            ', &
                           ((Q_critical-Q_bare)/Q_bare)*100, ' %'
    else
        print '(A,F10.2,A)', ' Reduction at r_cr:           ', &
                           ((Q_bare-Q_critical)/Q_bare)*100, ' %'
    end if
    print *, ''
    
    ! Generate data file for various insulation thicknesses
    filename = 'critical_radius_analysis.dat'
    open(unit=10, file=filename, status='replace')
    
    write(10, '(A)') '# Outer_Radius(m)  Thickness(m)  Q(W)  Q/Q_bare  Resistance(K/W)'
    
    n_points = 100
    delta_r = (3.0d0 * r_critical) / real(n_points, 8)
    
    do i = 0, n_points
        r_test = r_inner + i * delta_r
        
        if (geometry_type == 1) then
            if (r_test > r_inner) then
                R_total = log(r_test/r_inner)/(2.0d0*pi*k_insulation*length_cyl) + &
                         1.0d0/(h_outer*2.0d0*pi*r_test*length_cyl)
            else
                R_total = 1.0d0/(h_outer*2.0d0*pi*r_inner*length_cyl)
            end if
        else
            if (r_test > r_inner) then
                R_total = (1.0d0/r_inner - 1.0d0/r_test)/(4.0d0*pi*k_insulation) + &
                         1.0d0/(h_outer*4.0d0*pi*r_test**2)
            else
                R_total = 1.0d0/(h_outer*4.0d0*pi*r_inner**2)
            end if
        end if
        
        Q_test = (T_inner - T_ambient) / R_total
        
        write(10, '(F12.6,2X,F12.6,2X,F12.4,2X,F12.6,2X,F12.6)') &
              r_test, (r_test - r_inner), Q_test, Q_test/Q_bare, R_total
    end do
    
    close(10)
    
    print *, '========================================='
    print *, 'PRACTICAL RECOMMENDATIONS'
    print *, '========================================='
    print *, ''
    
    if (r_inner < r_critical) then
        print '(A,F10.4,A)', ' Minimum recommended thickness:  ', &
                           (r_critical*1.5 - r_inner)*1000, ' mm'
        print *, ' (r_o = 1.5 × r_cr for effective reduction)'
    else
        print *, ' Any insulation thickness is beneficial.'
        print *, ' Choose based on economic constraints.'
    end if
    print *, ''
    
    print *, '========================================='
    print *, 'DATA FILE'
    print *, '========================================='
    print *, ''
    print *, ' Complete analysis saved in:'
    print *, ' ', trim(filename)
    print *, ''
    print *, ' Format: Outer_Radius Thickness Q Q/Q_bare R_total'
    print *, ''
    
    print *, '========================================='
    print *, 'END OF CALCULATION'
    print *, '========================================='
    
end program critical_radius_insulation