program heat_mass_analogy
    implicit none

    ! Inputs
    integer :: geom_type ! 1 = External flow over flat plate, 2 = Internal flow inside tube
    double precision :: U_vel ! Free-stream or mean velocity, m/s
    double precision :: L_c ! Characteristic length (L for plate, D for tube), m
    double precision :: rho ! Density, kg/m3
    double precision :: mu ! Dynamic viscosity, Pa.s
    double precision :: k_f ! Thermal conductivity, W/m.K
    double precision :: Cp ! Specific heat, J/kg.K
    double precision :: D_AB ! Binary mass diffusivity, m2/s

    ! Physical properties & Dimensionless Numbers
    double precision :: nu, alpha
    double precision :: Re, Pr, Sc, Le
    double precision :: Cf, Nu_avg, Sh_avg
    double precision :: h_avg, hm_avg
    double precision :: ratio_actual, ratio_theory

    ! Profile variables
    integer :: i, n_points
    double precision :: x, dx, Re_x, Nu_x, Sh_x, h_x, hm_x, f_factor
    character(len=50) :: geom_name
    integer :: iostat_val

    ! Read inputs
    read(*,*,iostat=iostat_val) geom_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid geometry selection.'
        stop
    end if

    read(*,*,iostat=iostat_val) U_vel
    read(*,*,iostat=iostat_val) L_c
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu
    read(*,*,iostat=iostat_val) k_f
    read(*,*,iostat=iostat_val) Cp
    read(*,*,iostat=iostat_val) D_AB

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all analogy parameters.'
        stop
    end if

    ! Basic validations
    if (U_vel <= 0.0d0) then
        write(*,*) 'ERROR: Velocity must be positive.'
        stop
    end if
    if (L_c <= 0.0d0) then
        write(*,*) 'ERROR: Characteristic length must be positive.'
        stop
    end if
    if (rho <= 0.0d0 .or. mu <= 0.0d0 .or. k_f <= 0.0d0 .or. Cp <= 0.0d0 .or. D_AB <= 0.0d0) then
        write(*,*) 'ERROR: Physical properties must be positive.'
        stop
    end if

    ! Compute physical properties
    nu = mu / rho
    alpha = k_f / (rho * Cp)

    ! Compute dimensionless numbers
    Re = U_vel * L_c / nu
    Pr = nu / alpha
    Sc = nu / D_AB
    Le = alpha / D_AB ! Le = Sc / Pr

    ! Setup geometry cases
    if (geom_type == 1) then
        geom_name = "External Flow over Flat Plate"
        
        ! Transition Reynolds number is 5.0e5
        if (Re <= 5.0d5) then
            ! Fully laminar flow
            Cf = 1.328d0 / sqrt(Re)
            Nu_avg = 0.664d0 * sqrt(Re) * Pr**(1.0d0/3.0d0)
            Sh_avg = 0.664d0 * sqrt(Re) * Sc**(1.0d0/3.0d0)
        else
            ! Mixed laminar/turbulent flow
            Cf = 0.074d0 / (Re**0.2d0) - 1740.0d0 / Re
            Nu_avg = (0.037d0 * Re**0.8d0 - 871.0d0) * Pr**(1.0d0/3.0d0)
            Sh_avg = (0.037d0 * Re**0.8d0 - 871.0d0) * Sc**(1.0d0/3.0d0)
        end if

        h_avg = Nu_avg * k_f / L_c
        hm_avg = Sh_avg * D_AB / L_c
        ratio_theory = rho * Cp * Le**(2.0d0/3.0d0)

    else if (geom_type == 2) then
        geom_name = "Internal Flow in Circular Tube"
        
        if (Re <= 2300.0d0) then
            ! Laminar fully developed pipe flow (constant wall temp assumption)
            Cf = 16.0d0 / Re ! f = 64/Re, Cf = f/4
            Nu_avg = 3.66d0
            Sh_avg = 3.66d0
            ratio_theory = rho * Cp * Le**(2.0d0/3.0d0)
        else
            ! Turbulent flow - Petukhov correlation for friction factor
            ! f = (0.790 * ln(Re) - 1.64)**(-2)
            f_factor = (0.790d0 * log(Re) - 1.64d0)**(-2)
            Cf = f_factor / 4.0d0
            
            ! Dittus-Boelter correlation (standard heating: exponent = 0.4)
            Nu_avg = 0.023d0 * Re**0.8d0 * Pr**0.4d0
            Sh_avg = 0.023d0 * Re**0.8d0 * Sc**0.4d0
            ratio_theory = rho * Cp * Le**(0.6d0) ! Modified exponent due to Dittus-Boelter exponent 0.4
        end if

        h_avg = Nu_avg * k_f / L_c
        hm_avg = Sh_avg * D_AB / L_c

    else
        write(*,*) 'ERROR: Invalid geometry selection code.'
        stop
    end if

    ratio_actual = h_avg / hm_avg

    ! Output results
    write(*,*) '============================================================'
    write(*,*) '      HEAT & MASS TRANSFER ANALOGY SOLVER'
    write(*,*) '============================================================'
    write(*,*)
    write(*,'(A,A)')        '  Configuration Type       = ', trim(geom_name)
    write(*,'(A,F12.4,A)')  '  Velocity (U)             = ', U_vel, ' m/s'
    write(*,'(A,F12.4,A)')  '  Characteristic Length    = ', L_c, ' m'
    write(*,*)
    write(*,'(A)') '--- FLUID & DIFFUSION PROPERTIES ---------------------------'
    write(*,'(A,F12.4,A)')  '  Density (rho)            = ', rho, ' kg/m3'
    write(*,'(A,ES12.4,A)') '  Viscosity (mu)           = ', mu, ' Pa.s'
    write(*,'(A,F12.4,A)')  '  Thermal Cond (kf)        = ', k_f, ' W/m.K'
    write(*,'(A,F12.2,A)')  '  Specific Heat (Cp)       = ', Cp, ' J/kg.K'
    write(*,'(A,ES12.4,A)') '  Mass Diffusivity D_AB    = ', D_AB, ' m2/s'
    write(*,*)
    write(*,'(A)') '--- DIMENSIONLESS NUMBERS ----------------------------------'
    write(*,'(A,ES12.4)')   '  Reynolds Number (Re)     = ', Re
    write(*,'(A,F12.4)')    '  Prandtl Number (Pr)      = ', Pr
    write(*,'(A,F12.4)')    '  Schmidt Number (Sc)      = ', Sc
    write(*,'(A,F12.4)')    '  Lewis Number (Le)        = ', Le
    write(*,*)
    write(*,'(A)') '--- ANALOGY RESULTS (AVERAGE VALUES) -----------------------'
    write(*,'(A,F12.6)')    '  Skin Friction Coeff (Cf) = ', Cf
    write(*,'(A,F12.4)')    '  Average Nusselt Number   = ', Nu_avg
    write(*,'(A,F12.4)')    '  Average Sherwood Number  = ', Sh_avg
    write(*,'(A,F12.4,A)')  '  Average Heat HTC (h)     = ', h_avg, ' W/m2.K'
    write(*,'(A,ES12.4,A)') '  Average Mass MTC (hm)    = ', hm_avg, ' m/s'
    write(*,*)
    write(*,'(A)') '--- CHILTON-COLBURN RATIO VERIFICATION ---------------------'
    write(*,'(A,ES12.4,A)') '  Actual Ratio h/hm        = ', ratio_actual, ' J/m3.K'
    write(*,'(A,ES12.4,A)') '  Theoretical rho*Cp*Le^n  = ', ratio_theory, ' J/m3.K'
    write(*,'(A,F12.2,A)')  '  Analogy Deviation        = ', abs(ratio_actual - ratio_theory)/ratio_theory * 100.0d0, ' %'
    write(*,*)

    ! ============================================
    ! LOCAL PROFILE DATA (along flow path)
    ! ============================================
    write(*,'(A)') '--- LOCAL PROFILE ALONG SURFACE ----------------------------'
    write(*,'(A)') '  x [m]       Re_x          Nu_x         Sh_x         h_x [W/m2.K]  hm_x [m/s]'
    write(*,'(A)') '  ----------------------------------------------------------------------'

    n_points = 40
    dx = L_c / dble(n_points)

    do i = 1, n_points
        x = dble(i) * dx
        Re_x = U_vel * x / nu

        if (geom_type == 1) then
            ! Flat Plate boundary layer development
            if (Re_x <= 5.0d5) then
                ! Local laminar Nusselt & Sherwood
                Nu_x = 0.332d0 * sqrt(Re_x) * Pr**(1.0d0/3.0d0)
                Sh_x = 0.332d0 * sqrt(Re_x) * Sc**(1.0d0/3.0d0)
            else
                ! Local turbulent Nusselt & Sherwood
                Nu_x = 0.0296d0 * Re_x**0.8d0 * Pr**(1.0d0/3.0d0)
                Sh_x = 0.0296d0 * Re_x**0.8d0 * Sc**(1.0d0/3.0d0)
            end if
            h_x = Nu_x * k_f / x
            hm_x = Sh_x * D_AB / x
        else
            ! Internal Flow - tube entrance region thermal development
            ! Use simplified thermal entrance correlation for laminar/turbulent
            if (Re <= 2300.0d0) then
                ! Laminar entrance region
                Nu_x = 3.66d0 + (0.0668d0 * (L_c/x) * Re * Pr) / (1.0d0 + 0.04d0 * ((L_c/x) * Re * Pr)**(2.0d0/3.0d0))
                Sh_x = 3.66d0 + (0.0668d0 * (L_c/x) * Re * Sc) / (1.0d0 + 0.04d0 * ((L_c/x) * Re * Sc)**(2.0d0/3.0d0))
            else
                ! Turbulent entrance (approximate entrance region effect)
                Nu_x = Nu_avg * (1.0d0 + (L_c/x)**0.7d0) ! Simplistic entry-region factor
                Sh_x = Sh_avg * (1.0d0 + (L_c/x)**0.7d0)
            end if
            h_x = Nu_x * k_f / L_c
            hm_x = Sh_x * D_AB / L_c
        end if

        write(*,'(F8.4,2X,ES12.4,2X,F10.2,3X,F10.2,4X,F10.4,4X,ES12.4)') &
            x, Re_x, Nu_x, Sh_x, h_x, hm_x
    end do

    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    if (geom_type == 1) then
        write(*,'(A)') '  Flat Plate Laminar (Re <= 5e5): Nu_x = 0.332 Re_x^0.5 Pr^1/3, Sh_x = 0.332 Re_x^0.5 Sc^1/3.'
        write(*,'(A)') '  Flat Plate Turbulent (Re > 5e5): Nu_x = 0.0296 Re_x^0.8 Pr^1/3, Sh_x = 0.0296 Re_x^0.8 Sc^1/3.'
    else
        write(*,'(A)') '  Laminar Tube (Re <= 2300): Nu = 3.66 (Fully Developed Constant Temperature).'
        write(*,'(A)') '  Turbulent Tube (Re > 2300): Dittus-Boelter correlation Nu = 0.023 Re^0.8 Pr^n.'
    end if
    write(*,'(A)') '  Boundary layers assumed analogous under Chilton-Colburn definition.'

end program heat_mass_analogy
