program heat_exchanger_ntu
    implicit none

    double precision :: T_hi, T_ho, T_ci, T_co
    double precision :: m_h, m_c, Cp_h, Cp_c
    double precision :: U, A
    double precision :: Q_max, Q
    double precision :: dT1, dT2, LMTD
    double precision :: C_h, C_c, C_min, C_max, C_r
    double precision :: eff
    double precision :: NTU
    integer :: config
    integer :: i
    character(len=40) :: config_name
    double precision :: factor, term

    ! Read inputs from stdin
    read(*,*) config     ! 1=Parallel, 2=Counter, 3=Shell-and-Tube (1-2), 4=Cross flow (unmixed)
    read(*,*) T_hi       ! Hot fluid inlet temp [C]
    read(*,*) T_ci       ! Cold fluid inlet temp [C]
    read(*,*) m_h        ! Hot fluid mass flow rate [kg/s]
    read(*,*) m_c        ! Cold fluid mass flow rate [kg/s]
    read(*,*) Cp_h       ! Hot fluid specific heat [J/kg.K]
    read(*,*) Cp_c       ! Cold fluid specific heat [J/kg.K]
    read(*,*) U          ! Overall heat transfer coefficient [W/m2.K]
    read(*,*) A          ! Heat transfer area [m2]

    ! Input validation
    if (config < 1 .or. config > 4) then
        write(*,'(A)') 'ERROR: Invalid configuration selector.'
        stop
    end if
    if (T_hi <= T_ci) then
        write(*,'(A)') 'ERROR: Hot inlet temperature must be greater than cold inlet temperature.'
        stop
    end if
    if (m_h <= 0.0d0 .or. m_c <= 0.0d0 .or. Cp_h <= 0.0d0 .or. Cp_c <= 0.0d0) then
        write(*,'(A)') 'ERROR: Mass flow rates and specific heats must be positive.'
        stop
    end if
    if (U <= 0.0d0 .or. A <= 0.0d0) then
        write(*,'(A)') 'ERROR: Overall heat transfer coefficient and area must be positive.'
        stop
    end if

    ! Configuration name
    select case (config)
    case (1)
        config_name = 'Parallel Flow'
    case (2)
        config_name = 'Counter Flow'
    case (3)
        config_name = 'Shell-and-Tube (1-2 Pass)'
    case (4)
        config_name = 'Cross Flow (Unmixed)'
    end select

    ! Heat capacity rates
    C_h = m_h * Cp_h
    C_c = m_c * Cp_c

    if (C_h < C_c) then
        C_min = C_h
        C_max = C_c
    else
        C_min = C_c
        C_max = C_h
    end if
    C_r = C_min / C_max

    ! NTU
    NTU = U * A / C_min

    ! ============================================
    ! EFFECTIVENESS CALCULATION
    ! ============================================
    select case (config)
    case (1)
        ! Parallel Flow
        eff = (1.0d0 - exp(-NTU * (1.0d0 + C_r))) / (1.0d0 + C_r)

    case (2)
        ! Counter Flow
        if (C_r < 1.0d-6) then
            eff = 1.0d0 - exp(-NTU)
        else if (abs(C_r - 1.0d0) < 1.0d-6) then
            eff = NTU / (1.0d0 + NTU)
        else
            eff = (1.0d0 - exp(-NTU * (1.0d0 - C_r))) / (1.0d0 - C_r * exp(-NTU * (1.0d0 - C_r)))
        end if

    case (3)
        ! Shell-and-Tube (1 shell pass, 2 tube passes)
        factor = sqrt(1.0d0 + C_r**2)
        term = exp(-NTU * factor)
        if (abs(1.0d0 - term) < 1.0d-12) then
            eff = 0.0d0
        else
            eff = 2.0d0 / (1.0d0 + C_r + factor * (1.0d0 + term) / (1.0d0 - term))
        end if

    case (4)
        ! Cross Flow (unmixed)
        if (C_r < 1.0d-6) then
            eff = 1.0d0 - exp(-NTU)
        else
            eff = 1.0d0 - exp( (NTU**0.22d0) / C_r * (exp(-C_r * (NTU**0.78d0)) - 1.0d0) )
        end if
    end select

    ! Limit effectiveness to physical bounds [0, 1]
    if (eff < 0.0d0) eff = 0.0d0
    if (eff > 1.0d0) eff = 1.0d0

    ! Heat transfer rates
    Q_max = C_min * (T_hi - T_ci)
    Q = eff * Q_max

    ! Outlet temperatures
    T_ho = T_hi - Q / C_h
    T_co = T_ci + Q / C_c

    ! LMTD Calculation
    if (config == 1) then
        dT1 = T_hi - T_ci
        dT2 = T_ho - T_co
    else
        dT1 = T_hi - T_co
        dT2 = T_ho - T_ci
    end if

    if (dT1 <= 0.0d0 .or. dT2 <= 0.0d0) then
        LMTD = 0.0d0
    else if (abs(dT1 - dT2) < 1.0d-6) then
        LMTD = dT1
    else
        LMTD = (dT1 - dT2) / log(dT1 / dT2)
    end if

    ! ============================================
    ! OUTPUT
    ! ============================================
    write(*,'(A)') '============================================================='
    write(*,'(A)') '   HEAT EXCHANGER ANALYSIS - E-NTU METHOD'
    write(*,'(A)') '============================================================='
    write(*,*)
    write(*,'(A)') '--- CONFIGURATION -------------------------------------------'
    write(*,'(A,A)')        '  Flow Arrangement        = ', trim(config_name)
    write(*,*)
    write(*,'(A)') '--- INPUT PARAMETERS ----------------------------------------'
    write(*,'(A)') '  Hot Fluid:'
    write(*,'(A,F12.2,A)')  '    Inlet Temperature     = ', T_hi, ' C'
    write(*,'(A,F12.4,A)')  '    Mass Flow Rate        = ', m_h, ' kg/s'
    write(*,'(A,F12.2,A)')  '    Specific Heat         = ', Cp_h, ' J/kg.K'
    write(*,*)
    write(*,'(A)') '  Cold Fluid:'
    write(*,'(A,F12.2,A)')  '    Inlet Temperature     = ', T_ci, ' C'
    write(*,'(A,F12.4,A)')  '    Mass Flow Rate        = ', m_c, ' kg/s'
    write(*,'(A,F12.2,A)')  '    Specific Heat         = ', Cp_c, ' J/kg.K'
    write(*,*)
    write(*,'(A,F12.2,A)')  '  Overall U               = ', U, ' W/m2.K'
    write(*,'(A,F12.4,A)')  '  Surface Area A          = ', A, ' m2'
    write(*,*)
    write(*,'(A)') '--- HEAT CAPACITY RATES -------------------------------------'
    write(*,'(A,F12.4,A)')  '  C_hot  = m_h x Cp_h    = ', C_h, ' W/K'
    write(*,'(A,F12.4,A)')  '  C_cold = m_c x Cp_c    = ', C_c, ' W/K'
    write(*,'(A,F12.4,A)')  '  C_min                   = ', C_min, ' W/K'
    write(*,'(A,F12.4,A)')  '  C_max                   = ', C_max, ' W/K'
    write(*,'(A,F12.4)')    '  C_r = C_min/C_max       = ', C_r
    write(*,*)
    write(*,'(A)') '--- ANALYSIS RESULTS ----------------------------------------'
    write(*,'(A,F12.4)')    '  NTU = UA/C_min          = ', NTU
    write(*,'(A,F12.4,A)')  '  Effectiveness (e)       = ', eff*100.0d0, ' %'
    write(*,'(A,F12.4,A)')  '  Max Heat Transfer Qmax  = ', Q_max, ' W'
    write(*,'(A,F12.4,A)')  '  Actual Heat Transfer Q  = ', Q, ' W'
    write(*,'(A,F12.2,A)')  '  Hot Outlet Temperature  = ', T_ho, ' C'
    write(*,'(A,F12.2,A)')  '  Cold Outlet Temperature = ', T_co, ' C'
    write(*,'(A,F12.4,A)')  '  LMTD                    = ', LMTD, ' deg-C'
    write(*,*)

    ! Temperature profiles for chart
    write(*,'(A)') '--- TEMPERATURE PROFILE ALONG EXCHANGER --------------------'
    write(*,'(A)') '  Position %    T_hot [C]    T_cold [C]'
    write(*,'(A)') '  --------------------------------------------------'

    do i = 0, 20
        if (config == 1) then
            ! Parallel flow: both fluids travel same direction
            write(*,'(F10.1,4X,F10.2,4X,F10.2)') dble(i)/20.0d0*100, &
                T_hi - (T_hi - T_ho) * dble(i)/20.0d0, &
                T_ci + (T_co - T_ci) * dble(i)/20.0d0
        else
            ! Counter/Shell/Cross flow: cold fluid travels opposite direction
            write(*,'(F10.1,4X,F10.2,4X,F10.2)') dble(i)/20.0d0*100, &
                T_hi - (T_hi - T_ho) * dble(i)/20.0d0, &
                T_co - (T_co - T_ci) * dble(i)/20.0d0
        end if
    end do

    write(*,*)
    write(*,'(A)') '--- FORMULAS USED -------------------------------------------'
    write(*,'(A)') '  Q = e x Q_max'
    write(*,'(A)') '  Q_max = C_min x (T_hi - T_ci)'
    write(*,'(A)') '  NTU = U x A / C_min'
    write(*,'(A)') '  C_r = C_min / C_max'

end program heat_exchanger_ntu
