function Distillation_sfcn2(block)
% This is a level 2 sfcn in simulink describing the production of biodiesel
% through reactive distillation 
% The compounds are always listed in the following order: metOH/TG/DG/MG/GL/BD


%% Establishing the number of inputs and outputs:
nu = 6;
ny = 7;
Ntrays = 20;
NC = 6;

% block definitions
block.NumInputPorts  = nu;
block.NumOutputPorts = ny;

% Setup port properties to be inherited or dynamic
block.SetPreCompInpPortInfoToDynamic;
block.SetPreCompOutPortInfoToDynamic;

for i = 1:nu
    block.InputPort(i).Dimensions = 1;
    block.InputPort(i).DirectFeedthrough = true;
end
block.InputPort(2).Dimensions = NC; %The second input port is the feed molar fraction


for i = 1:ny
    block.OutputPort(i).Dimensions = 1;
end

block.OutputPort(1).Dimensions = NC;
block.OutputPort(2).Dimensions = NC;
block.OutputPort(3).Dimensions = (Ntrays+2)*(NC-1);
block.OutputPort(4).Dimensions = Ntrays+2;
block.OutputPort(5).Dimensions = Ntrays+2;
block.OutputPort(6).Dimensions = Ntrays+1;
block.OutputPort(7).Dimensions = Ntrays+1;

block.SampleTimes = [0 0];

% Setup Dwork
block.NumContStates = (Ntrays+2)*(NC-1)+Ntrays+2;

block.RegBlockMethod('InitializeConditions',    @InitializeConditions);
block.RegBlockMethod('Outputs',                 @Outputs);
block.RegBlockMethod('Derivatives',             @Derivatives);
block.RegBlockMethod('SetInputPortSamplingMode',@SetInpPortFrameData);


function SetInpPortFrameData(block, idx, fd)
  
block.InputPort(idx).SamplingMode = fd;
for i = 1:block.NumOutputPorts
    block.OutputPort(i).SamplingMode = fd;
end


function InitializeConditions(block)

load steadystate.mat
block.ContStates.Data = steadystate;

function Outputs(block)
global l XU
% Important Parameters
Ntrays = 20;                        %Number of trays in column
NC = 6;                             %Number of compounds
Mw = [32.042 878 616 354 92.095 294]; %Molecular weights [kg/kmol]
alpha = [1.45*10^3 1 2.35*10 3.06*10 4.48*10 2.81*10];
Dc = 1.500;                         %Diameter of column [m]
lw = 0.77*Dc;                       %Length of weir [m]
hw = 0.20;                          %Heigth of weir [m]
q = 1.028;                          %Fraction of liquid in feed

%% Key Variables
% State variables
x = block.ContStates.Data(1:(Ntrays+2)*(NC-1)); %Molar fractions
H = block.ContStates.Data(((Ntrays+2)*(NC-1)+1):NC*(Ntrays+2)); %Molar holdup
xL = reshape(x,Ntrays+2,5);  %Reshaping into a matrix
HL = reshape(H,Ntrays+2,1);  %Reshaping into a vector
%We still use the sequence: MetOH/TG/DG/MG/GL/BD

% Disturbances
F0   = block.InputPort(1).Data;
x_in  = block.InputPort(2).Data;

% Manipulation
D   = block.InputPort(3).Data;
Vb  = block.InputPort(4).Data;
Ref = block.InputPort(5).Data;
B   = block.InputPort(6).Data;

%% Calculation of compositions and temperature by method 1
%Setting the initial temperature
Te = 364.15*ones(Ntrays+2,1); 
Te(1,1) = 337.63; Te(Ntrays+2,1) = 397.15; 
XX = zeros(Ntrays+2,NC);  x_BD = zeros(Ntrays+2,1);
y = zeros(Ntrays+2,NC-1); yy = zeros(Ntrays+2,NC);

for i = 1:Ntrays+2
    %Molar fraction of BD:
    x_BD(i,1) = 1 - sum(xL(i,:));
    XX(i,:) = [xL(i,:) x_BD(i,1)];
    
    %Vapour composition by relative volatilities:
    yy(i,:) = (alpha.*XX(i,:))/sum(alpha.*XX(i,:));
    y(i,1) = yy(1,1); y(i,2) = yy(1,2); y(i,3) = yy(1,3); y(i,4) = yy(1,4); y(i,5) = yy(1,5);  
end

%% Calculations of densities in the liquid phase:
%The constants
%       Ar        Br    Cr     Dr
AR = [2.3267   0.27073 512.5 0.24713    % MetOH
      0.026085 0.14259 934.6 0.28571    % TG
      0.92382  0.24386 850   0.22114    % GL
      0.20469  0.23737 767.4 0.28571];  % BD
 
i = 1:Ntrays+2;
rho_MetOH = AR(1,1)./(AR(1,2).^(1+(1-(Te(i,1)/AR(1,3))).^AR(1,4)));
rho_TG = AR(2,1)./(AR(2,2).^(1+(1-(Te(i,1)/AR(2,3))).^AR(2,4)));
rho_GL = AR(3,1)./(AR(3,2).^(1+(1-(Te(i,1)/AR(3,3))).^AR(3,4)));
rho_BD = AR(4,1)./(AR(4,2).^(1+(1-(Te(i,1)/AR(4,3))).^AR(4,4)));
rho_DG = rho_TG;
rho_MG = rho_TG;
rho = [rho_MetOH rho_TG rho_DG rho_MG rho_GL rho_BD];

% kmol of components per tray [kmol]:
N(:,1) = HL.*XX(:,1);   % MetOH
N(:,2) = HL.*XX(:,2);   % TG
N(:,3) = HL.*XX(:,3);   % DG
N(:,4) = HL.*XX(:,4);   % MG
N(:,5) = HL.*XX(:,5);   % GL
N(:,6) = HL.*XX(:,6);   % BD

V = N./rho;             % volume of component per tray [m3]
Vtot = sum(V,2);        % total volume per tray [m3]

% Volume fraction [m3/m3]
ve(:,1) = V(:,1)./Vtot(:,1);    % MetOH
ve(:,2) = V(:,2)./Vtot(:,1);    % TG
ve(:,3) = V(:,3)./Vtot(:,1);    % DG
ve(:,4) = V(:,4)./Vtot(:,1);    % MG
ve(:,5) = V(:,5)./Vtot(:,1);    % GL
ve(:,6) = V(:,6)./Vtot(:,1);    % BD

% Common liquid density for the tray [kmol/m3]
rho_l = sum(rho.*ve,2);   % common liquid density for the tray [kmol/m3]

% Weights per tray [kg]:
wt(:,1) = XX(:,1).*HL*Mw(1,1);  % MetOH
wt(:,2) = XX(:,2).*HL*Mw(1,2);  % TG
wt(:,3) = XX(:,3).*HL*Mw(1,3);  % DG
wt(:,4) = XX(:,4).*HL*Mw(1,4);  % MG
wt(:,5) = XX(:,5).*HL*Mw(1,5);  % GL
wt(:,6) = XX(:,6).*HL*Mw(1,6);  % BD
wt_tot = sum(wt,2);

% Weight Fractions
x_wt(:,1) = wt(:,1)./wt_tot;    % MetOH
x_wt(:,2) = wt(:,2)./wt_tot;    % TG
x_wt(:,3) = wt(:,3)./wt_tot;    % DG
x_wt(:,4) = wt(:,4)./wt_tot;    % MG
x_wt(:,5) = wt(:,5)./wt_tot;    % GL
x_wt(:,6) = wt(:,6)./wt_tot;    % BD

% A common liquid density for the tray [kg/m3]:
rho_Li(:,1) = rho(:,1).*x_wt(:,1)*Mw(1,1);  % MetOH
rho_Li(:,2) = rho(:,2).*x_wt(:,2)*Mw(1,2);  % TG
rho_Li(:,3) = rho(:,3).*x_wt(:,3)*Mw(1,3);  % DG
rho_Li(:,4) = rho(:,4).*x_wt(:,4)*Mw(1,4);  % MG
rho_Li(:,5) = rho(:,5).*x_wt(:,5)*Mw(1,5);  % GL
rho_Li(:,6) = rho(:,6).*x_wt(:,6)*Mw(1,6);  % BD
rho_L(i,1) = sum(rho_Li,2);    

%% Temperature Estimation

% Initializing sizes
T = Te;

% Temperature Estimation by UNIQUAC:
for l = 1:Ntrays+2
    XU = [xL(:,2) xL(:,3) xL(:,4) x_BD(:,1) xL(:,5) xL(:,1)];
    T(l,1) = fzero('sumY_UNIQUAC2',Te(l,1));        % Equilibrium temperature [K]
end


%% Calculating the flowrates in the column
% Initializing sizes
V = ones(Ntrays+1,1);  %V is a vector of the molar vapor flows in the column
L = ones(Ntrays+1,1);  %L is a vector of the different trays in the column.

% Molar balances over the top and the bottom: 
L(1,1) = Ref;

% Applying Francis weir for the liquid flowrates
how = ones(Ntrays,1); Lw = ones(Ntrays,1);
for i = 1:Ntrays
      how(i,1) = (HL(i+1,1)-(pi*hw*rho_l(i+1,1)*0.88*(Dc/2)^2))/(pi*0.88*rho_l(i+1,1)*(Dc/2)^2);    %Calculation of how
      if how(i,1) <0
          Lw(i+1,1) = 0;
          L(i+1,1) = 0;
      else          
          Lw(i+1,1) = ((how(i,1)*10^3)/750)^1.5*rho_L(i+1,1)*lw;    %Francis weir [kg/s]
          L(i+1,1) = (Lw(i+1,1)*3600)/sum(XX(i+1,:).*Mw);           %Liquid flow [kmol/h]
      end
end

% Vapour flowrates
V(Ntrays+1,1) = Vb;

% Assuming constant molar overflow for vapour flowrates
for i = 2:Ntrays
    V(i,1) = V(Ntrays+1,1);
end

% Correction for the feed:
V(1,1) = V(2,1) + ((1-q)*F0); 
       

%% Determining distillate and bottoms flow
%Determination of the distillate flow
vD = D*XX(1,:);       % partial molar flows (kmol/h)

%Determination of the bottoms flow
vB = B*XX(Ntrays+2,:); % partial molar flows (kmol/h)


%% Output variables
block.OutputPort(1).Data = vD;
block.OutputPort(2).Data = vB;
block.OutputPort(3).Data = x;
block.OutputPort(4).Data = H;
block.OutputPort(5).Data = T;
block.OutputPort(6).Data = L;
block.OutputPort(7).Data = V;


function Derivatives(block)
%% Important parameters
global l XU

Ntrays = 20;                %Number of trays
NC = 6;                     %Number of compounds
R = 8.314472*10^(-3);       %Gas constant [kJ/K*mol]
Mw = [32.042 878 616 354 92.095 294]; %Molecular weights [kg/kmol]
q = 1.028;                  %Liquid feed fraction
Dc = 1.500;                 %Diameter of column [m]
lw = 0.77*Dc;               %Length of weir [m]
hw = 0.20;                  %Heigth of weir [m]
alpha = [1.44*10^3 1 2.35*10 3.06*10 4.48*10 2.81*10];


%% Variables
% State variables
  x = block.ContStates.Data(1:(Ntrays+2)*(NC-1));
  H = block.ContStates.Data(((Ntrays+2)*(NC-1)+1):NC*(Ntrays+2));
  xL = reshape(x,Ntrays+2,NC-1);
  HL = reshape(H,Ntrays+2,1);

% Disturbances
  F0   = block.InputPort(1).Data;
  x_in = block.InputPort(2).Data;
%The input F0 is the flowrate [kmol/h] and x_in is the composition.
%The composition is given in the following order: MetOH/TG/DG/MG/GL/BD
  
% Manipulation 
  D   = block.InputPort(3).Data;
  Vb = block.InputPort(4).Data;
  Ref = block.InputPort(5).Data;
  B = block.InputPort(6).Data;
  
%% Calculation of temperatures and compositions

% Calculating the temperature on each tray:
Te = 364.15*ones(Ntrays+2,1); % Setting the initial temperature
Te(1,1) = 337.63; Te(Ntrays+2,1) = 397.15; % Initial temperatures of reboiler and condenser.
XX = zeros(Ntrays+2,NC);  x_BD = zeros(Ntrays+2,1);
y = zeros(Ntrays+2,NC-1); yy = zeros(Ntrays+2,NC);

for i = 1:Ntrays+2
    x_BD(i,1) = 1 - sum(xL(i,:));
    XX(i,:) = [xL(i,:) x_BD(i,1)];
    
    %Vapour phase composition by relative volatilities:
    yy(i,:) = (alpha.*XX(i,:))/sum(alpha.*XX(i,:));
    y(i,1) = yy(1,1); y(i,2) = yy(1,2); y(i,3) = yy(1,3); y(i,4) = yy(1,4); y(i,5) = yy(1,5); 
end

%% Calculations of densities in the liquid phase:
%The constants
%       Ar        Br    Cr     Dr
AR = [2.3267   0.27073 512.5 0.24713    % MetOH
      0.026085 0.14259 934.6 0.28571    % TG
      0.92382  0.24386 850   0.22114    % GL
      0.20469  0.23737 767.4 0.28571];  % BD
 
i = 1:Ntrays+2;
rho_MetOH = AR(1,1)./(AR(1,2).^(1+(1-(Te(i,1)/AR(1,3))).^AR(1,4)));
rho_TG = AR(2,1)./(AR(2,2).^(1+(1-(Te(i,1)/AR(2,3))).^AR(2,4)));
rho_GL = AR(3,1)./(AR(3,2).^(1+(1-(Te(i,1)/AR(3,3))).^AR(3,4)));
rho_BD = AR(4,1)./(AR(4,2).^(1+(1-(Te(i,1)/AR(4,3))).^AR(4,4)));
rho_DG = rho_TG;
rho_MG = rho_TG;
rho = [rho_MetOH rho_TG rho_DG rho_MG rho_GL rho_BD];

% kmol of components per tray [kmol]:
N(:,1) = HL.*XX(:,1);   % MetOH
N(:,2) = HL.*XX(:,2);   % TG
N(:,3) = HL.*XX(:,3);   % DG
N(:,4) = HL.*XX(:,4);   % MG
N(:,5) = HL.*XX(:,5);   % GL
N(:,6) = HL.*XX(:,6);   % BD

V = N./rho;             % volume of component per tray [m3]
Vtot = sum(V,2);        % total volume per tray [m3]

% Volume fraction [m3/m3]
ve(:,1) = V(:,1)./Vtot(:,1);    % MetOH
ve(:,2) = V(:,2)./Vtot(:,1);    % TG
ve(:,3) = V(:,3)./Vtot(:,1);    % DG
ve(:,4) = V(:,4)./Vtot(:,1);    % MG
ve(:,5) = V(:,5)./Vtot(:,1);    % GL
ve(:,6) = V(:,6)./Vtot(:,1);    % BD

% Common liquid density for the tray [kmol/m3]
rho_l = sum(rho.*ve,2);   % common liquid density for the tray [kmol/m3]

% Weights per tray [kg]:
wt(:,1) = XX(:,1).*HL*Mw(1,1);  % MetOH
wt(:,2) = XX(:,2).*HL*Mw(1,2);  % TG
wt(:,3) = XX(:,3).*HL*Mw(1,3);  % DG
wt(:,4) = XX(:,4).*HL*Mw(1,4);  % MG
wt(:,5) = XX(:,5).*HL*Mw(1,5);  % GL
wt(:,6) = XX(:,6).*HL*Mw(1,6);  % BD
wt_tot = sum(wt,2);

% Weight Fractions
x_wt(:,1) = wt(:,1)./wt_tot;    % MetOH
x_wt(:,2) = wt(:,2)./wt_tot;    % TG
x_wt(:,3) = wt(:,3)./wt_tot;    % DG
x_wt(:,4) = wt(:,4)./wt_tot;    % MG
x_wt(:,5) = wt(:,5)./wt_tot;    % GL
x_wt(:,6) = wt(:,6)./wt_tot;    % BD

% A common liquid density for the tray [kg/m3]:
rho_Li(:,1) = rho(:,1).*x_wt(:,1)*Mw(1,1);  % MetOH
rho_Li(:,2) = rho(:,2).*x_wt(:,2)*Mw(1,2);  % TG
rho_Li(:,3) = rho(:,3).*x_wt(:,3)*Mw(1,3);  % DG
rho_Li(:,4) = rho(:,4).*x_wt(:,4)*Mw(1,4);  % MG
rho_Li(:,5) = rho(:,5).*x_wt(:,5)*Mw(1,5);  % GL
rho_Li(:,6) = rho(:,6).*x_wt(:,6)*Mw(1,6);  % BD
rho_L(i,1) = sum(rho_Li,2);    


%% Calculating the flowrates in the column
% Initializing sizes
V = ones(Ntrays+1,1);  %V is a vector of the molar vapor flows in the column
L = ones(Ntrays+1,1);  %L is a vector of the different trays in the column.
% Molar balances over the top: 
L(1,1) = Ref;   %The first liquid molar flow is equal to the reflux.

% Applying Francis weir for the liquid flowrates
how = ones(Ntrays,1); Lw = ones(Ntrays,1);
for i = 2:Ntrays+1
    %Calculation of height over weir [m]:
      how(i,1) = (HL(i,1)-(pi*hw*rho_l(i,1)*0.88*(Dc/2)^2))/(pi*0.88*rho_l(i,1)*(Dc/2)^2); 
      if how(i,1) <0
          %Ensures no negative flowrates within column
          Lw(i,1) = 0;
          L(i,1) = 0;
      else
          %Francis weir [kg/s]
          Lw(i,1) = ((how(i,1)*10^3)/750)^1.5*rho_L(i,1)*lw;    
          %Liquid flowrate [kmol/h]
          L(i,1) = (Lw(i,1)*3600)/sum(XX(i,:).*Mw);           
      end
end

% Vapour flowrates
V(Ntrays+1,1) = Vb;          

%Assume constant molar flowrate:
for i = 2:Ntrays
    V(i,1) = V(Ntrays+1,1);
end

% Correction for the feed:
V(1,1) = V(2,1) + ((1-q)*F0); 

%% Calculating the vapour compositions by vapour pressures + temperature by method 2
T = Te;
%      TG       DG      MG      BD        GL    MetOH  
XU = [xL(:,2) xL(:,3) xL(:,4) x_BD(:,1) xL(:,5) xL(:,1)];

for l = 1:Ntrays+2 
   % Function ''fzero('sumY',value) returns the T value (equilibrium temperature) that gives 1-sumY=0
  T(l,1) = fzero('sumY_UNIQUAC2',Te(l,1));        % Equilibrium temperature [K]
end
 
%% Calculation of the kinetics in the liquid phase:

%The sizes of the kinetic constants:
k1 = ones(Ntrays+2,1); k2 = ones(Ntrays+2,1); k3 = ones(Ntrays+2,1);
k4 = ones(Ntrays+2,1); k5 = ones(Ntrays+2,1); k6 = ones(Ntrays+2,1);

for i = 1:Ntrays+2
%Kinetic data for the reaction, rate constants [m3/kmol*h]
k1(i) = 3.9*3600*(10^7)*exp(-54.9987/(R*T(i,1)));    
k2(i) = 5.78*3600*(10^5)*exp(-41.5555/(R*T(i,1)));
k3(i) = 5.906*3600*(10^12)*exp(-83.0942/(R*T(i,1)));
k4(i) = 9.888*3600*(10^9)*exp(-61.2496/(R*T(i,1)));
k5(i) = 5.335*3600*(10^3)*exp(-26.8655/(R*T(i,1)));
k6(i) = 2.1*3600*(10^4)*exp(-40.1162/(R*T(i,1)));
end

%We need to convert the molar input to concentrations in order to calculate
%the reaction outputs:
C = zeros(Ntrays+2,6); 
R = zeros(Ntrays+2,6);
M = zeros(size(C));
Msum = zeros(Ntrays+2,1);
for i = 1:Ntrays+2
    vH = HL(i,1)*XX(i,:);       % Molar holdup per component
    vr = sum(vH./rho(i,:));       % volume holdup
    C(i,:) = vH/vr;                 % concentration of each component
    % The reactions. R is the gain/loss in [kmol/h]
    R_TG = -(k1(i)*C(i,2)*C(i,1)) + (k2(i)*C(i,3)*C(i,6));
    R_DG = (k1(i)*C(i,2)*C(i,1)) - (k2(i)*C(i,3)*C(i,6)) - (k3(i)*C(i,3)*C(i,1)) + (k4(i)*C(i,4)*C(i,6));
    R_MG = (k3(i)*C(i,3)*C(i,1)) - (k4(i)*C(i,4)*C(i,6)) - (k5(i)*C(i,4)*C(i,1)) + (k6(i)*C(i,5)*C(i,6));
    R_BD = (k1(i)*C(i,2)*C(i,1))-(k2(i)*C(i,3)*C(i,6))+(k3(i)*C(i,3)*C(i,1))-(k4(i)*C(i,4)*C(i,6))+(k5(i)*C(i,4)*C(i,1))-(k6(i)*C(i,5)*C(i,6));
    R_GL = (k5(i)*C(i,4)*C(i,1)) - (k6(i)*C(i,5)*C(i,6));
    R_MetOH = -(k1(i)*C(i,2)*C(i,1))+(k2(i)*C(i,3)*C(i,6))-(k3(i)*C(i,3)*C(i,1))+(k4(i)*C(i,4)*C(i,6))-(k5(i)*C(i,4)*C(i,1))+(k6(i)*C(i,5)*C(i,6));
    % We store the molar gain of each component in an R matrix:
    R(i,1) = R_MetOH; R(i,2) = R_TG; R(i,3) = R_DG; R(i,4) = R_MG; 
    R(i,5) = R_GL; R(i,6) = R_BD;
    M(i,:) = R(i,:)*vr; %The change in molar holdup of each component saved in a vector
    Msum(i,1) = sum(M(i,:));
end

%% Molar balances in the column (the states)
dHxdt = zeros(Ntrays+2,5); 
dHdt = zeros(Ntrays+2,1);

% Column
i = 2:Ntrays+1;
dHdt(i,1) = L(i-1,1) - L(i,1) + V(i,1) - V(i-1,1) + Msum(i); 
% Correction for the feed:
dHdt(2,1) = dHdt(2,1) + F0;

%The molar holdup composition balance for the feed tray:
for j = 1:5;
    for i = 2
        dHxdt(i,j) = ((F0*x_in(j,1)) + (V(i,1)*y(i,j)) + L(i-1,1)*xL(i-1,j) - (V(i-1,1)*y(i-1,j)) - (L(i,1)*xL(i,j)) + M(i,j));
    end
end

%The molar holdup composition balance for the other reactive trays:
for j = 1:5;
    for i = 3:Ntrays+1
        dHxdt(i,j) = (V(i,1)*y(i,j)) + (L(i-1,1)*xL(i-1,j) - (V(i-1,1)*y(i-1,j)) - (L(i,1)*xL(i,j)) + M(i,j));   
    end
end

%Reboiler (assumed to be an equilibrium stage)
i = Ntrays+2;
dHdt(i,1) = L(Ntrays+1,1) - V(Ntrays+1,1) - B;

j=1:5;
dHxdt(i,j)= (L(Ntrays+1,1)*xL(Ntrays+1,j) - V(Ntrays+1,1)*y(Ntrays+2,j) - B*xL(Ntrays+2,j));

%Total condenser (not an equilibrium stage)
i = 1;
dHdt(i,1) = V(1,1) - L(1,1) - D;

j = 1:5;
dHxdt(i,j)= V(1,1)*y(i,j) - L(1,1)*xL(i,j) - D*xL(i,j);

%Computing the derivative for the mole fractions from d(Mx) = x dM + M dx
dxdt=(dHxdt - (xL.*(dHdt*ones(1,5))))./(HL*ones(1,5));

dxdtout = reshape(dxdt,(Ntrays+2)*5,1);
dxdtout = [dxdtout; dHdt];
     
     
  block.Derivatives.Data = dxdtout;

