function [func,ode,d,n,par] = AmmonReac_ODE_Ver1(t,KS,KM)

% This function defines the differential system of equation ODE 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
%
%       ode:        Structure containing the dae with the following fields
%                   .x      Differential variables (Temperatures, 
%                           Concentrations)
%                   .p      Parameters (1:3, Split ratios to reactor beds)
%                           and disturbances (4:7)
%                   .ode    Differential equations
%                   .quad   Cost function (Ammonia extent of reaction [kg/s]
%
%       d:          Variables for the input disturbance, necessary for
%                   substitution
%
%       n:          Structure containing the number of variables with the
%                   following fields
%                   .x      Differential variables
%                   .z     Auxilary vector
%                   .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);

% Import of CASADI
addpath([dirmain '\Util1. Casadi3.1.0'])
addpath([dirmain '\Util2. Functions'])
addpath(genpath([dirmain '\Util4. MATLAB2tikz']))
import casadi.*

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

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

% Import of CASADI
import casadi.*

%% 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 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]
par.rhoG = 50;              % Typical Density of the Gas [kg/m3]
par.rhoC = 2200;            % Typical Density of the Catalyst [kg/m3]
par.alph = 0.33;            % Bed Void Fraction [-]

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

% 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(6,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 [-]
d_nom(6) = 1;               % Catalyst activity 

% Definition of the system states
T = MX.sym('T',N);          % Temperature in each of the discretization units
wNH3 = MX.sym('wNH3',N);    % Mass fraction of ammonia in each of the discretization units [-]
d = MX.sym('d',6);          % Disturbances, inlet values (m, p, T, wNH3, H2/N2) and catalyst activity
in = d(1:6)+d_nom;          % Inlet values (m, p, T, wNH3, H2/N2)

eSint = MX.sym('eSint',3);  % Integrated error of the slave
deM_intdt = MX.sym('deM_intdt',3);  % Integrated error of the master
eMint = MX.sym('eMint',3);  % Integrated error of the master

if exist('KS','var')
    % Calculation of the controller action
    u = max(p0(1:3) + KS.I.*eSint,0.0001);
 
    % SOC Master controller is used
    if exist('KM','var')
        % Calculation of the controller action
        cSset = (z0(31:33) + KM.p .*deM_intdt + KM.I .*eMint);
    else
        cSset = MX.sym('TSset',3);  % Setpoint of the inlet temperatures of each bed
    end
else
    u = MX.sym('u',3);          % Split ratio to Reac1, Reac2, and Reac3
end

%% Setup of the model
% Redefinition of massflows
mS1 = u(1)*in(1);
mS2 = u(2)*in(1);
mS3 = u(3)*in(1);
mSHEx = (1-u(1)-u(2)-u(3))*in(1);

% Predefinition of the cells
dTdt    = cell(N,1);
dwNH3dt = cell(N,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)));
TexHEx = (eps*T(N)+(1-eps)*in(3));

% Reactor 1
m1 = mS1+mSHEx;
Tin1 = (mSHEx*TexHEx + mS1*in(3))/m1;
wNH3in1 = in(4);

    % Bed 1
    [r] = Rrate_NewDist_MESC(T(1),wNH3(1),in,kval,d(6));
    dTdt{1}    = (m1*par.Cp*(Tin1-T(1)) + r*dm1*par.dHrx)/(dm1*par.Cpc);
    dwNH3dt{1} = (m1*(wNH3in1-wNH3(1)) + r*dm1)/(par.alph*par.rhoG*V1/N1);

    % Bed 2:N1x
    for j = 2:N1
        [r]        = Rrate_NewDist_MESC(T(j),wNH3(j),in,kval,d(6));
        dTdt{j}    = (m1*par.Cp*(T(j-1)-T(j)) + r*dm1*par.dHrx)/(dm1*par.Cpc);
        dwNH3dt{j} = (m1*(wNH3(j-1)-wNH3(j)) + r*dm1)/(par.alph*par.rhoG*V1/N1);
    end

% Reactor 2
m2 = mS2+m1;
Tin2 = (m1*T(N1) + mS2*in(3))/m2;
wNH3in2 = (m1*wNH3(N1) + mS2*in(4))/m2;

    % Bed 1
    [r] = Rrate_NewDist_MESC(T(N1+1),wNH3(N1+1),in,kval,d(6));
    dTdt{N1+1}    = (m2*par.Cp*(Tin2-T(N1+1)) + r*dm2*par.dHrx)/(dm2*par.Cpc);
    dwNH3dt{N1+1} = (m2*(wNH3in2-wNH3(N1+1)) + r*dm2)/(par.alph*par.rhoG*V2/N2);

    % Bed 2:N1
    for j = N1+2:N1+N2
        [r] = Rrate_NewDist_MESC(T(j),wNH3(j),in,kval,d(6));
        dTdt{j}    = (m2*par.Cp*(T(j-1)-T(j)) + r*dm2*par.dHrx)/(dm2*par.Cpc);
        dwNH3dt{j} = (m2*(wNH3(j-1)-wNH3(j)) + r*dm2)/(par.alph*par.rhoG*V2/N2);
    end

% Reactor 3
m3 = in(1);
Tin3 = (m2*T(N1+N2) + mS3*in(3))/m3;
wNH3in3 = (m2*wNH3(N1+N2) + mS3*in(4))/m3;

    % Bed 1
    [r] = Rrate_NewDist_MESC(T(N1+N2+1),wNH3(N1+N2+1),in,kval,d(6));
    dTdt{N1+N2+1}    = (m3*par.Cp*(Tin3-T(N1+N2+1)) + r*dm3*par.dHrx)/(dm3*par.Cpc);
    dwNH3dt{N1+N2+1} = (m3*(wNH3in3-wNH3(N1+N2+1)) + r*dm3)/(par.alph*par.rhoG*V3/N3);
    
    % Bed 2:N1
    for j = N1+N2+2:N
        [r] = Rrate_NewDist_MESC(T(j),wNH3(j),in,kval,d(6));
        dTdt{j}	   = (m3*par.Cp*(T(j-1)-T(j)) + r*dm3*par.dHrx)/(dm3*par.Cpc);
        dwNH3dt{j} = (m3*(wNH3(j-1)-wNH3(j)) + r*dm3)/(par.alph*par.rhoG*V3/N3);
    end

Tin    = [Tin1;    Tin2;    Tin3];
wNH3in = [wNH3in1; wNH3in2; wNH3in3];

%% Combination of the variables of the system

% Preassignment of the algebraic states
z = [wNH3; 3*ones(N,1); Tin; TexHEx; wNH3in; 3*ones(3,1)];

% 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)];

% Temperature controller is used (with ot without a master controller)
if exist('KS','var')
    
    % Assignment of the differential and algebraic states
    x = [T; eSint; wNH3];
    z = [z; u];
    p = [cSset; d];
        
    % Calculation of the integral error
    deS_intdt = cSset-Tin;
    
    % Rearrangement of dTdt
    ddt = vertcat(dTdt{:},deS_intdt,dwNH3dt{:});
    
    % SOC Master controller is used
    if exist('KM','var')
        % Assignment of the differential and algebraic states
        x = [T; eSint; eMint; wNH3];
        z = [z; cSset; d];
        p = [cMset];
        
        % Calculation of the integral error
        deM_intdt = 0-H_Sel*y;
        
        % Rearrangement of dTdt
        ddt = vertcat(dTdt{:}, deS_intdt, deM_intdt, dwNH3dt{:});
    end
    
else
       % Assignment of the differential and algebraic states
    x = [T; wNH3];
    p = [u; d];
            
    % Rearrangement of dTdt
    ddt = vertcat(dTdt{:}, dwNH3dt{:});
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 the DAE system
L = -(wNH3(end)-wNH3in1)*in(1);
ode = struct('x',x,'p',p,'ode',ddt,'quad',L);

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

