program Rectangulaire_Fin_Heat_Conduction
    implicit none
    
    ! Variable declarations
    real(8) :: L, t, w, k, h, T0, Tinf
    real(8) :: P, A, m, mL, Q, eta, epsilon, A_surface, A_total
    real(8) :: theta_base, x, theta_x, T_x
    integer :: i, n_points
    
    ! Read input parameters
    write(*,*) '=============================================='
    write(*,*) '   RECTANGULAR FIN HEAT CONDUCTION ANALYSIS'
    write(*,*) '=============================================='
    write(*,*)
    
    write(*,*) 'Reading input parameters...'
    read(*,*) L          ! Fin length [m]
    read(*,*) t          ! Fin thickness [m]
    read(*,*) w          ! Fin width [m]
    read(*,*) k          ! Thermal conductivity [W/m·K]
    read(*,*) h          ! Convection coefficient [W/m²·K]
    read(*,*) T0         ! Base temperature [°C]
    read(*,*) Tinf       ! Ambient temperature [°C]
    
    write(*,*)
    write(*,*) '=============================================='
    write(*,*) '              INPUT PARAMETERS'
    write(*,*) '=============================================='
    write(*,*)
    write(*,'(A,F10.6,A)') ' Fin Length (L):              ', L, ' m'
    write(*,'(A,F10.6,A)') ' Fin Thickness (t):           ', t, ' m'
    write(*,'(A,F10.6,A)') ' Fin Width (w):               ', w, ' m'
    write(*,'(A,F10.2,A)') ' Thermal Conductivity (k):    ', k, ' W/m·K'
    write(*,'(A,F10.2,A)') ' Convection Coefficient (h):  ', h, ' W/m²·K'
    write(*,'(A,F10.2,A)') ' Base Temperature (T₀):       ', T0, ' °C'
    write(*,'(A,F10.2,A)') ' Ambient Temperature (T∞):    ', Tinf, ' °C'
    write(*,*)
    
    ! Calculate geometric parameters
    P = 2.0d0 * (w + t)              ! Perimeter [m]
    A = w * t                         ! Cross-sectional area [m²]
    A_surface = 2.0d0 * w * L         ! Surface area (top and bottom) [m²]
    A_total = A_surface + 2.0d0 * t * L  ! Total surface area including edges [m²]
    
    write(*,*) '=============================================='
    write(*,*) '         GEOMETRIC CALCULATIONS'
    write(*,*) '=============================================='
    write(*,*)
    write(*,'(A,F12.6,A)') ' Perimeter (P):               ', P, ' m'
    write(*,'(A,F12.9,A)') ' Cross-sectional Area (A):    ', A, ' m²'
    write(*,'(A,F12.6,A)') ' Surface Area (A_surface):    ', A_surface, ' m²'
    write(*,'(A,F12.6,A)') ' Total Surface Area:          ', A_total, ' m²'
    write(*,*)
    
    ! Calculate fin parameter
    m = sqrt((h * P) / (k * A))       ! Fin parameter [1/m]
    mL = m * L                         ! Dimensionless fin length
    
    write(*,*) '=============================================='
    write(*,*) '           FIN PARAMETERS'
    write(*,*) '=============================================='
    write(*,*)
    write(*,'(A,F12.4,A)') ' Fin Parameter (m):           ', m, ' 1/m'
    write(*,'(A,F12.4)')   ' Dimensionless Length (mL):   ', mL
    write(*,*)
    
    ! Calculate heat transfer rate for infinitely long fin approximation
    theta_base = T0 - Tinf
    
    if (mL > 3.0d0) then
        ! For large mL, use infinitely long fin approximation
        Q = sqrt(h * P * k * A) * theta_base
        eta = 1.0d0 / mL
        write(*,*) 'Note: Using infinitely long fin approximation (mL > 3)'
    else
        ! For finite length fin with adiabatic tip
        Q = sqrt(h * P * k * A) * theta_base * tanh(mL)
        eta = tanh(mL) / mL
    end if
    
    ! Calculate fin efficiency and effectiveness
    epsilon = Q / (h * A_surface * theta_base)
    
    write(*,*) '=============================================='
    write(*,*) '           PERFORMANCE RESULTS'
    write(*,*) '=============================================='
    write(*,*)
    write(*,'(A,F12.4,A)') ' Heat Transfer Rate (Q):      ', Q, ' W'
    write(*,'(A,F10.4,A)') ' Fin Efficiency (η):          ', eta * 100.0d0, ' %'
    write(*,'(A,F10.4)')   ' Fin Effectiveness (ε):       ', epsilon
    write(*,*)
    
    if (epsilon < 2.0d0) then
        write(*,*) 'WARNING: Fin effectiveness < 2 (fin may not be efficient)'
    else if (epsilon > 10.0d0) then
        write(*,*) 'EXCELLENT: Fin effectiveness > 10 (highly efficient fin)'
    else
        write(*,*) 'GOOD: Fin effectiveness is acceptable'
    end if
    write(*,*)
    
    ! Calculate temperature distribution along the fin
    write(*,*) '=============================================='
    write(*,*) '      TEMPERATURE DISTRIBUTION'
    write(*,*) '=============================================='
    write(*,*)
    write(*,*) '  Position (m)  |  Temperature (°C)  |  θ/θ₀'
    write(*,*) '----------------------------------------------'
    
    n_points = 20
    do i = 0, n_points
        x = (dble(i) / dble(n_points)) * L
        
        ! Temperature distribution for adiabatic tip
        theta_x = theta_base * (cosh(m * (L - x)) / cosh(mL))
        T_x = Tinf + theta_x
        
        write(*,'(2X,F10.6,4X,A,2X,F10.3,8X,A,2X,F8.5)') &
            x, '|', T_x, '|', theta_x / theta_base
    end do
    
    write(*,*)
    write(*,*) '=============================================='
    write(*,*) '           HEAT FLUX ANALYSIS'
    write(*,*) '=============================================='
    write(*,*)
    
    ! Heat flux at the base
    write(*,'(A,F12.4,A)') ' Heat flux at base (q₀):      ', Q / A, ' W/m²'
    write(*,'(A,F12.4,A)') ' Average heat flux:           ', Q / A_surface, ' W/m²'
    write(*,*)
    
    ! Tip temperature
    theta_x = theta_base / cosh(mL)
    T_x = Tinf + theta_x
    
    write(*,*) '=============================================='
    write(*,*) '         TEMPERATURE EXTREMES'
    write(*,*) '=============================================='
    write(*,*)
    write(*,'(A,F10.3,A)') ' Base Temperature (x=0):      ', T0, ' °C'
    write(*,'(A,F10.3,A)') ' Tip Temperature (x=L):       ', T_x, ' °C'
    write(*,'(A,F10.3,A)') ' Temperature Drop:            ', T0 - T_x, ' °C'
    write(*,*)
    
    ! Design recommendations
    write(*,*) '=============================================='
    write(*,*) '       DESIGN RECOMMENDATIONS'
    write(*,*) '=============================================='
    write(*,*)
    
    if (mL < 1.5d0) then
        write(*,*) '• Fin is relatively short (mL < 1.5)'
        write(*,*) '  Consider increasing fin length for better performance'
    else if (mL > 5.0d0) then
        write(*,*) '• Fin is very long (mL > 5.0)'
        write(*,*) '  Additional length provides minimal benefit'
        write(*,*) '  Consider reducing length or using multiple shorter fins'
    else
        write(*,*) '• Fin length is well optimized (1.5 < mL < 5.0)'
    end if
    write(*,*)
    
    if (eta < 0.6d0) then
        write(*,*) '• Low fin efficiency (η < 60%)'
        write(*,*) '  Recommendations:'
        write(*,*) '  - Reduce fin length'
        write(*,*) '  - Increase fin thickness'
        write(*,*) '  - Use material with higher thermal conductivity'
    else if (eta > 0.9d0) then
        write(*,*) '• Excellent fin efficiency (η > 90%)'
        write(*,*) '  Fin is performing optimally'
    else
        write(*,*) '• Good fin efficiency (60% < η < 90%)'
    end if
    write(*,*)
    
    ! Calculate optimal fin length
    write(*,*) '=============================================='
    write(*,*) '          OPTIMIZATION ANALYSIS'
    write(*,*) '=============================================='
    write(*,*)
    
    ! For maximum heat transfer per unit volume, mL ≈ 1.419
    write(*,'(A,F10.6,A)') ' Current mL value:            ', mL, ''
    write(*,'(A,F10.6,A)') ' Optimal mL (max Q/volume):   ', 1.419d0, ''
    
    if (abs(mL - 1.419d0) < 0.3d0) then
        write(*,*) ' Status: Near optimal design'
    else if (mL < 1.419d0) then
        write(*,'(A,F10.6,A)') ' Suggested length increase:   ', &
            (1.419d0 / m - L) * 1000.0d0, ' mm'
    else
        write(*,'(A,F10.6,A)') ' Suggested length decrease:   ', &
            (L - 1.419d0 / m) * 1000.0d0, ' mm'
    end if
    write(*,*)
    
    ! Material comparison
    write(*,*) '=============================================='
    write(*,*) '        MATERIAL COMPARISON'
    write(*,*) '=============================================='
    write(*,*)
    write(*,*) 'Heat transfer with different materials:'
    write(*,*) '(keeping all other parameters constant)'
    write(*,*)
    
    call material_comparison(h, P, A, theta_base, L, 'Aluminum', 200.0d0)
    call material_comparison(h, P, A, theta_base, L, 'Copper', 385.0d0)
    call material_comparison(h, P, A, theta_base, L, 'Steel', 50.0d0)
    call material_comparison(h, P, A, theta_base, L, 'Brass', 110.0d0)
    call material_comparison(h, P, A, theta_base, L, 'Cast Iron', 52.0d0)
    write(*,*)
    
    ! Energy savings calculation
    write(*,*) '=============================================='
    write(*,*) '         ENERGY ANALYSIS'
    write(*,*) '=============================================='
    write(*,*)
    write(*,'(A,F12.4,A)') ' Heat dissipation per fin:    ', Q, ' W'
    write(*,'(A,F12.4,A)') ' Heat per hour:               ', Q * 3600.0d0 / 1000.0d0, ' kJ/h'
    write(*,'(A,F12.4,A)') ' Heat per day:                ', Q * 86400.0d0 / 1000.0d0, ' kJ/day'
    write(*,*)
    
    ! Comparison without fin
    write(*,'(A,F12.4,A)') ' Heat transfer without fin:   ', h * A * theta_base, ' W'
    write(*,'(A,F12.4,A)') ' Improvement factor:          ', Q / (h * A * theta_base), ' times'
    write(*,*)
    
    write(*,*) '=============================================='
    write(*,*) '          CALCULATION COMPLETE'
    write(*,*) '=============================================='
    write(*,*)
    
contains

    subroutine material_comparison(h, P, A, theta_base, L, material_name, k_mat)
        implicit none
        real(8), intent(in) :: h, P, A, theta_base, L, k_mat
        character(len=*), intent(in) :: material_name
        real(8) :: m_mat, mL_mat, Q_mat, eta_mat
        
        m_mat = sqrt((h * P) / (k_mat * A))
        mL_mat = m_mat * L
        
        if (mL_mat > 3.0d0) then
            Q_mat = sqrt(h * P * k_mat * A) * theta_base
        else
            Q_mat = sqrt(h * P * k_mat * A) * theta_base * tanh(mL_mat)
        end if
        
        eta_mat = tanh(mL_mat) / mL_mat
        
        write(*,'(A,A12,A,F10.2,A,F8.2,A)') '  ', material_name, ':  Q = ', &
            Q_mat, ' W  (η = ', eta_mat * 100.0d0, '%)'
    end subroutine material_comparison

end program Rectangulaire_Fin_Heat_Conduction