% The distillation model

% Liquid control in the refulx drum and reboiler: flow equations for D and B
switch config
    case('LV')% This flag signals the inclusion of equations for liquid levels. If flag=1, Simulink is being used and no liquid level is calculated in the code.
    D = Ds + (M(NT) - MDs)*KcD;      % Distillate flow
    B = Bs + (M(1) - MBs)*KcB;       % Bottoms flow
    confflag = 0 ;
    case('BV')
    D = Ds + (M(NT) - MDs)*KcD;      % Distillate flow
    LT = Bs + (M(1) - MBs)*KcB;      % Reflux flow
    confflag = 0 ;
    case('XdXb')
    D = Ds + (M(NT) - MDs)*KcD;      % Distillate flow
    B = Bs + (M(1) - MBs)*KcB;       % Bottoms flow
    if number == 1 
        LT = LT0_1 - KcL1*e(1) - (KcL1/tauiL1)*eint(1);       % Reflux flow
        VB = VB0_1 - KcV1*e(2) - (KcV1/tauiV1)*eint(2);       % Boilup
    elseif number == 2
        LT = LT0_2 - KcL2*e(1) - (KcL2/tauiL2)*eint(1);       % Reflux flow
        VB = VB0_2 - KcV2*e(2) - (KcV2/tauiV2)*eint(2);       % Boilup
    end 
    confflag = 1 ;
end

% Vapour-liquid equilibria (multicomponent ideal VLE, Stichlmair-Fair, 'Distillation', p. 36, 1998)
y = (x*diag(alpha(1:NC-1)))./((x*(alpha(1:NC-1) - 1)' + 1)*ones(1,NC-1));

% Vapor flows assuming constant molar flows
V = VB*ones(NT-1,1); V(NF:NT-1) = V(NF:NT-1) + (1-qF)*F;

% Liquid flows are given by Franci's Weir Formula L(i)=K*Mow(i)^1.5 
% Liquid flow L(i) dependent only on the holdup over the weir Mow(i)
% M(i)= Mow(i) + Muw(i) (Total holdup = holdup over weir + holdup below weir) 
L(2:NF,1)      = Kbf*(M(2:NF) - Muw).^1.5;         % Liquid flows below feed (Kmol/min)
L(NF+1:NT-1,1) = Kuf*(M(NF+1:NT-1) - Muw).^1.5;    % Liquid flows above feed (Kmol/min)
L(NT,1)        = LT;                               % Condenser's liquid flow (Kmol/min)

% Time derivatives from material balances for: (1) total holdup and (2) component holdup
% Column
j = 2:NT-1; dMdt(j,1)  = L(j+1) - L(j) + V(j-1) - V(j);
i = 1:NC-1; dMxdt(j,i) = diag(L(j+1))*x(j+1,i) - diag(L(j))*x(j,i) + diag(V(j-1))*y(j-1,i) - diag(V(j))*y(j,i);

% Correction for feed at the feed stage: The feed is assumed to be mixed into the feed stage
dMdt(NF)    = dMdt(NF) + F;
dMxdt(NF,:) = dMxdt(NF,:) + F*zF(1,1:NC-1);

% Reboiler (assumed to be an equilibrium stage)
dMdt(1)     = L(2) - V(1) - B;
dMxdt(1,i)  = L(2)*x(2,i) - V(1)*y(1,i) - B*x(1,i);

% Total condenser (no equilibrium stage)
dMdt(NT)    = V(NT-1) - LT - D;
dMxdt(NT,i) = V(NT-1)*y(NT-1,i) - L(NT)*x(NT,i) - D*x(NT,i);

% Compute the derivative for the mole fractions from d(Mx) = x dM + M dx
%dxdt = inv(diag(M))*(dMxdt - diag(dMdt)*x);
dxdt =(diag(M))\(dMxdt - diag(dMdt)*x);

% Rearrange elements of composition vector derivatives (dxdt)
dxdt = dxdt(:);