function [func,dae,d_c2,lb,ub,n,Cineq] = RecLoop_Ver2(t,dec,d_nom,KS,KM,H,H_Ver)

% This function defines the differential-algebraic system of equation DAE
% for the ammonia reactor defined in Morud 1998. As addition, a simple 
% recycle loop was incorporated to investigate the influence of the recycle
% on the SOC variables.
% It is a stand alone code which requires CASADI to be loaded beforehand.
%
% The disturbance is for simplicity incorporated into the integrator
% function as a parameter.
%
% If desired, the controller parameters can be directly incorporated. It is
% possible to operate with only a slave controller or only a master
% controller.
%
% Ver2 uses the mass fraction of ammonia and nitrogen as independent variables.
% It requires the following files:
%       MassFlow_Ver2.m
%       Rrate_NewDist_Ver2.m
%
% The algebraic variables are:
%       wNH3:       Ammonia mass fraction [-]
%       wN2:        Nitrogen mass fraction [-]
%       in:         Inlet variables
%
% The input is given as
%       t:          Time step for the integrator
%
%       dec:        Decision variable structure containing the following fields
%                   .Loop   True/False for a closed loop system
%                   .dp     True/False for a pressure drop in the system
%
%       d_nom:      Nominal inlet conditions
%       KS:         Slave controller tuning parameters (not needed)
%       KM:         Master controller tuning parameters (not needed)
%       H:          Structure containing the different H matrices (not
%                   needed)
%       H_Ver:      Integer for the selected H matrix (not needed)
%
% The output is given as
%       func:       Integrator function for CASADI
%
%       dae:        Structure containing the dae with the following fields
%                   .x      Differential variables (Temperatures [C])
%                   .z      Algebraic variables (Reactor and 
%                           inlet concentrations [-], and inlet temperatures [C])
%                   .p      Parameters (1:3, Split ratios to reactor beds)
%                           and disturbances (4:10)
%                   .ode    Differential equations for the temperatures
%                   .alg    Algebraic equations for the ceoncentrations and
%                           inlet temperatures
%                   .quad   Cost function (Ammonia extent of reaction [kg/s]
%
%       d:          Variables for the input disturbance, necessary for
%                   substitution
%
%       lb:         Structure containing the lower bounds with the
%                   following fields
%                   .x      Differential variables
%                   .z      Algebraic variables
%                   .u      Parameters
%
%       ub:         Structure containing the upper bounds with the
%                   following fields
%                   .x      Differential variables
%                   .z      Algebraic variables
%                   .u      Parameters
%
%       n:          Structure containing the number of variables with the
%                   following fields
%                   .x      Differential variables
%                   .z      Algebraic variables
%                   .y      Measurement variables
%                   .d      Disturbance variables
%                   .u      Manipulated variables
%                   .p      Parameters (MV+DV)
%
%       Cineq:      Inequality constraint, if optimization should be
%                   performed.

% Definition of the current path folder
FileName = mfilename('fullpath');
[directory,~,~] = fileparts(FileName);
[dirmain,~,~] = fileparts(directory);

% Load of the data
load([dirmain '\Data\Init_opt.mat'])

% Import of CASADI
import casadi.*

if nargin == 6
    H_Ver = 1;
    H_Sel = H;
elseif nargin == 7
    H_Sel = H{H_Ver};
end
    

%% Definition of the parameters of the model

% Definition of the discretization on each reactor
N1 = 10;
N2 = 10;
N3 = 10;
N = N1+N2+N3;

% Definition of the volume of each reactor [m3]
V1 = 6.69;
V2 = 9.63;
V3 = 15.2;

% Definition of the bed mass of the catalyst in the 3 reactors
dm1 = V1/N1*2200;
dm2 = V2/N2*2200;
dm3 = V3/N3*2200;

% Definition of parameters used
par.Cp = 3500;              % Gas heat capacity [J/kg/K]
par.Cpc = 1100;             % Catalyst heat capacity [J/kg/K]
par.dHrx = 2.7e6;           % Heat of Reaction [kJ/kg] Note: 111370 * 28 = 3.118 e6
par.UA   = 151688;          % UA value of the heat exchanger [W/K]

% Definition of the reaction rate values
kval.Eafor = 87090;         % Forward activation energy [J/mol]
kval.Eabac = 198464;        % Backward activation energy [J/mol]
kval.Afor  = 1.79e+4;       % Forward preexp factor
kval.Abac  = 2.57e+16;      % Backward preexp factor
kval.R     = 8.31;          % Ideal gas constant

% Definition of variables for the pressure drop
pd.kd = 1.5;                % Pressure drop coefficient [kmol/s/sqrt(bar)]
pd.n  = 1.5;                % Polytropic factor [-]
pd.mu = 0.8;                % Compression Efficiency [-]

% Definition of the system states
T = MX.sym('T',N);          % Temperature in each of the discretization units [c]
wNH3 = MX.sym('wNH3',N);  	% Mass fraction of ammonia in each of the discretization units [-]
wN2 = MX.sym('wN2',N);     	% Mass fraction of nitrogen in each of the discretization units [-]

Tin = MX.sym('Tin',3);      % Temperature in the inlet flow of the beds [C]
wNH3in = MX.sym('wNH3in',3);% Mass fraction of ammonia in the inlet flow of the beds [-]
wN2in = MX.sym('wN2in',3);  % Mass fraction of nitrogen in the inlet flow of the beds [-]
TexHEx = MX.sym('TexHEx',1);% Outlet temperature of the Heat Exchanger

Rea = MX.sym('Rea',5);      % Inlet variabels reactor (m, p, T, wNH3, wN2) 
Rec = MX.sym('Rec',5);      % Variabels Recycle pre split (m, p, T, wNH3, wN2)
W   = MX.sym('W',1);        % Compressor Duty [kW]

d = MX.sym('d',6);          % Disturbances: Inlet values (m, p, T, wNH3), catalyst activity, and NH3 recovery
d_wN2 = MX.sym('d_wN2',1);  % Disturbance/Optimization variable inlet value wN2
d_c1 = vertcat(d(1:4), d_wN2, d(5:6));
u = MX.sym('u',4);          % Split ratio to Reac1, Reac2, Reac3, and Purge

eSint = MX.sym('eSint',4);  % Integrated error of the slave
cSset = MX.sym('TSset',4);  % Setpoint of the inlet temperatures of each bed

eMint = MX.sym('eMint',4);  % Integrated error of the master
cMset = MX.sym('TMset',4);  % Setpoint of the master controller
cM    = MX.sym('cM',4);     % Controlled variable of the master controller

%% Setup of the model
% Redefinition of the feed flows
Fin = d_nom(1:5) + d_c1(1:5);

% Adjustment of the mass fractions
[m] = MassFrac(Fin(4),Fin(5));
Fin(5) = m(2);
ndot.Fin = MassFlow_Ver2(Fin(4),Fin(5),Fin);

% Redefinition of massflows
m1 = (1-u(2)-u(3))*Rea(1);
m2 = (1-u(3))*Rea(1);
m3 = Rea(1);
mSHEx = (1-u(1)-u(2)-u(3))*Rea(1);

% Calculation of the reaction rates and flowrates within the system
[r1,ndot.sys1] = Rrate_NewDist_Ver2(T(1:10), wNH3(1:10), wN2(1:10), [m1 Rea(2)],kval,d_c1(6));
[r2,ndot.sys2] = Rrate_NewDist_Ver2(T(11:20),wNH3(11:20),wN2(11:20),[m2 Rea(2)],kval,d_c1(6));
[r3,ndot.sys3] = Rrate_NewDist_Ver2(T(21:30),wNH3(21:30),wN2(21:30),[m3 Rea(2)],kval,d_c1(6));
r = horzcat(r1,r2,r3);
ndot.sys = horzcat(ndot.sys1,ndot.sys2,ndot.sys3);
ndot.Rea = MassFlow_Ver2(Rea(4),Rea(5),Rea);
ndot.Rec = MassFlow_Ver2(Rec(4),Rec(5),Rec);

% Predefinition of the cells
Ceq_NH3 = cell(N,1);
Ceq_N2 = cell(N,1);
dTdt = cell(N,1);
Ceqin = cell(10,1);
Ceq_Rea = cell(5,1);
Ceq_Rec = cell(5,1);

% Calculation of the Heat Exchanger outlet temperature
CStar = (1-u(1)-u(2)-u(3));
NTU   = par.UA/(mSHEx*par.Cp);
eps   = (1-exp(-NTU*(1-CStar)))/(1-CStar*exp(-NTU*(1-CStar)));
Ceqin{4} = TexHEx - (eps*T(N)+(1-eps)*Rea(3));

%% Reactor 1

% Calculation of the inlet variables
Ceqin{1} = m1*Tin(1) - (mSHEx*TexHEx + u(1)*Rea(1)*Rea(3));    % Temperature [K]
Ceqin{5} = wNH3in(1) - Rea(4);                                 % Mass fraction ammonia [-]
Ceqin{8} = wN2in(1)  - Rea(5);                                 % Mass fraction nitrogen [-]

% Bed 1
dTdt{1}    = (m1*par.Cp*(Tin(1)-T(1)) - r(3,1)*dm1*par.dHrx)/(dm1*par.Cpc);
Ceq_NH3{1} = m1*(wNH3in(1)-wNH3(1)) - r(3,1)*dm1;
Ceq_N2{1}  = m1*(wN2in(1)-wN2(1))   - r(2,1)*dm1;

% Bed 2:N1
for j = 2:N1
    dTdt{j}    = (m1*par.Cp*(T(j-1)-T(j)) - r(3,j)*dm1*par.dHrx)/(dm1*par.Cpc);
    Ceq_NH3{j} = m1*(wNH3(j-1)-wNH3(j)) - r(3,j)*dm1;
    Ceq_N2{j}  = m1*(wN2(j-1)-wN2(j))   - r(2,j)*dm1;
end

%% Reactor 2

% Calculation of the inlet variables
Ceqin{2} = m2*Tin(2) - (m1*T(N1) + (u(2)*Rea(1))*Rea(3));       % Temperature [K]
Ceqin{6} = m2*wNH3in(2) - (m1*wNH3(N1) + (u(2)*Rea(1))*Rea(4)); % Mass fraction ammonia [-]
Ceqin{9} = m2*wN2in(2)  - (m1*wN2(N1)  + (u(2)*Rea(1))*Rea(5)); % Mass fraction nitrogen [-];                                         % H2/N2 ratio [-]

% Bed 1
dTdt{N1+1}    = (m2*par.Cp*(Tin(2)-T(N1+1)) - r(3,N1+1)*dm2*par.dHrx)/(dm2*par.Cpc);
Ceq_N2{N1+1}  = m2*(wN2in(2)-wN2(N1+1))   - r(2,N1+1)*dm2;
Ceq_NH3{N1+1} = m2*(wNH3in(2)-wNH3(N1+1)) - r(3,N1+1)*dm2;

% Bed 2:N1
for j = N1+2:N1+N2
    dTdt{j}    = (m2*par.Cp*(T(j-1)-T(j)) - r(3,j)*dm2*par.dHrx)/(dm2*par.Cpc);
    Ceq_N2{j}  = m2*(wN2(j-1)-wN2(j))   - r(2,j)*dm2;
    Ceq_NH3{j} = m2*(wNH3(j-1)-wNH3(j)) - r(3,j)*dm2;
end

%% Reactor 3

% Calculation of the inlet variables
Ceqin{3}  = m3*Tin(3) - (m2*T(N1+N2) + u(3)*Rea(1)*Rea(3));         % Temperature [K]
Ceqin{7}  = m3*wNH3in(3) - (m2*wNH3(N1+N2) + u(3)*Rea(1)*Rea(4));   % Mass fraction [-]
Ceqin{10} = m3*wN2in(3)  - (m2*wN2(N1+N2)  + u(3)*Rea(1)*Rea(5));   % Mass fraction [-]

% Bed 1
dTdt{N1+N2+1}    = (m3*par.Cp*(Tin(3)-T(N1+N2+1)) - r(3,N1+N2+1)*dm3*par.dHrx)/(dm3*par.Cpc);
Ceq_NH3{N1+N2+1} = m3*(wNH3in(3)-wNH3(N1+N2+1)) - r(3,N1+N2+1)*dm3;
Ceq_N2{N1+N2+1}  = m3*(wN2in(3)-wN2(N1+N2+1)) - r(2,N1+N2+1)*dm3;

% Bed 2:N1
for j = N1+N2+2:N
    dTdt{j}    = (m3*par.Cp*(T(j-1)-T(j)) - r(3,j)*dm3*par.dHrx)/(dm3*par.Cpc);
    Ceq_N2{j}  = m3*(wN2(j-1)-wN2(j))   - r(2,j)*dm3;
    Ceq_NH3{j} = m3*(wNH3(j-1)-wNH3(j)) - r(3,j)*dm3;
end

%% Recycle loop

% Reactor inlet variables
Ceq_Rea{1} = Fin(1) + (1-u(4))*Rec(1) - Rea(1);	% Flow rate 
Ceq_Rea{2} = Fin(2) - Rea(2);                   % Pressure drop correlation
Ceq_Rea{3} = Fin(3) - Rea(3);                   % Temperature correlation
Ceq_Rea{4} = Fin(1)*Fin(4) + ...
             Rec(1)*Rec(4)*(1-u(4)) - ...
             Rea(1)*Rea(4);                     % Mass fraction ammonia [-]
Ceq_Rea{5} = Fin(1)*Fin(5) + ...
             Rec(1)*Rec(5)*(1-u(4)) - ...
             Rea(1)*Rea(5);                     % Mass fraction nitrogen [-]

% Recycle variables pre Split
if dec.Loop
    Ceq_Rec{1} = Rea(1)*(1-wNH3(N)*(d_nom(7)+d_c1(7))) - Rec(1);
    Ceq_Rec{2} = pd.kd*sqrt(Rea(2)-Rec(2)) - sum(ndot.Rec);                         % Pressure drop correlation
    Ceq_Rec{3} = (Rea(3)+273.15) - (Rec(3)+273.15)*(Rea(2)/Rec(2))^((pd.n-1)/pd.n); % Temperature correlation
    Ceq_Rec{4} = Rea(1)*wNH3(N)*(1-(d_nom(7)+d_c1(7))) - Rec(1)*Rec(4);
    Ceq_Rec{5} = Rea(1)*wN2(N) - Rec(1)*Rec(5);
else
    Ceq_Rec{1} = Rec(1);
    Ceq_Rec{2} = Rec(2);                        % Pressure drop correlation
    Ceq_Rec{3} = Rec(3);                        % Temperature correlation
    Ceq_Rec{4} = Rec(4);
    Ceq_Rec{5} = Rec(5);
end

% Additional equations required
Ceq_W = W;                                      % Compressor duty

if dec.dp
    Ceq_W      = W - 1/pd.mu*sum(ndot.Rec)*(1-u(4))*(pd.n*kval.R*(Rec(3)+273.15))/(pd.n-1)*...
                     ((Rea(2)/Rec(2))^((pd.n-1)/pd.n)-1);
end

% Definition of the inequality constraint on the input
Cineq = u(1)+u(2)+u(3)-1;

% Definition of the cost function
if dec.Loop
    L = -0.2*Rea(1)*wNH3(N)*(d_nom(7)+d_c1(7)) - 0.01*u(4)*Rec(1) + 0.040*W/3600 + .08*Fin(1);
%     L = sum(ndot.Rec); % Rec(1); % 
%     Cineq = [Cineq;
%              -Rea(1)*wNH3(end)*(d_nom(7)+d_c1(7))+16.16]; % Fin(1)]; % 
else
    L = -(wNH3(end)-wNH3in(1))*Rea(1);
end

%% Combination of the variables of the system
    
% Preassignment of the differential states
z = [wNH3; wN2; Tin; TexHEx; wNH3in; wN2in; Rea; Rec; W];
d_c2 = [d_wN2; d];

% Definition of the measurement vector
y = [z(N+4); z(N+1); T(N1); z(N+2); T(N1+N2); z(N+3); T(N)];

% Rearrangement of differential and algebraic equations
ddt = vertcat(dTdt{:});
Ceq  = vertcat(Ceq_NH3{:}, Ceq_N2{:}, Ceqin{:}, Ceq_Rea{:}, Ceq_Rec{:}, Ceq_W);
        
% Assignment of initital states based on whether Slave and/or Master
% controller are used
if exist('KS','var')
    % Assignment of the differential and algebraic states
    x = [T; eSint];
    z = [z; u];
    p = [cSset; d_c2];
    
    % Calculation of the integral error
    deS_intdt = cSset-Tin;
    
    % Calculation of the controller action
    Ceq_uS = u - max(p0(1:3) + KS.I.*eSint,0.0001);

    % Rearrangement of dTdt and Ceq
    ddt = vertcat(ddt{:}, deS_intdt);
    Ceq = vertcat(Ceq{:}, Ceq_uS);

    % SOC Master controller is used
    if exist('KM','var')
        % Assignment of the differential and algebraic states
        x = [T; eSint; eMint];
        z = [z; cSset; cM];
        p = [cMset; d_c2];
                
        % Calculation of the controlled variables
        Ceq_cM = cM - H_Sel*y;
        
        % Calculation of the integral error
        deM_intdt = cMset-cM;

        % Calculation of the controller action
        Ceq_uM = cSset - (z0(31:33) + KM.p(:,H_Ver).*deM_intdt + KM.I(:,H_Ver).*eMint);

        % Rearrangement of dTdt and Ceq
        ddt = vertcat(ddt{:}, deM_intdt);
        Ceq = vertcat(Ceq{:}, Ceq_uM, Ceq_cM);
    end
    
elseif exist('KM','var')
    % Assignment of the differential and algebraic states
    x = [T; eMint];
    z = [z; u];
    p = [cMset; d_c2];
    
    % Calculation of the controlled variables
    Ceq_cM = cM - H_Sel*y;

    % Calculation of the integral error
    deM_intdt = cMset-cM;

    % Calculation of the controller action
    Ceq_uM = u - max(z0(31:33) + KM.p(:,H_Ver).*deM_intdt + KM.I(:,H_Ver).*eMint,0.0001);

    % Rearrangement of dTdt and Ceq
    ddt = vertcat(ddt{:}, deM_intdt);
    Ceq = vertcat(Ceq{:}, Ceq_uM, Ceq_cM);
else
    % Assignment of the differential and algebraic states
    x = T;
    p = [u; d_c2];
end

% Definition of the number of states, algebraic variables and inputs
n.x = length(x);
n.z = length(z);
n.y = length(y);
n.d = length(d_c1);
n.u = length(u);
n.p = n.u + n.d;

% Definition of boundary conditions
lb.x = 250*ones(n.x,1);   ub.x = 800*ones(n.x,1);
lb.z = zeros(n.z,1);      ub.z = [ones(2*N,1); 800*ones(n.u,1); ones(2*(n.u-1),1); kron([1; 1],[100; 250; 800; 1; 1]); 1e6];
lb.p = [0;0;0;5e-3;-.2];  ub.p = [ones(n.u,1); .2];

%  Defnition of  the DAE system
dae = struct('x',x,'z',z,'p',p,'ode',ddt,'alg',Ceq,'quad',L);

% Definition of the integrator
opts = struct('tf',t);
func = integrator('func','idas',dae,opts);