function Tprime=HEX(t,X,U)

% The function is called by the s-function HEXs.m and calculates the 
% derivatives of the temperatures in the hot side, cold side and wall 
% elements, plus the hot side outlet temperature after bypass
% in a countercurrent gas-gas heat exchanger. 

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

global par

%  Collecting parameters

%%
Nhex    =   par.Nhex;
hh      =   par.hh;
hc      =   par.hc;
cph     =   par.cph;
cpc     =   par.cpc;
cpw     =   par.cpw;
mw      =   par.mw;
Vh      =   par.Vh;
Vc      =   par.Vc;
Ah      =   par.Ah;
Ac      =   par.Ac;
Aw      =   par.Aw;
P       =   par.P;
R       =   par.R;
M       =   par.M;
Hvap    =   par.Hvap;



%%

% Splitting the states
Th      =   X(1:Nhex);                          % Gas temperatures in elements 1 to Nhex in hot side                         
Tc      =   X(Nhex+1:2*Nhex);                   % Gas temperatures in elements 1 to Nhex in cold side    
Tw      =   X(2*Nhex+1:3*Nhex);                 % Wall temperature in elements 1 to Nhex        
Thout   =   X(3*Nhex+1);                        % Hot outlet temperature after mixing bypass and heat exchanged stream

% Inputs and disturbances
Tc0         = U(1);                             % Inlet temperature cold side  [C]
Th0         = U(2);                             % Inlet temperature hot side (outlet from reactor) [C]
m           = U(3);                             % Mass flow of gas [kg/s]
mair        = U(4);                             % Mass flow of air fed to process [kg/s]
a           = U(5);                             % Hot side bypass fraction (a=0 -> no bypass)

% Defining length of output vectors
dThdt       =   zeros(Nhex,1);
dTcdt       =   zeros(Nhex,1);
dTwdt       =   zeros(Nhex,1);
dThoutdt    =   zeros(1,1);

if Th0 < 100 
    wH2O    =   0;
elseif a == 1
    wH2O    =   0;
else
    wH2O    =   par.wH2O;
end 


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


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

    dTcdt(1) =     (Tc0-Tc(1))*(((m+mair)*R*(Tc(1)+273.15))/(P*M*Vc)) ...
                    - (Tc(1)-Tw(Nhex))*(hc*Ac*R*(Tc(1)+273.15)/(cpc*P*M*Vc));
    
    dThdt(1) =     (Th0-Th(1))*(((m+mair)*(1-a)*R*(Th(1)+273.15))/(P*M*Vh))...
                    - (Th(1)-Tw(1))*(hh*Ah*R*(Th(1)+273.15)/(cph*P*M*Vh));


    dTwdt(1) =     ((Th(1)-Tw(1))*Ah*hh  +  (Tc(Nhex)-Tw(1))*Ac*hc)/(mw*cpw);

   

% ----------------------- Remaining elements ------------------------------

for i=2:Nhex
    
    
    
    dTcdt(i) =     (Tc(i-1)-Tc(i))*(((m+mair)*R*(Tc(i)+273.15))/(P*M*Vc))...
                   -  (Tc(i)-Tw(Nhex+1-i))*(hc*Ac*R*(Tc(i)+273.15)/(cpc*P*M*Vc));
   

    dTwdt(i) =      ((Th(i)-Tw(i))*Ah*hh  +  (Tc(Nhex+1-i)-Tw(i))*Ac*hc)/(mw*cpw);

end

for i=2:Nhex-1
    
        dThdt(i) =     (Th(i-1)-Th(i))*(((m+mair)*(1-a)*R*(Th(i)+273.15))/(P*M*Vh))...
                    - (Th(i)-Tw(i))*(hh*Ah*R*(Th(i)+273.15)/(cph*P*M*Vh));
end

dThdt(Nhex) =   (Th(Nhex-1)-Th(Nhex))*(((m+mair)*(1-a)*R*(Th(Nhex)+273.15))/(P*M*Vh))...
                    - (Th(Nhex)-Tw(Nhex))*(hh*Ah*R*(Th(Nhex)+273.15)/(cph*P*M*Vh)) - (wH2O*m*Hvap*R*(Th(Nhex)+273.15)/(cph*P*M*Vc));
% - Hot outlet temperature after mixing bypass and heat exchanged stream - 

    dThoutdt = ((m+mair)*a*Th0  + (m+mair)*(1-a)*Th(Nhex) - (m+mair)*Thout)*((R*(Thout+273.15))/(P*M*Vh)); 
    

% Defines output vector 

    Tprime=[dThdt; dTcdt; dTwdt; dThoutdt];

end
