program multicomponent_distillation_fug
    implicit none
    integer :: iostat_val, iter
    double precision :: alpha_LK, alpha_HK, zF_LK, zF_HK, xD_LK, xB_LK
    double precision :: q_feed, reflux_factor, F_kmol_h
    double precision :: N_min, theta, f_val, df_val, R_min, R_op, X_gil, Y_gil, N_theor
    double precision :: D_kmol_h, B_kmol_h, N_rect, N_strip, feed_stage
    double precision, parameter :: TOL = 1.0d-8

    ! Read inputs
    read(*,*,iostat=iostat_val) alpha_LK       ! Relative volatility of Light Key (LK) wrt HK (e.g. 2.45)
    read(*,*,iostat=iostat_val) zF_LK          ! Feed mole fraction of LK (e.g. 0.40)
    read(*,*,iostat=iostat_val) zF_HK          ! Feed mole fraction of Heavy Key (HK) (e.g. 0.60)
    read(*,*,iostat=iostat_val) xD_LK          ! Distillate target purity of LK (e.g. 0.98)
    read(*,*,iostat=iostat_val) xB_LK          ! Bottoms residual fraction of LK (e.g. 0.02)
    read(*,*,iostat=iostat_val) q_feed         ! Thermal state of feed (1.0=Sat liquid, 0.0=Sat vapor, e.g. 1.0)
    read(*,*,iostat=iostat_val) reflux_factor  ! Reflux multiplier R / Rmin (e.g. 1.30)
    read(*,*,iostat=iostat_val) F_kmol_h       ! Total Feed Flow Rate [kmol/h] (e.g. 100.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for FUG distillation calculation.'
        stop
    end if

    if (alpha_LK <= 1.0d0 .or. xD_LK <= xB_LK .or. reflux_factor <= 1.0d0) then
        write(*,*) 'ERROR: alpha_LK must be > 1.0, xD_LK > xB_LK, and reflux_factor > 1.0'
        stop
    end if

    alpha_HK = 1.0d0

    ! 1. Fenske Equation: Minimum stages N_min
    N_min = log((xD_LK / (1.0d0 - xD_LK)) * ((1.0d0 - xB_LK) / xB_LK)) / log(alpha_LK)

    ! Component overall mass balance (Key components)
    D_kmol_h = F_kmol_h * (zF_LK - xB_LK) / (xD_LK - xB_LK)
    B_kmol_h = F_kmol_h - D_kmol_h

    ! 2. Underwood Equation: Solve root theta by Newton-Raphson
    theta = 1.0d0 + (alpha_LK - 1.0d0) * 0.5d0 ! initial guess between 1.0 and alpha_LK
    do iter = 1, 100
        f_val = (alpha_LK * zF_LK / (alpha_LK - theta)) + (alpha_HK * zF_HK / (alpha_HK - theta)) - (1.0d0 - q_feed)
        df_val = (alpha_LK * zF_LK / ((alpha_LK - theta)**2)) + (alpha_HK * zF_HK / ((alpha_HK - theta)**2))
        theta = theta - f_val / df_val
        if (theta <= 1.0d0) theta = 1.0001d0
        if (theta >= alpha_LK) theta = alpha_LK - 0.0001d0
        if (abs(f_val) < TOL) exit
    end do

    ! Calculate R_min
    R_min = (alpha_LK * xD_LK / (alpha_LK - theta)) + (alpha_HK * (1.0d0 - xD_LK) / (alpha_HK - theta)) - 1.0d0
    if (R_min < 0.1d0) R_min = 0.5d0

    R_op = reflux_factor * R_min

    ! 3. Gilliland (Molokanov) Correlation for Actual Stages N_theor
    X_gil = (R_op - R_min) / (R_op + 1.0d0)
    Y_gil = 1.0d0 - exp(((1.0d0 + 54.4d0 * X_gil) / (11.0d0 + 117.2d0 * X_gil)) * ((X_gil - 1.0d0) / sqrt(X_gil)))
    N_theor = (N_min + Y_gil) / (1.0d0 - Y_gil)

    ! 4. Kirkbride Equation for Feed Stage Location
    N_rect = N_theor * (((zF_HK / zF_LK) * (xB_LK / (1.0d0 - xD_LK))**2 * (B_kmol_h / D_kmol_h))**0.206d0) / &
             (1.0d0 + (((zF_HK / zF_LK) * (xB_LK / (1.0d0 - xD_LK))**2 * (B_kmol_h / D_kmol_h))**0.206d0))
    N_strip = N_theor - N_rect
    feed_stage = N_rect + 1.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — MULTICOMPONENT DISTILLATION (FUG)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2)')    'Relative Volatility alpha(LK/HK) = ', alpha_LK
    write(*,'(A,F10.3,A,F10.3)') 'Feed zF(LK) / Distillate xD(LK)  = ', zF_LK, ' / ', xD_LK
    write(*,'(A,F10.3,A,F10.2)') 'Feed Thermal State (q) / Reflux x= ', q_feed, ' / ', reflux_factor
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Minimum Theoretical Stages (Nmin)= ', N_min, ' stages (Total Reflux)'
    write(*,'(A,F10.3)')    'Underwood Root (theta)           = ', theta
    write(*,'(A,F10.3)')    'Minimum Reflux Ratio (Rmin)      = ', R_min
    write(*,'(A,F10.3)')    'Operating Reflux Ratio (R)       = ', R_op
    write(*,'(A,F10.2,A)')  'Actual Theoretical Stages (N)    = ', N_theor, ' stages'
    write(*,'(A,F10.1,A)')  'Optimal Feed Tray Location       = Stage ', feed_stage, ' from top'
    write(*,'(A,F10.2,A,F10.2,A)') 'Distillate / Bottoms Flow Rates  = ', D_kmol_h, ' kmol/h / ', B_kmol_h, ' kmol/h'
    write(*,'(A)') '============================================================'

end program multicomponent_distillation_fug