function [CeqC,JacC,Ceq,Jac,dCeqdx,into,x,colind] = ConReact(a,par,scl,con,~,def,itg)

% This function defines the flowsheet of the reactor section of the
% simplified ammonia cycle. It can be used to define the surrogate model of
% this part of the plant. It consists of the recycle heat exchanger, the
% heat exchanger network for heating the reactor and the two reactors
% itself. It can be easily extended to more complicated models if needed.
%
% Written by Julian Straus, Trondheim 09.12.2018
% 
% The input of this function is given by
%       a:      States of the system, length(a) is depending on the
%               investigated model.
%
% The main output is defined as
%       CeqC:   Equality constraints of the non-linear system of equations.
%
%       JacC:   Jacobian of the system.
%
% The seconday output for post analysis is given by the following
% structures. Each structure includes values for each subunit.
%       Ceq:    Cell structure containing the equality contraints of the
%               UOs.
%
%       Jac:    Cell structure containing the full Jacobians of the UOs at
%               the right position.
%
%       dCeqdx: Cell structure containing the full jacobians of the UOs.
%
%       into:   Cell structure containing the Input vector to the UOs.
%
%       x:      Cell structure containing the reassigned states of the UOs.
%
%       colind: Column indices for the Jacobians of the UOs.
%
% The necessary structures are defined as
%       par:    Cell structure containing the parameters used in the 
%               calculation, output of the system.
%
%       scl:    Cell structure containing the scaling factors for each unit
%               operations. The assignment is given in par.ReactSec.
%
%       con:    Physical constants, must contain
%               .R      Ideal Gas constant
%
%       def:    Defined stream values, e.g pressure and composition for all
%               inlets and the pressure of the outlets or the flow of the
%               inlets and:
%               .sizeMat The size of the cell array par.ReactSec.
%
%       itg:    Intiger values, which should include
%               .k      Total number of variables per stream
%               .j      Total number of components in each stream

% Reassignment of the solution vector
x = cell(def.sizeMat(1)+1,1);
into = cell(def.sizeMat);
Ceq = into;
dCeqdx = into;
rowlen = into;
colind = into;
ctvar = into;
x{end} = a;

% Predefinition of the Jacobian of the System
Jac = into;%sparse(length(a),length(a));

% Definition of the positions of HEx2, Split1, HEx3, and HEx4
l1 = 1;
l2 = 3;
l3 = 7;
l4 = 9;

%% Heat Exchanger 2
% Reassignment of the variables for heat exchanger 2
k = 1;
x{k} = a(1:32);
a = a(33:end);
ctvar{k} = itg.R;

% Definition of the input values for heat exchanger 2
into{k}= [x{k}(1:8);...         Hot inlet flow, pressure, and temperature 
          x{k}(9:16);...        Hot outlet flow, pressure, and temperature 
          x{k}(17:24);...       Cold inlet flow, pressure, and temperature 
          x{k}(25:32)];       % Cold outlet flow, pressure, and temperature 

% Calculation of the equality constraints given by heat exchanger 2
[Ceq{k},dCeqdx{k}] = HExeNTU(into{k},par{k},scl{k},con,[],itg);

% Calculation of the parts of the Jacobian defined by heat exchanger 2
rowlen{k} = length(Ceq{k});                         % Number of the rows
colind{k} = ctvar{k}+1:ctvar{k}+length(x{k});       % Indexes of the all inlet and outlet columns
Jac{k} = sparse(rowlen{k},itg.nvar);
Jac{k}(:,colind{k}) = dCeqdx{k}(:,1:32);

%% Split 1
% Reassignment of the variables for split 1
x{k+2} = a(1:26);
a = a(27:end);
ctvar{k+2} = ctvar{k}+numel(x{k});
k = k+2;

% Definition of the input values for split 1
into{k} = [into{k-2}(25:32);...	Inlet flow, pressure, and temperature
           x{k}(1:8);...        Outlet 1 flow (Valve), pressure, and temperature
           x{k}(9:16);...       Outlet 2 flow (Post Reactor), pressure, and temperature
           x{k}(17:24);...      Outlet 3 flow (Interstage), pressure, and temperature
           x{k}(25:26); 0];...  Split Ratios
             
% Calculation of the equality constraints given by split 1
[Ceq{k},dCeqdx{k}] = Split(into{k},par{k},scl{k},[],[],itg);

% Calculation of the parts of the Jacobian defined by split 1
rowlen{k} = length(Ceq{k});                         % Number of the rows
colind{k} = [colind{k-2}(25:32) ...                   Indexes of the inlet columns
             ctvar{k}+1:ctvar{k}+length(x{k})];     % Indexes of the outlet columns
Jac{k} = sparse(rowlen{k},itg.nvar);
Jac{k}(:,colind{k}) = dCeqdx{k}(:,1:34);

%% Valve 1 (Outlet 1 of the Split)
% Reassignment of the variables for the Valve 1
x{k+2} = a(1:9);
a = a(10:end);
ctvar{k+2} = ctvar{k}+numel(x{k});
k = k+2;

% Definition of the input values for the first valve
into{k} = [into{l2}(9:16);...	Inlet flow, pressure, and temperature
           x{k}(1:8)];        % Outlet flow, pressure, and temperature

par{k}.K = x{k}(9)*scl{k}.Var(end); % Manipulated variable for the flowrate

% Calculation of the equality conditions given by Valve 1
[Ceq{k},dCeqdx{k}] = ValveSimp(into{k},par{k},scl{k},[],[],itg);

% Calculation of the parts of the Jacobian defined by valve 1
rowlen{k} = length(Ceq{k});                         % Number of the rows
colind{k} = [colind{l2}(9:16) ...                     Indexes of the inlet columns
             ctvar{k}+1:ctvar{k}+length(x{k})];     % Indexes of the outlet columns
Jac{k} = sparse(rowlen{k},itg.nvar);
Jac{k}(:,colind{k}) = dCeqdx{k};
        
%% Reactor Heat Exchanger 3 Post Reactor 2 (Outlet 2 of the Split)
% Reassignment of the variables for heat exchanger 3
x{k+2} = a(1:25);
a = a(26:end);
ctvar{k+2} = ctvar{k}+numel(x{k});
k = k+2;

% Definition of the input values for heat exchanger 3
into{k} = [x{k}(1:8);...       	Hot inlet flow, pressure, and temperature 
           x{k}(9:16);...       Hot outlet flow, pressure, and temperature 
           into{l2}(17:24);...  Cold inlet flow, pressure, and temperature 
           x{k}(17:24)];      % Cold outlet flow, pressure, and temperature 

par{k}.kCold = x{k}(end)*scl{k}.Var(end);

% Calculation of the equality constraints given by heat exchanger 3
[Ceq{k},dCeqdx{k}] = HExeNTU(into{k},par{k},scl{k},con,[],itg);

% Calculation of the parts of the Jacobian defined by heat exchanger 3
rowlen{k} = length(Ceq{k});                         % Number of the rows
colind{k} = [ctvar{k}+1:ctvar{k}+16 ...               Indexes of the hot columns
             colind{l2}(17:24) ...                    Indexes of the cold inlet columns
             ctvar{k}+17:ctvar{k}+25];              % Indexes of the cold outlet columns
Jac{k} = sparse(rowlen{k},itg.nvar);
Jac{k}(:,colind{k}) = dCeqdx{k}(:,[1:32 34]);

%% Interstage Heat exchanger 4 (Outlet 3 of the Split)
% Reassignment of the variables for heat exchanger 4
x{k+2} = a(1:24);
a = a(25:end);
ctvar{k+2} = ctvar{k}+numel(x{k});
k = k+2;

% Definition of the input values for heat exchanger 4
into{k} = [x{k}(1:8);...       	Hot inlet flow, pressure, and temperature 
           x{k}(9:16);...       Hot outlet flow, pressure, and temperature 
           into{l2}(25:32);...  Cold inlet flow, pressure, and temperature 
           x{k}(17:24)];      % Cold outlet flow, pressure, and temperature 
         
% Calculation of the equality constraints given by heat exchanger 4
[Ceq{k},dCeqdx{k}] = HExeNTU(into{k},par{k},scl{k},con,[],itg);

% Calculation of the parts of the Jacobian defined by heat exchanger 4
rowlen{k} = length(Ceq{k});                         % Number of the rows
colind{k} = [ctvar{k}+1:ctvar{k}+16 ...               Indexes of the hot columns
             colind{l2}(25:32) ...                    Indexes of the cold inlet columns
             ctvar{k}+17:ctvar{k}+24];             	% Indexes of the cold outlet columns
Jac{k} = sparse(rowlen{k},itg.nvar);
Jac{k}(:,colind{k}) = dCeqdx{k}(:,1:32);

%% Mixer 2
% Reassignment of the variables for mixer 2
x{k+2} = a(1:8);
a = a(9:end);
ctvar{k+2} = ctvar{k}+numel(x{k});
k = k+2;

% Definition of the input values for mixer 2
into{k} = [into{k-6}(9:16);...  Inlet 1 flow, pressure, temperature, and composition
           into{l3}(25:32);...  Inlet 2 flow, pressure, temperature, and composition
           into{l4}(25:32);...  Inlet 3 flow, pressure, temperature, and composition
           x{k}(1:8)];        % Outlet flow, pressure, temperature, and composition

% Calculation of the equality conditions given by mixer 2
[Ceq{k},dCeqdx{k}] = Mix(into{k},par{k},scl{k},[],[],itg);

% Calculation of the parts of the Jacobian defined by mixer 2
rowlen{k} = length(Ceq{k});                         % Number of the rows
colind{k} = [colind{k-6}(9:16) ...                    Indexes of the inlet 1 columns
             colind{l3}(25:32) ...                    Indexes of the inlet 2 columns
             colind{l4}(25:32) ...                    Indexes of the inlet 3 columns
             ctvar{k}+1:ctvar{k}+length(x{k})];    	% Indexes of the outlet columns
Jac{k} = sparse(rowlen{k},itg.nvar);
Jac{k}(:,colind{k}) = dCeqdx{k};

%% Reactor 1
% Reassignment of the variables for reactor 1
x{k+2} = a(1:(par{k+2}.n-1)*itg.k);
a = a((par{k+2}.n-1)*itg.k+1:end);
ctvar{k+2} = ctvar{k}+numel(x{k});
k = k+2;

% Definition of the input values for reactor 1
into{k} = [into{k-2}(25:32);...	Inlet reactor 1
           x{k};...             Outlet reactor 1-n-1
           into{l4}(1:8)];    % Outlet reactor n

% Calculation of the equality constraints given by reactor 1
[Ceq{k},dCeqdx{k}] = CSTR(into{k},par{k},scl{k},con,[],itg);

% Calculation of the parts of the Jacobian defined by reactor 1
rowlen{k} = length(Ceq{k});                         % Number of the rows
colind{k} = [colind{k-2}(25:32) ...                   Indexes of the inlet of reactor 1 columns
             ctvar{k}+1:ctvar{k}+length(x{k}) ...     Indexes of the outlet of reactor 1-n columns
             colind{l4}(1:8)];                      % Indexes of the outlet of reactor n columns
Jac{k} = sparse(rowlen{k},itg.nvar);
Jac{k}(:,colind{k}) = dCeqdx{k};

%% Reactor 2
% Reassignment of the variables for reactor 2
x{k+2} = a(1:(par{k+2}.n-1)*itg.k);
a = a((par{k+2}.n-1)*itg.k+1:end);
ctvar{k+2} = ctvar{k}+numel(x{k});
k = k+2;

% Definition of the input values for reactor 2
into{k} = [into{l4}(9:16);...	Inlet reactor 1
           x{k};...          	Outlet reactor 1-n-1
           into{l3}(1:8)];    % Outlet reactor n

% Calculation of the equality constraints given by reactor 2
[Ceq{k},dCeqdx{k}] = CSTR(into{k},par{k},scl{k},con,[],itg);

% Calculation of the parts of the Jacobian defined by reactor 2
rowlen{k} = length(Ceq{k});                         % Number of the rows
colind{k} = [colind{l4}(9:16) ...                     Indexes of the inlet of reactor 1 columns
             ctvar{k}+1:ctvar{k}+length(x{k}) ...     Indexes of the outlet of reactor 1-n columns
             colind{l3}(1:8)];                      % Indexes of the outlet of reactor n columns
Jac{k} = sparse(rowlen{k},itg.nvar);
Jac{k}(:,colind{k}) = dCeqdx{k};

%% Primary cooling before heat exchanger 2, HEx 5
% Reassignment of the variables for heat exchanger 5
x{k+2} = a(1:18);
a = a(19:end);
ctvar{k+2} = ctvar{k}+numel(x{k});
k = k+2;

% Definition of the input values for heat exchanger 5
into{k} = [into{l3}(9:16);...	Hot inlet flow, pressure, and temperature 
           x{k}(1);...          Reference temperature
           into{l1}(1:8);...    Hot outlet flow, pressure, and temperature
           x{k}(2:9);...        Cold inlet flow, pressure, and temperature 
           x{k}(10:17);...      Cold outlet flow, pressure, and temperature 
           x{k}(18)];...        Transferred heat

% Calculation of the equality constraints given by the first heat exchanger
[Ceq{k},dCeqdx{k}] = HExSimpH(into{k},par{k},scl{k},con,[],itg);

% Calculation of the parts of the Jacobian defined by heat exchanger 5
rowlen{k} = length(Ceq{k});...                      Number of the rows
colind{k} = [colind{l3}(9:16) ...                   Indexes of the inlet columns
             ctvar{k}+1 ...                         Transferred heat
             colind{l1}(1:8) ...                    Indexes of the inlet columns
             ctvar{k}+2:ctvar{k}+length(x{k})];...  Indexes of the outlet columns
Jac{k} = sparse(rowlen{k},itg.nvar);
Jac{k}(:,colind{k}) = dCeqdx{k}(:,1:34);

%% Call and reassignment of the definitions within the submodels
[Ceqdef,Jacdef] = def.fun(itg.nvar,into,colind);

Ceq(2:2:end) = Ceqdef(2:2:end);
Jac(2:2:end) = Jacdef(2:2:end);

%% Combination of the equality constraints
CeqC = cell2mat(Ceq);
JacC = cell2mat(Jac);