function dTdt=incin(t,X,U)

% This function is called by the s-function incins.m and calculates the
% derivative of the temperature of the gas leaving the electric heater.

%            X(1:N)    - Gas temperature in element 1 to N [C]
%            X(N+1:2*N)- Catalyst temperature in element 1 to N [C]

% The parameters are given in "parameters.m".

global par 

% Collecting parameters

%%
N           =   par.N;
cpcat       =   par.cpcat ;
cpgas       =   par.cpgas;
UAcat       =   par.UAcat;
%UAcatinf    =   par.UAcatinf;
UAwall      =   par.UAwall;
kAcat       =   par.kAcat;
Tamb        =   par.Tamb;
mcat        =   par.mcat;
deltaHc     =   par.deltaHc;
Vg          =   par.Vg;
P       =   par.P;
R       =   par.R;
M       =   par.M;
%MwMIBK  = par.MwMIBK;

%%


% Defining the states
T = X(1:2*N)';                              % X(1:N) - Gas temperature in element 1-N [C] 
                                            % X(N+1:2*N)- Catalyst temperature in element 1-N [C]


% Inputs and disturbances
T0      = U(1);                             % Inlet temperature to incinerator [C]
m       = U(2);                             % Mass flow of process gas [kg/s]
mair    = U(3);                             % Mass flow of air fed to process [kg/s]
w       = U(4);                             % Mass fraction of VOC

beta=1;




%%%%%%%%%%%%%%   Calculation of state derivatives   %%%%%%%%%%%%%%%%%%%%%%

% -------------------------  Initial elements -----------------------------



% Gas phase
dTdt(1)=    (((N*(m+mair)*R*(T(1)+273.15))/(P*M*Vg))*(T0-T(1)))...           
            -((UAcat*R*(T(1)+273.15)/(P*M*Vg*cpgas))*(T(1)-(T(N+1))));

% Catalyst phase        
dTdt(N+1)=  (UAcat/(mcat*cpcat)*(T(1)-T(N+1)))...
            +((kAcat*N/(mcat*cpcat))*(T(N+2)-T(N+1))) -(UAwall/(mcat*cpcat))*(T(N+1)-(Tamb));

% -------------------------- Interior elements ----------------------------
% Gas phase
for i=2
        
                        %UAcat=par.UAcatinf;
                        
            dTdt(i)=    (((N*(m+mair)*R*(T(i)+273.15))/(P*M*Vg))*(T(i-1)-T(i)))...
                        -((UAcat*R*(T(i)+273.15)/(P*M*Vg*cpgas))*(T(i)-(T(N+i))))+(1-beta)*(m*w*deltaHc*N*R*(T(i)+273.15)/(P*M*Vg*cpgas));
% Catalyst phase   



            dTdt(N+i)=  (UAcat/(mcat*cpcat)*(T(i)-T(N+i))) + (kAcat*N*N/(mcat*cpcat))*(T(N+i-1)+T(i+N+1)-2*T(N+i))+beta*(m*w*N*deltaHc/(mcat*cpcat)) ;
                        w=0; % Sets mass fraction of VOC to 0 after combustion
                        %UAcat=par.UAcat;
  for i=3:(N-1)                  
                   % Below ignition temperature
             dTdt(i)=    (((N*(m+mair)*R*(T(i)+273.15))/(P*M*Vg))*(T(i-1)-T(i)))...
                         -((UAcat*R*(T(i)+273.15)/(P*M*Vg*cpgas))*(T(i)-(T(N+i))));
            dTdt(N+i)=  (UAcat/(mcat*cpcat)*(T(i)-T(N+i))) + (kAcat*N*N/(mcat*cpcat))*(T(N+i-1)+T(i+N+1)-2*T(N+i))...
                        -(UAwall/(mcat*cpcat))*(T(N+i)-(Tamb));
        end



                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        end
        




%----------------------------  Final elements -----------------------------
% Gas phase
    dTdt(N)=    (((N*(m+mair)*R*(T(N)+273.15))/(P*M*Vg))*(T(N-1)-T(N)))...
                -((UAcat*R*(T(N)+273.15)/(P*M*Vg*cpgas))*(T(N)-T(2*N)));     
%Catalyst phase
    dTdt(2*N)=  (UAcat/(mcat*cpcat)*(T(N)-T(2*N)))    +((kAcat*N/(mcat*cpcat))*(T((2*N)-1)-T(2*N)))...
                -(UAwall/(mcat*cpcat))*(T(2*N)-(Tamb));   % Catalyst


 % Defining output vector
 
 dTdt = [dTdt];
 