McCabe-Thiele Distillation Column

program distillation_mccabe_thiele
    implicit none
    integer :: i, n, maxStages, feedStage
    double precision :: alpha, xD, xB, zF, q, Rfactor, F, Rmin, R, Nmin
    double precision :: xq, yq, mR, bR, mS, bS, x, y, xeq, ynew, D, B
    double precision :: xprev, eps
    character(len=20) :: section
    integer :: iostat_val

    read(*,*,iostat=iostat_val) alpha
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid relative volatility input.'
        stop
    end if
    read(*,*,iostat=iostat_val) xD
    read(*,*,iostat=iostat_val) xB
    read(*,*,iostat=iostat_val) zF
    read(*,*,iostat=iostat_val) q
    read(*,*,iostat=iostat_val) Rfactor
    read(*,*,iostat=iostat_val) F
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all distillation inputs.'
        stop
    end if
    if (alpha <= 1.0d0) then
        write(*,*) 'ERROR: Relative volatility alpha must be greater than 1.'
        stop
    end if
    if (.not.(xB > 0.0d0 .and. xB < zF .and. zF < xD .and. xD < 1.0d0)) then
        write(*,*) 'ERROR: Require 0 < xB < zF < xD < 1.'
        stop
    end if
    if (Rfactor <= 1.0d0) then
        write(*,*) 'ERROR: R/Rmin must be greater than 1.'
        stop
    end if

    eps = 1.0d-10
    if (abs(q-1.0d0) < 1.0d-8) then
        xq = zF
        yq = y_eq(xq, alpha)
    else
        call find_q_intersection(alpha,zF,q,xq,yq)
    end if
    Rmin = (xD - yq)/max(eps, yq - xq)
    if (Rmin < 0.05d0) Rmin = 0.05d0
    R = Rfactor*Rmin
    mR = R/(R+1.0d0)
    bR = xD/(R+1.0d0)
    mS = (yq - xB)/max(eps, xq - xB)
    bS = yq - mS*xq
    Nmin = log((xD/(1.0d0-xD))*((1.0d0-xB)/xB))/log(alpha)
    D = F*(zF-xB)/(xD-xB)
    B = F-D

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   DISTILLATION COLUMN - MCCABE-THIELE ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,F12.5)') '  Relative Volatility alpha = ', alpha
    write(*,'(A,F12.5)') '  Distillate xD             = ', xD
    write(*,'(A,F12.5)') '  Bottoms xB                = ', xB
    write(*,'(A,F12.5)') '  Feed zF                   = ', zF
    write(*,'(A,F12.5)') '  Feed q                    = ', q
    write(*,'(A,ES12.4,A)') '  Feed Flow                 = ', F, ' mol/s'
    write(*,*)
    write(*,'(A)') '--- REFLUX / PINCH RESULTS ----------------------------------'
    write(*,'(A,F12.5)') '  q-line Pinch xq           = ', xq
    write(*,'(A,F12.5)') '  q-line Pinch yq           = ', yq
    write(*,'(A,ES12.4)') '  Minimum Reflux Ratio      = ', Rmin
    write(*,'(A,ES12.4)') '  Actual Reflux Ratio       = ', R
    write(*,'(A,ES12.4)') '  Fenske Minimum Stages     = ', Nmin
    write(*,*)

    write(*,'(A)') '--- VLE CURVE -----------------------------------------------'
    write(*,'(A)') '  x             y_eq'
    do i=0,100
        x = dble(i)/100.0d0
        write(*,'(F10.6,2X,F10.6)') x, y_eq(x,alpha)
    end do
    write(*,*)

    write(*,'(A)') '--- STAGE STEPPING ------------------------------------------'
    write(*,'(A)') '  stage         x             y             section'
    x = xD
    y = xD
    n = 0
    feedStage = 0
    maxStages = 200
    do while (x > xB .and. n < maxStages)
        n = n + 1
        xeq = x_from_y(y, alpha)
        if (xeq >= xq) then
            ynew = mR*xeq + bR
            section = 'RECT'
        else
            if (feedStage == 0) feedStage = n
            ynew = mS*xeq + bS
            section = 'STRIP'
        end if
        write(*,'(I8,2X,F12.6,2X,F12.6,2X,A)') n, xeq, ynew, trim(section)
        xprev = x
        x = xeq
        y = ynew
        if (abs(x-xprev) < 1.0d-12) exit
    end do
    if (feedStage == 0) feedStage = n
    write(*,*)
    write(*,'(A)') '--- DESIGN RESULTS ------------------------------------------'
    write(*,'(A,ES12.4)') '  Number of Stages          = ', dble(n)
    write(*,'(A,I8)')      '  Feed Plate Location       = ', feedStage
    write(*,'(A,ES12.4)') '  Minimum Reflux Ratio      = ', Rmin
    write(*,'(A,ES12.4)') '  Actual Reflux Ratio       = ', R
    write(*,'(A,ES12.4)') '  Fenske Minimum Stages     = ', Nmin
    write(*,'(A,ES12.4,A)') '  Distillate Flow           = ', D, ' mol/s'
    write(*,'(A,ES12.4,A)') '  Bottoms Flow              = ', B, ' mol/s'
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  VLE: y = alpha*x/(1+(alpha-1)*x).'
    write(*,'(A)') '  Fenske equation for total-reflux minimum stages.'
    write(*,'(A)') '  McCabe-Thiele stepping between operating lines and equilibrium curve.'

contains
    double precision function y_eq(x, a)
        double precision, intent(in) :: x, a
        y_eq = a*x/(1.0d0 + (a-1.0d0)*x)
    end function y_eq
    double precision function x_from_y(y, a)
        double precision, intent(in) :: y, a
        x_from_y = y/(a - y*(a-1.0d0))
    end function x_from_y
    subroutine find_q_intersection(a,z,q,xq,yq)
        double precision, intent(in) :: a,z,q
        double precision, intent(out) :: xq,yq
        integer :: it
        double precision :: lo,hi,mid,fmid, yline
        lo=1.0d-9; hi=0.999999d0
        do it=1,120
            mid=0.5d0*(lo+hi)
            yline = q/(q-1.0d0)*mid - z/(q-1.0d0)
            fmid = y_eq(mid,a) - yline
            if (fmid > 0.0d0) then
                lo=mid
            else
                hi=mid
            end if
        end do
        xq=0.5d0*(lo+hi)
        yq=y_eq(xq,a)
    end subroutine find_q_intersection
end program distillation_mccabe_thiele