program liquid_liquid_extraction
    implicit none
    integer :: mode, N, i, reqStages, iostat_val
    double precision :: F, xF, S, yS, KD, targetRecovery, E, solventRatio, residualFrac
    double precision :: xR, yE, recovery, totalSolute, raffSolute, extSolute, extractYavg

    read(*,*,iostat=iostat_val) mode
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid mode input.'
        stop
    end if
    read(*,*,iostat=iostat_val) F
    read(*,*,iostat=iostat_val) xF
    read(*,*,iostat=iostat_val) S
    read(*,*,iostat=iostat_val) yS
    read(*,*,iostat=iostat_val) KD
    read(*,*,iostat=iostat_val) N
    read(*,*,iostat=iostat_val) targetRecovery
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all extraction inputs.'
        stop
    end if
    if (F <= 0.0d0 .or. S < 0.0d0 .or. KD <= 0.0d0) then
        write(*,*) 'ERROR: Feed flow, solvent flow and KD must be valid.'
        stop
    end if
    if (xF < 0.0d0 .or. xF >= 1.0d0 .or. yS < 0.0d0 .or. yS >= 1.0d0) then
        write(*,*) 'ERROR: Compositions must be between 0 and less than 1.'
        stop
    end if
    if (N < 1) N = 1
    if (targetRecovery < 0.0d0 .or. targetRecovery >= 100.0d0) then
        write(*,*) 'ERROR: Target recovery must be between 0 and less than 100 percent.'
        stop
    end if

    totalSolute = F*xF + S*yS
    solventRatio = S/F
    E = KD*S/F

    if (mode == 1) then
        xR = (F*xF + S*yS)/(F + KD*S)
        yE = KD*xR
        raffSolute = F*xR
        extSolute = totalSolute - raffSolute
        recovery = 100.0d0*extSolute/max(1.0d-30, F*xF)
    else if (mode == 2) then
        residualFrac = (1.0d0/(1.0d0 + E/dble(N)))**N
        xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
        yE = KD*xR
        raffSolute = F*xR
        extSolute = max(0.0d0, F*xF - raffSolute)
        recovery = 100.0d0*extSolute/max(1.0d-30, F*xF)
    else
        residualFrac = kremser_residual(E,N)
        xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
        yE = KD*xR
        raffSolute = F*xR
        extSolute = max(0.0d0, F*xF - raffSolute)
        recovery = 100.0d0*extSolute/max(1.0d-30, F*xF)
    end if
    if (S > 0.0d0) then
        extractYavg = (S*yS + extSolute)/S
    else
        extractYavg = 0.0d0
    end if

    reqStages = 1
    do i=1,200
        if (mode == 1) then
            xR = (F*xF + S*yS)/(F + KD*S)
        else if (mode == 2) then
            residualFrac = (1.0d0/(1.0d0 + E/dble(i)))**i
            xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
        else
            residualFrac = kremser_residual(E,i)
            xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
        end if
        recovery = 100.0d0*(F*xF - F*xR)/max(1.0d-30,F*xF)
        if (recovery >= targetRecovery) then
            reqStages = i
            exit
        end if
        reqStages = i
    end do

    ! Recompute specified case after required-stage loop
    if (mode == 1) then
        xR = (F*xF + S*yS)/(F + KD*S)
        yE = KD*xR
    else if (mode == 2) then
        residualFrac = (1.0d0/(1.0d0 + E/dble(N)))**N
        xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
        yE = KD*xR
    else
        residualFrac = kremser_residual(E,N)
        xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
        yE = KD*xR
    end if
    raffSolute = F*xR
    extSolute = max(0.0d0, F*xF - raffSolute)
    recovery = 100.0d0*extSolute/max(1.0d-30,F*xF)
    if (S > 0.0d0) extractYavg = (S*yS + extSolute)/S

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   LIQUID-LIQUID EXTRACTION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,I8)')       '  Mode                     = ', mode
    write(*,'(A,ES12.4)')   '  Feed Flow (F)            = ', F
    write(*,'(A,ES12.4)')   '  Feed Composition x_F     = ', xF
    write(*,'(A,ES12.4)')   '  Solvent Flow (S)         = ', S
    write(*,'(A,ES12.4)')   '  Solvent Composition y_S  = ', yS
    write(*,'(A,ES12.4)')   '  Distribution Coefficient (K_D) = ', KD
    write(*,'(A,I8)')       '  Specified Stages         = ', N
    write(*,*)
    write(*,'(A)') '--- DESIGN RESULTS ------------------------------------------'
    write(*,'(A,I8)')       '  Number of Stages         = ', reqStages
    write(*,'(A,ES12.4)')   '  Solvent-to-Feed Ratio    = ', solventRatio
    write(*,'(A,ES12.4)')   '  Extraction Factor        = ', E
    write(*,'(A,ES12.4,A)') '  Solute Recovery          = ', recovery, ' percent'
    write(*,'(A,ES12.4)')   '  Final Raffinate x_R      = ', xR
    write(*,'(A,ES12.4)')   '  Final Extract y_E        = ', extractYavg
    write(*,'(A,ES12.4)')   '  Equilibrium y*_E         = ', yE
    write(*,'(A,ES12.4)')   '  Solute Extracted         = ', extSolute
    write(*,'(A,ES12.4)')   '  Raffinate Solute         = ', raffSolute
    write(*,*)
    write(*,'(A)') '--- STAGE PROFILE -------------------------------------------'
    write(*,'(A)') '  stage        x_R           y_E           recovery[%]   raff_solute   ext_solute'
    write(*,'(A)') '  -----------------------------------------------------------------------------'
    do i=1,max(N,reqStages)
        if (mode == 1) then
            xR = (F*xF + S*yS)/(F + KD*S)
        else if (mode == 2) then
            residualFrac = (1.0d0/(1.0d0 + E/dble(i)))**i
            xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
        else
            residualFrac = kremser_residual(E,i)
            xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
        end if
        yE = KD*xR
        raffSolute = F*xR
        extSolute = max(0.0d0,F*xF-raffSolute)
        recovery = 100.0d0*extSolute/max(1.0d-30,F*xF)
        if (S > 0.0d0) then
            extractYavg = (S*yS + extSolute)/S
        else
            extractYavg = 0.0d0
        end if
        write(*,'(I8,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
            i, xR, extractYavg, recovery, raffSolute, extSolute
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Distribution equilibrium: y* = K_D x.'
    write(*,'(A)') '  Single stage balance: x1 = (F xF + S yS)/(F + KD S).'
    write(*,'(A)') '  Crosscurrent approximation uses equal solvent split across stages.'
    write(*,'(A)') '  Countercurrent Kremser/Hunter-Nash style residual fraction = (E-1)/(E^(N+1)-1).'

contains
    double precision function kremser_residual(Ein,Nin)
        implicit none
        double precision, intent(in) :: Ein
        integer, intent(in) :: Nin
        if (abs(Ein-1.0d0) < 1.0d-10) then
            kremser_residual = 1.0d0/dble(Nin+1)
        else if (Ein > 0.0d0) then
            kremser_residual = (Ein-1.0d0)/(Ein**dble(Nin+1)-1.0d0)
            if (kremser_residual < 0.0d0) kremser_residual = abs(kremser_residual)
            if (kremser_residual > 1.0d0) kremser_residual = 1.0d0
        else
            kremser_residual = 1.0d0
        end if
    end function kremser_residual
end program liquid_liquid_extraction
