function mct_main(xD, xF, xW, q, R, VLE, isInfo)
    // Scilab implementation of McCabe-Thiele method
    //
    //  xD     - distillate mole fraction              (scalar)
    //  xF     - feed mole fraction                    (scalar)
    //  xW     - bottoms mole fraction                 (scalar)
    //  q      - mole fraction of liquid in the feed   (scalar)
    //  R      - reflux ratio                          (scalar)
    //  VLE    - vapour-liquid data                    (2 x N matrix)
    //           VLE(1,:) = x values
    //           VLE(2,:) = y values
    //  isInfo - 0 = no info printed
    //  `        1 = calculated values printed
    //           2 = 1 + fsolve informations printed (if you have problems,
    //               control the termination indicator of fsolve)
    //
    // ABOUT
    // McCabe-Thiele method (rev.4)
    // author:  Jakub Kopáč
    // e-mail:  kopac.jakub@gmail.com
    //
    // Changes in revision 4
    // - improved detection of nonfeasible definition of reflux ratio
    // - improved detection of incorrect input data

    // ######### DEFINITION OF SUBFUNCTIONS #########
    function abort_for_error(str)
        mprintf("\n  McCabe-Thiele error: %s", str)
        mprintf("\n  # PROGRAM ABORTED #\n")
        abort
    endfunction
    function xiyi = sub_minreflux (xy)
        //find minimal reflux as 2 NAEs for intersection point of q-line and equilibrium line
        x       = xy(1)
        y       = xy(2)
        yeq    = interpln(VLE, x)
        xiyi(1) = ( y*(q-1) - (q*x-xF) )   //q-line
        xiyi(2) = y-yeq                    //equlibrium line
    endfunction
    function xiyi = sub_intersection(xy)
        x       = xy(1)
        y       = xy(2)
        xiyi(1) = y - ( x*R/(R+1) + xD/(R+1) ) //operating line for striping section
        xiyi(2) = y*(q-1) + xF - x*q           //q-line
    endfunction
    function fx = sub_xstep(x, y)
        //hladam x-ovu suradnicu pre vodorovnu ciaru
        fx = y - interpln(VLE, x)
    endfunction
    function fy = sub_ystep(y, x)
        //hladam y-ovu suradnicu pre zvilsu ciaru
        fy = y - interpln([x_work_lines; y_work_lines], x)
    endfunction
    function sub_graphics(xeq,yeq,xwl,ywl,steps,n_steps)
        drawlater()
        plot([xF xyRmin(1)], [xF xyRmin(2)])  //q-line
        plot(xwl, ywl, "r")         //working lines
        plot([xW xW], [0 xW], "k")  //vertical line at xW
        plot([xF xF], [0 xF], "k")  //vertical line at xF
        plot([xD xD], [0 xD], "k")  //vertical line at xD
        plot([0 1], [0 1], "k")     //diagonal line
        plot(xeq, yeq, "b")         //equilibrium line
        plot(steps(:,1), steps(:,2), "m")
        title(msprintf("Number of theoretical plates = %g",n_steps),"fontsize",3)
        xlabel("$x$", "fontsize", 5)
        ylabel("$y$", "fontsize", 5)
        drawnow()
    endfunction
    function x = sub_bisekcia(fun, funpar)
        a   = 0
        b   = 1
        tol = 1e-8
        fa  = fun(a, funpar)
        fb  = fun(b, funpar)

        if (fa*fb)>0 then
            abort_for_error("Unexpected error in bisection algorithm")
        end
        while (b-a) > tol
            x  = (a+b)/2
            fx = fun(x, funpar)
            if fx==0 then
                return
            end
            if fa*fx<0 then
                b  = x
            else
                a  = x
                fa = fx
            end
        end
    endfunction
    function [stepspoints,steps] = sub_stepskernel()
        steps_limit = 300
        x(1)  = xD
        y(1)  = xD
        ii    = 1
        steps = 0
        while x($)>xW
            ii    = ii+1
            x(ii) = sub_bisekcia(sub_xstep, y(ii-1))
            y(ii) = y(ii-1)

            ii    = ii+1
            x(ii) = x(ii-1)
            y(ii) = sub_bisekcia(sub_ystep, x(ii-1))
            steps = steps+1

            if steps>steps_limit then
                mprintf("\nWARNING: Current number of plates is %g",steps)
                WhatToDo = input("Do you want to continue ? (0/1)..")
                if WhatToDo==0 then
                    mprintf("\nCalculation CANCELLED - actual results will be ploted")
                    stepspoints = [x,y]
                    return
                elseif WhatToDo==1
                    steps_limit = steps_limit+100
                end
            end
        end

        y($) = 0
        stepspoints = [x,y]
    endfunction
    function check_scalar(x, xName)
        if ~isscalar(x)
            abort_for_error(xName+" must be scalar")
        end
    endfunction

    // ######### INPUT CHECK #########
    check_scalar(xD, "xD")
    check_scalar(xF, "xF")
    check_scalar(xW, "xW")
    check_scalar(q, "q")
    check_scalar(R, "R")
    check_scalar(isInfo, "isInfo")
    if or([xD<xF, xD<xW, xF<xW]) then
        if xD<xF then
            mprintf("\n  xD can not be lower than xF")
        end
        if xD<xW then
            mprintf("\n  xD can not be lower than xW")
        end
        if xF<xW then
            mprintf("\n  xF can not be lower than xW")
        end
        abort_for_error("Incorrect definition of mole fractions")
    end
    if or([xW<=0, xW>=1]) then
        mprintf("\n  xW must be in range (0,1)")
        abort_for_error("Incorrect definition of mole fractions")
    end
    if or([xF<=0, xF>=1]) then
        mprintf("\n  xF must be in range (0,1)")
        abort_for_error("Incorrect definition of mole fractions")
    end
    if or([xD<=0, xD>=1]) then
        mprintf("\n  xD must be in range (0,1)")
        abort_for_error("Incorrect definition of mole fractions")
    end
    nr = size(VLE, "r");
    if nr<>2 then
        abort_for_error("VLE data wrong size - it must be 2 x N matrix")
    end
    if interpln(VLE, xD)<xD then
        abort_for_error("Equilibrium y for xD is lower than xD")
    end
    if or([VLE(1,:)<0, VLE(1,:)>1]) then
        abort_for_error("x data must be in range [0,1]")
    end
    if or([VLE(2,:)<0, VLE(2,:)>1]) then
        abort_for_error("y data must be in range [0,1]")
    end
    dxeq = diff(VLE(1,:))
    if or(dxeq<=0) then
        abort_for_error("x data must be strictly increasing")
    end
    // ######### INPUT CHECK end #########

    //find minimal reflux ratio intersection points
    [xyRmin,fx1,info1] = fsolve([xF; xF], sub_minreflux)
    //calculate minimal reflux from intersection point (I)
    // R_min = (xD-yI)/(yI-xI)
    R_min = ( xD-xyRmin(2) ) / ( xyRmin(2)-xyRmin(1) )

    // detect infeasible Rmin calculation and update Rmin if necessary
    imin = find(VLE(1,:) <= xyRmin(1))
    imax = find(VLE(1,:) >= xD)
    data = (imin($)+1):1:(imax(1)-1)
    xVLE = VLE(1, data)
    yVLE = VLE(2, data)
    for ii=1:size(xVLE, "*")
        temp       = (xD-yVLE(ii)) / (xD-xVLE(ii))
        R_min2(ii) = temp / (1 - temp)
    end
    R_min = max([R_min, max(R_min2)])

    if R<=R_min then
        mprintf("\n  Calculated R_min = %g", R_min)
        abort_for_error("Reflux lower than calculated minimal reflux")
    end
    //working lines intersection
    [xyInter, fx2, info2] = fsolve([xF; xF], sub_intersection)
    //define x-y coordinates for working lines
    x_work_lines = [0 xW xyInter(1) xD 1]
    y_work_lines = [0 xW xyInter(2) xD 1]

    [stepspoints, n_steps] = sub_stepskernel()

    sub_graphics(VLE(1,:),VLE(2,:),x_work_lines,y_work_lines,stepspoints,n_steps)

    if isInfo>1 then
        mprintf("\n *** fsolve informations ***\n")
        mprintf("\nsub_minreflux")
        mprintf("\n   fx(1) = %1.2e\n   fx(2) = %1.2e",fx1(1),fx1(2))
        if info1==1 then
            mprintf("\n   fsolve termination indicator = %d (OK)\n",info1)
        else
            mprintf("\n   fsolve termination indicator = %d (PROBLEM)\n",info1)
        end
        mprintf("\nsub_intersection")
        mprintf("\n   fx(1) = %1.2e\n   fx(2) = %1.2e",fx2(1),fx2(2))
        if info2==1 then
            mprintf("\n   fsolve termination indicator = %d (OK)",info2)
        else
            mprintf("\n   fsolve termination indicator = %d (PROBLEM)",info2)
        end
        mprintf("\n\n")
    end
    if isInfo>0 then
        xW_calculated = stepspoints($,1)
        nw_nf_ratio   = (xD-xF) / (xD-xW_calculated)
        mprintf("\n *** results ***\n")
        mprintf("\n Calculated minimal reflux ratio = %1.3f\n", R_min)
        mprintf("\n  R / R_min = %1.3f", R/R_min)
        mprintf("\n    nW / nF = %7.5f", nw_nf_ratio)
        mprintf("\n    nD / nF = %7.5f", 1-nw_nf_ratio)
        mprintf("\n         xD = %7.4e", xD)
        mprintf("\n         xW = %7.4e", xW_calculated)
    end
endfunction
