function [func,dae,d,lb,ub,n,Cineq,ToutFun] = AmmonReac_NewDist(t,KS,KM,H,H_Ver)

% This function defines the differential-algebraic system of equation DAE
% for the ammonia reactor defined in Morud 1998. 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.
%
% The input is given as
%       t:          Time step for the integrator
%       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:9)
%                   .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 == 4
    H_Ver = 1;
    H_Sel = H;
elseif nargin == 5
    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.118e6 2.7e6
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 feed variables
d_nom = zeros(4,1);
d_nom(1) = 252/3.6;         % Feed flow (kg/s)
d_nom(2) = 200;             % Pressure [bar]
d_nom(3) = 250;             % Temperature [C]
d_nom(4) = 0.08;            % Mass fraction ammonia [-]
d_nom(5) = 3;               % Hydrogen/nitrogen molar ratio [-]

% 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 [-]
R = MX.sym('R',N);          % H2/N2 Ratio in the system in each of the discretization units [-]
Tin = MX.sym('Tin',3);      % Temperature in the inlet flow of the beds [C]
cin = MX.sym('cin',3);      % Mass fraction of ammonia in the inlet flow of the beds [-]
Rin = MX.sym('Rin',3);      % H2/N2 Ratio in the system in the inlet flow of the beds [-]
TexHEx = MX.sym('TexHEx',1);    % Outlet temperature of the Heat Exchanger

d = MX.sym('d',6);          % Disturbances, inlet values (m, p, T, c, H2/N2) and catalyst activity
% Vd = MX.sym('Vd',3);        % Volume disturbance
in = d(1:5)+d_nom;          % Inlet values (m, p, T, c, H2/N2)

u = MX.sym('u',3);          % Split ratio to Reac1, Reac2, and Reac3

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

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

%% Setup of the model
% Redefinition of massflows and calculation of the molar flows
m1 = (1-u(2)-u(3))*in(1);
m2 = (1-u(3))*in(1);
m3 = in(1);
mSHEx = (1-u(1)-u(2)-u(3))*in(1);
[ndot.Fin] = MassFlow(in(4),in(5),in);

% Predefinition of the cells
Ceq_wNH3 = cell(N,1);
Ceq_R = cell(N,1);
dTdt = cell(N,1);
Ceqin = cell(10,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)*in(3));
Tout     = (1-eps*CStar)*T(N)+CStar*eps*in(3);

% Calculation of the reaction rates and flowrates
[r,ndot.sys] = Rrate_NewDist(T,wNH3,R,in,kval,d(6));

%% Reactor 1

% Calculation of the inlet variables
Ceqin{1} = m1*Tin(1) - (mSHEx*TexHEx + u(1)*in(1)*in(3));   % Temperature [K]
Ceqin{5} = cin(1) - in(4);                                  % Mass fraction ammonia [-]
Ceqin{8} = Rin(1) - in(5);                                  % H2/N2 ratio [-]

% Bed 1
dTdt{1}     = (m1*par.Cp*(Tin(1)-T(1)) - r(3,1)*dm1*par.dHrx)/(dm1*par.Cpc);
Ceq_wNH3{1} = m1*(cin(1)-wNH3(1)) - r(3,1)*dm1;
Ceq_R{1}    = R(1) - ((1-u(2)+u(3))*ndot.Fin(1) - r(1,1)*dm1/2)/...
                  ((1-u(2)+u(3))*ndot.Fin(2) - r(2,1)*dm1/28);

% 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_wNH3{j} = m1*(wNH3(j-1)-wNH3(j)) - r(3,j)*dm1;
    Ceq_R{j}    = R(j) - (ndot.sys(1,j-1) - r(1,j)*dm1/2)/...
                         (ndot.sys(2,j-1) - r(2,j)*dm1/28);
end

%% Reactor 2

% Calculation of the inlet variables
Ceqin{2} = m2*Tin(2) - (m1*T(N1) + (u(2)*in(1))*in(3));         % Temperature [K]
Ceqin{6} = m2*cin(2) - (m1*wNH3(N1) + (u(2)*in(1))*in(4));      % Mass fraction ammonia [-]
Ceqin{9} = (ndot.sys(1,N)+u(2)*ndot.Fin(1)) - ...
           (ndot.sys(2,N)+u(2)*ndot.Fin(2)) * ...
            Rin(2);                                             % 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_wNH3{N1+1} = m2*(cin(2)-wNH3(N1+1)) - r(3,N1+1)*dm2;
Ceq_R{N1+1}     = R(N1+1) - (ndot.sys(1,N1)+u(2)*ndot.Fin(1) - r(1,N1+1)*dm2/2)/...
                            (ndot.sys(2,N1)+u(2)*ndot.Fin(2) - r(2,N1+1)*dm2/28);

% 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_wNH3{j} = m2*(wNH3(j-1)-wNH3(j)) - r(3,j)*dm2;
    Ceq_R{j}    = R(j) - (ndot.sys(1,j-1) - r(1,j)*dm2/2)/...
                         (ndot.sys(2,j-1) - r(2,j)*dm2/28);
end

%% Reactor 3

% Calculation of the inlet variables
Ceqin{3}  = m3*Tin(3) - (m2*T(N1+N2) + u(3)*in(1)*in(3));       % Temperature [K]
Ceqin{7}  = m3*cin(3) - (m2*wNH3(N1+N2) + u(3)*in(1)*in(4));    % Mass fraction  ammonia[-]
Ceqin{10} = (ndot.sys(1,N)+u(3)*ndot.Fin(1)) - ...
            (ndot.sys(2,N)+u(3)*ndot.Fin(2)) * ...
             Rin(3);                                            % H2/N2 ratio [-]

% 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_wNH3{N1+N2+1} = m3*(cin(3)-wNH3(N1+N2+1)) - r(3,N1+N2+1)*dm3;
Ceq_R{N1+N2+1}    = R(N1+N2+1) - (ndot.sys(1,N1+N2)+u(3)*ndot.Fin(1) - r(1,N1+N2+1)*dm2/2)/...
                              (ndot.sys(2,N1+N2)+u(3)*ndot.Fin(2) - r(2,N1+N2+1)*dm2/28);

% 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_wNH3{j} = m3*(wNH3(j-1)-wNH3(j)) - r(3,j)*dm3;
    Ceq_R{j}    = R(j) - (ndot.sys(1,j-1) - r(1,j)*dm3/2)/...
                         (ndot.sys(2,j-1) - r(2,j)*dm3/28);
end

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

%% Combination of the variables of the system
    
% Preassignment of the algebraic states
z = [wNH3; R; Tin; TexHEx; cin; Rin];

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

% Rearrangement of differential and algebraic equations
Ceq = vertcat(Ceq_wNH3{:}, Ceq_R{:}, Ceqin{:});
ddt = vertcat(dTdt{:});
        
% 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];
    
    % 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];
                
        % 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(2*N+1:2*N+3) + 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];
    
    % 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(p0(1:3) + 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];
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);
n.u = length(u);
n.p = n.u + n.d;

% Definition of boundary conditions
lb.x = 150*ones(n.x,1);   ub.x = 800*ones(n.x,1);
lb.z = [zeros(N,1); 1.5*ones(N,1); 150*ones(n.u+1,1); zeros(n.u,1); 1.5*ones(n.u,1)];
ub.z = [ones(N,1);    5*ones(N,1); 800*ones(n.u+1,1); ones(n.u,1);    5*ones(n.u,1)];
lb.p = zeros(n.p,1);      ub.p = ones(n.p,1);
lb.u = 1e-4*ones(n.u,1);      ub.u = ones(n.u,1);

% Definition of the DAE system
L = -(wNH3(end)-cin(1))*in(1);
dae = struct('x',x,'z',z,'p',p,'ode',ddt,'alg',Ceq,'quad',L);
ToutFun = Function('ToutFun',{[x;z;p]},{Tout});

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