program double_pipe_design
    implicit none

    double precision :: T_hi, T_ho, T_ci, T_co
    double precision :: m_h, m_c, Cp_h, Cp_c
    double precision :: D_i_in, D_i_out, D_o_in
    double precision :: k_w
    double precision :: h_i, h_o
    double precision :: R_fi, R_fo
    
    double precision :: Q_h, Q_c, Q_avg
    double precision :: dT1, dT2, LMTD
    double precision :: R_clean_sum, R_fouled_sum
    double precision :: U_clean, U_fouled
    double precision :: A_clean, A_fouled
    double precision :: L_clean, L_fouled
    double precision :: overdesign
    double precision :: discrepancy
    
    integer :: config
    integer :: i
    character(len=40) :: config_name
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs from stdin
    read(*,*) config     ! 1=Parallel, 2=Counter
    read(*,*) T_hi       ! Hot fluid inlet temp [C]
    read(*,*) T_ho       ! Hot fluid outlet temp [C]
    read(*,*) T_ci       ! Cold fluid inlet temp [C]
    read(*,*) T_co       ! Cold fluid outlet 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(*,*) D_i_in     ! Inner pipe inner diameter [m]
    read(*,*) D_i_out    ! Inner pipe outer diameter [m]
    read(*,*) D_o_in     ! Outer pipe inner diameter [m]
    read(*,*) k_w        ! Wall thermal conductivity [W/m.K]
    read(*,*) h_i        ! Inner pipe convection coefficient [W/m2.K]
    read(*,*) h_o        ! Outer annulus convection coefficient [W/m2.K]
    read(*,*) R_fi       ! Inner fouling factor [m2.K/W]
    read(*,*) R_fo       ! Outer fouling factor [m2.K/W]

    ! Input validation
    if (config < 1 .or. config > 2) 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 (T_hi <= T_ho) then
        write(*,'(A)') 'ERROR: Hot inlet temperature must be greater than hot outlet temperature.'
        stop
    end if
    if (T_co <= T_ci) then
        write(*,'(A)') 'ERROR: Cold outlet 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 (D_i_in <= 0.0d0 .or. D_i_out <= 0.0d0 .or. D_o_in <= 0.0d0) then
        write(*,'(A)') 'ERROR: Diameters must be positive.'
        stop
    end if
    if (D_i_in >= D_i_out) then
        write(*,'(A)') 'ERROR: Inner pipe inner diameter must be less than its outer diameter.'
        stop
    end if
    if (D_i_out >= D_o_in) then
        write(*,'(A)') 'ERROR: Inner pipe outer diameter must be less than outer pipe inner diameter.'
        stop
    end if
    if (k_w <= 0.0d0 .or. h_i <= 0.0d0 .or. h_o <= 0.0d0) then
        write(*,'(A)') 'ERROR: Conductivity and heat transfer coefficients must be positive.'
        stop
    end if
    if (R_fi < 0.0d0 .or. R_fo < 0.0d0) then
        write(*,'(A)') 'ERROR: Fouling factors cannot be negative.'
        stop
    end if

    ! Configuration name
    if (config == 1) then
        config_name = 'Parallel Flow (Concentric)'
    else
        config_name = 'Counter Flow (Concentric)'
    end if

    ! Heat duties
    Q_h = m_h * Cp_h * (T_hi - T_ho)
    Q_c = m_c * Cp_c * (T_co - T_ci)
    Q_avg = Q_h ! Use hot side duty as primary design load
    discrepancy = abs(Q_h - Q_c) / Q_h * 100.0d0

    ! 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
        write(*,'(A)') 'ERROR: Temperature cross occurred. Sizing impossible for parallel flow.'
        stop
    else if (abs(dT1 - dT2) < 1.0d-6) then
        LMTD = dT1
    else
        LMTD = (dT1 - dT2) / log(dT1 / dT2)
    end if

    ! Thermal Resistances and Overall Heat Transfer Coefficients (based on outer surface of inner tube, Ao)
    ! R_clean = 1/ho + (D_i_out * log(D_i_out/D_i_in)) / (2*kw) + (D_i_out/D_i_in) * (1/hi)
    R_clean_sum = (1.0d0 / h_o) + &
                  (D_i_out * log(D_i_out / D_i_in)) / (2.0d0 * k_w) + &
                  (D_i_out / D_i_in) * (1.0d0 / h_i)
                  
    ! R_fouled = R_clean + R_fo + R_fi * (D_i_out/D_i_in)
    R_fouled_sum = R_clean_sum + R_fo + R_fi * (D_i_out / D_i_in)

    U_clean = 1.0d0 / R_clean_sum
    U_fouled = 1.0d0 / R_fouled_sum

    ! Required surface areas (based on Ao)
    A_clean = Q_avg / (U_clean * LMTD)
    A_fouled = Q_avg / (U_fouled * LMTD)

    ! Required tube lengths
    L_clean = A_clean / (PI * D_i_out)
    L_fouled = A_fouled / (PI * D_i_out)

    ! Overdesign percentage
    overdesign = (A_fouled - A_clean) / A_clean * 100.0d0

    ! ============================================
    ! OUTPUT
    ! ============================================
    write(*,'(A)') '============================================================='
    write(*,'(A)') '   DOUBLE-PIPE HEAT EXCHANGER SIZING & RATING'
    write(*,'(A)') '============================================================='
    write(*,*)
    write(*,'(A)') '--- CONFIGURATION & FLOW TYPE ------------------------------'
    write(*,'(A,A)')        '  Flow Arrangement        = ', trim(config_name)
    write(*,*)
    write(*,'(A)') '--- GEOMETRY SPECIFICATIONS ---------------------------------'
    write(*,'(A,F12.2,A)')  '  Inner Pipe Inner Dia   = ', D_i_in*1000.0d0, ' mm'
    write(*,'(A,F12.2,A)')  '  Inner Pipe Outer Dia   = ', D_i_out*1000.0d0, ' mm'
    write(*,'(A,F12.2,A)')  '  Outer Shell Inner Dia  = ', D_o_in*1000.0d0, ' mm'
    write(*,'(A,F12.2,A)')  '  Wall Thermal Cond (kw)  = ', k_w, ' W/m.K'
    write(*,*)
    write(*,'(A)') '--- HEAT TRANSFER COEFFICIENTS & FOULING ---------------------'
    write(*,'(A,F12.2,A)')  '  Tube Convection (hi)    = ', h_i, ' W/m2.K'
    write(*,'(A,F12.2,A)')  '  Annulus Convection (ho) = ', h_o, ' W/m2.K'
    write(*,'(A,F12.6,A)')  '  Tube Fouling (Rfi)      = ', R_fi, ' m2.K/W'
    write(*,'(A,F12.6,A)')  '  Annulus Fouling (Rfo)   = ', R_fo, ' m2.K/W'
    write(*,*)
    write(*,'(A)') '--- THERMAL ANALYSIS & DUTY ---------------------------------'
    write(*,'(A,F12.2,A)')  '  Hot Fluid Inlet Temp    = ', T_hi, ' C'
    write(*,'(A,F12.2,A)')  '  Hot Fluid Outlet Temp   = ', T_ho, ' C'
    write(*,'(A,F12.2,A)')  '  Cold Fluid Inlet Temp   = ', T_ci, ' C'
    write(*,'(A,F12.2,A)')  '  Cold Fluid Outlet Temp  = ', T_co, ' C'
    write(*,'(A,F12.2,A)')  '  Calculated Heat Duty Q  = ', Q_avg, ' W'
    write(*,'(A,F12.2,A)')  '  Cold Side Absorbed Qc   = ', Q_c, ' W'
    write(*,'(A,F12.3,A)')  '  Energy Mismatch         = ', discrepancy, ' %'
    write(*,'(A,F12.2,A)')  '  LMTD                    = ', LMTD, ' deg-C'
    write(*,*)
    write(*,'(A)') '--- OVERALL HEAT TRANSFER COEFFICIENTS -----------------------'
    write(*,'(A,F12.2,A)')  '  U_clean (based on Ao)   = ', U_clean, ' W/m2.K'
    write(*,'(A,F12.2,A)')  '  U_fouled (based on Ao)  = ', U_fouled, ' W/m2.K'
    write(*,*)
    write(*,'(A)') '--- SIZING RESULTS ------------------------------------------'
    write(*,'(A,F12.4,A)')  '  Required Clean Area     = ', A_clean, ' m2'
    write(*,'(A,F12.4,A)')  '  Required Fouled Area    = ', A_fouled, ' m2'
    write(*,'(A,F12.2,A)')  '  Required Clean Length   = ', L_clean, ' m'
    write(*,'(A,F12.2,A)')  '  Required Fouled Length  = ', L_fouled, ' m'
    write(*,'(A,F12.2,A)')  '  Overdesign Factor       = ', overdesign, ' %'
    write(*,*)

    ! Temperature profiles for chart plotting
    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
            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 flow
            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 = m_dot * Cp * delta_T'
    write(*,'(A)') '  LMTD = (dT1 - dT2) / ln(dT1/dT2)'
    write(*,'(A)') '  1/U_clean = 1/ho + D_out*ln(D_out/D_in)/(2*kw) + (D_out/D_in)*(1/hi)'
    write(*,'(A)') '  1/U_fouled = 1/U_clean + Rfo + Rfi*(D_out/D_in)'
    write(*,'(A)') '  A = Q / (U * LMTD)'
    write(*,'(A)') '  L = A / (pi * D_out)'

end program double_pipe_design
