﻿program internal_flow_convection
    implicit none
    
    double precision :: D, L, T_s, T_in, m_dot, rho, mu, nu, k_f, Pr, Cp
    double precision :: A_c, P_w, U_m, Re_D, f
    double precision :: x_fd_h, x_fd_t
    double precision :: Nu_DB, Nu_Gn, Nu_used, h_avg
    double precision :: T_out, T_lm, Q_total
    double precision :: dP, Pump_power
    character(len=30) :: flow_regime, correlation_used
    double precision :: pi
    integer :: i, n_points
    double precision :: x, dx, T_x, T_mx
    
    pi = 3.14159265358979d0
    
    ! Read inputs
    read(*,*) D        ! Tube inner diameter [m]
    read(*,*) L        ! Tube length [m]
    read(*,*) T_s      ! Surface temperature [°C] (constant wall temp)
    read(*,*) T_in     ! Inlet fluid temperature [°C]
    read(*,*) m_dot    ! Mass flow rate [kg/s]
    read(*,*) rho      ! Fluid density [kg/m³]
    read(*,*) mu       ! Dynamic viscosity [Pa·s]
    read(*,*) k_f      ! Thermal conductivity [W/m·K]
    read(*,*) Pr       ! Prandtl number
    read(*,*) Cp       ! Specific heat [J/kg·K]
    
    ! Computed properties
    nu = mu / rho
    A_c = pi * D**2 / 4.0d0
    P_w = pi * D
    U_m = m_dot / (rho * A_c)
    Re_D = rho * U_m * D / mu
    
    ! Hydrodynamic entry length
    if (Re_D < 2300.0d0) then
        flow_regime = 'Laminar'
        x_fd_h = 0.05d0 * Re_D * D
        x_fd_t = 0.05d0 * Re_D * Pr * D
    else if (Re_D < 4000.0d0) then
        flow_regime = 'Transitional'
        x_fd_h = 10.0d0 * D
        x_fd_t = 10.0d0 * D
    else
        flow_regime = 'Turbulent'
        x_fd_h = 10.0d0 * D
        x_fd_t = 10.0d0 * D
    end if
    
    ! ============================================
    ! FRICTION FACTOR
    ! ============================================
    if (Re_D < 2300.0d0) then
        ! Laminar: f = 64/Re
        f = 64.0d0 / Re_D
    else
        ! Turbulent: Petukhov correlation
        f = (0.790d0 * log(Re_D) - 1.64d0)**(-2.0d0)
    end if
    
    ! ============================================
    ! NUSSELT NUMBER
    ! ============================================
    if (Re_D < 2300.0d0) then
        ! Laminar, fully developed, constant wall temperature
        Nu_DB = 3.66d0
        Nu_Gn = 3.66d0
        Nu_used = 3.66d0
        correlation_used = 'Nu = 3.66 (const T_s)'
    else
        ! Dittus-Boelter: Nu = 0.023 Re^0.8 Pr^n (n=0.4 for heating)
        if (T_s > T_in) then
            Nu_DB = 0.023d0 * Re_D**0.8d0 * Pr**0.4d0
        else
            Nu_DB = 0.023d0 * Re_D**0.8d0 * Pr**0.3d0
        end if
        
        ! Gnielinski: Nu = (f/8)(Re-1000)Pr / [1 + 12.7(f/8)^0.5 (Pr^(2/3)-1)]
        Nu_Gn = ((f/8.0d0) * (Re_D - 1000.0d0) * Pr) / &
                (1.0d0 + 12.7d0 * sqrt(f/8.0d0) * (Pr**(2.0d0/3.0d0) - 1.0d0))
        
        ! Use Gnielinski (more accurate for transition)
        if (Re_D >= 3000.0d0) then
            Nu_used = Nu_Gn
            correlation_used = 'Gnielinski'
        else
            Nu_used = Nu_DB
            correlation_used = 'Dittus-Boelter'
        end if
    end if
    
    ! Average heat transfer coefficient
    h_avg = Nu_used * k_f / D
    
    ! ============================================
    ! OUTLET TEMPERATURE (constant wall temp)
    ! ============================================
    ! T_out = T_s - (T_s - T_in) * exp(-P_w * L * h_avg / (m_dot * Cp))
    T_out = T_s - (T_s - T_in) * exp(-pi * D * L * h_avg / (m_dot * Cp))
    
    ! Log-mean temperature difference
    if (abs(T_s - T_out) > 1.0d-6 .and. abs(T_s - T_in) > 1.0d-6) then
        T_lm = ((T_s - T_in) - (T_s - T_out)) / log((T_s - T_in) / (T_s - T_out))
    else
        T_lm = abs(T_s - T_in)
    end if
    
    ! Total heat transfer
    Q_total = m_dot * Cp * (T_out - T_in)
    
    ! Pressure drop
    dP = f * (L / D) * rho * U_m**2 / 2.0d0
    
    ! Pumping power
    Pump_power = dP * m_dot / rho
    
    ! ============================================
    ! OUTPUT
    ! ============================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   FORCED CONVECTION - INTERNAL FLOW IN CIRCULAR TUBE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUT PARAMETERS ---'
    write(*,'(A,F12.4,A)') '  Tube Diameter (D)       = ', D*1000, ' mm'
    write(*,'(A,F12.4,A)') '  Tube Length (L)         = ', L, ' m'
    write(*,'(A,F12.2,A)') '  Wall Temperature (T_s)  = ', T_s, ' °C'
    write(*,'(A,F12.2,A)') '  Inlet Temperature       = ', T_in, ' °C'
    write(*,'(A,ES12.4,A)') '  Mass Flow Rate          = ', m_dot, ' kg/s'
    write(*,*)
    write(*,'(A)') '--- FLUID PROPERTIES ---'
    write(*,'(A,F12.2,A)')  '  Density                 = ', rho, ' kg/m³'
    write(*,'(A,ES12.4,A)') '  Dynamic Viscosity       = ', mu, ' Pa·s'
    write(*,'(A,F12.4,A)')  '  Thermal Conductivity    = ', k_f, ' W/m·K'
    write(*,'(A,F12.4)')    '  Prandtl Number          = ', Pr
    write(*,'(A,F12.2,A)')  '  Specific Heat           = ', Cp, ' J/kg·K'
    write(*,*)
    write(*,'(A)') '--- FLOW ANALYSIS ---'
    write(*,'(A,F12.4,A)')  '  Mean Velocity           = ', U_m, ' m/s'
    write(*,'(A,ES12.4)')   '  Reynolds Number (Re_D)  = ', Re_D
    write(*,'(A,A)')        '  Flow Regime             = ', trim(flow_regime)
    write(*,'(A,ES12.4)')   '  Friction Factor (f)     = ', f
    write(*,'(A,F12.4,A)')  '  Entry Length (hydro)    = ', x_fd_h, ' m'
    write(*,'(A,F12.4,A)')  '  Entry Length (thermal)  = ', x_fd_t, ' m'
    write(*,*)
    write(*,'(A)') '--- HEAT TRANSFER RESULTS ---'
    write(*,'(A,A)')        '  Correlation Used        = ', trim(correlation_used)
    if (Re_D >= 2300.0d0) then
        write(*,'(A,F12.4)')    '  Nu (Dittus-Boelter)     = ', Nu_DB
        write(*,'(A,F12.4)')    '  Nu (Gnielinski)         = ', Nu_Gn
    end if
    write(*,'(A,F12.4)')    '  Nusselt Number Used     = ', Nu_used
    write(*,'(A,F12.4,A)')  '  Average h               = ', h_avg, ' W/m²·K'
    write(*,'(A,F12.2,A)')  '  Outlet Temperature      = ', T_out, ' °C'
    write(*,'(A,F12.4,A)')  '  LMTD                    = ', T_lm, ' °C'
    write(*,'(A,F12.4,A)')  '  Total Heat Transfer (Q) = ', Q_total, ' W'
    write(*,*)
    write(*,'(A)') '--- PRESSURE DROP ---'
    write(*,'(A,F12.4,A)')  '  Pressure Drop (ΔP)     = ', dP, ' Pa'
    write(*,'(A,F12.4,A)')  '  Pressure Drop           = ', dP/1000.0d0, ' kPa'
    write(*,'(A,F12.6,A)')  '  Pumping Power           = ', Pump_power, ' W'
    write(*,*)
    
    ! Temperature profile along tube
    write(*,'(A)') '--- TEMPERATURE PROFILE ALONG TUBE ---'
    write(*,'(A)') '  x [m]       T_bulk [°C]   T_wall [°C]   ΔT [°C]'
    write(*,'(A)') '  ---'
    
    n_points = 30
    dx = L / dble(n_points)
    
    do i = 0, n_points
        x = dble(i) * dx
        ! T_m(x) = T_s - (T_s - T_in) * exp(-P_w * x * h_avg / (m_dot * Cp))
        T_mx = T_s - (T_s - T_in) * exp(-pi * D * x * h_avg / (m_dot * Cp))
        write(*,'(F8.3,4X,F10.2,6X,F10.2,6X,F10.2)') x, T_mx, T_s, abs(T_s - T_mx)
    end do
    
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---'
    write(*,'(A)') '  Laminar:     Nu = 3.66 (const T_s, fully developed)'
    write(*,'(A)') '  Dittus-B:    Nu = 0.023 Re^0.8 Pr^n (n=0.4 heat, 0.3 cool)'
    write(*,'(A)') '  Gnielinski:  Nu = (f/8)(Re-1000)Pr/[1+12.7√(f/8)(Pr^(2/3)-1)]'
    write(*,'(A)') '  Valid for:   0.5 < Pr < 2000, 3000 < Re < 5×10⁶'
    
end program internal_flow_convection
