function [sys,x0,str,ts] = kurihara_model(t,x,u,flag)
% fcc_model is an S-function which implement a simplified model of FCC unit
% (simplified Kurihara model)
% Réf:
% [1] Mc Farlane & Bacon (Queen's University)
% [2] Lee & Kugelman (MOBIL Research)

%
% Dispatch the flag.
%
switch flag
  case 0
    [sys,x0,str,ts]=mdlInitializeSizes; % Initialization
  case 1
    sys = mdlDerivatives(t,x,u); % Calculate derivatives
  case 3
    sys = mdlOutputs(t,x,u); % Calculate outputs  
  case { 2, 4, 9 } % Unused flags
    sys = [];
  otherwise
    error(['Unhandled flag = ',num2str(flag)]); % Error handling
end
% End of csfunc.
%==============================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the 
% S-function.
%==============================================================
%
function [sys,x0,str,ts] = mdlInitializeSizes()
global  Tra Trg  Crc Csc Ccat

fcc_init;
%
% Call simsizes for a sizes structure, fill it in and convert it 
% to a sizes array.
%
sizes = simsizes;
sizes.NumContStates  = 5;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 5;
sizes.NumInputs      = 2;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);
%
% Initialize the initial conditions.
%
x0 = [Tra;Trg;Crc;Csc;Ccat];
%load('xfinal.mat');
%x0=x;
%
% str is an empty matrix.
%
str = [];
%
% Initialize the array of sample times; in this example the sample 
% time is continuous, so set ts to 0 and its offset to 0.
%
ts = [0 0];
% End of mdlInitializeSizes.
%
%==============================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%==============================================================
function sys = mdlDerivatives(~,x,u)
global R Hra Pra Tra DHcr Rtf Ttf Stf Acr kcr Acc kcc Fcf...
    Hrg Prg Trg DHfv DHrg Rai Tai Sai Rrc Sc Aor Crc Csc Ccat
Tra=x(1);
Trg=x(2);
Crc=x(3);
Csc=x(4);
Ccat=x(5);
Rai=u(1);
Rrc=u(2);

Kcc = kcc/(Ccat*Crc^0.06)*exp( -Acc/(R*Tra) );
Rcc = Kcc*Pra*Hra;
Rad = Fcf*Rtf;
Rcf = Rcc + Rad;
Rsc = Rrc + Rcf;

Kcr = kcr/(Ccat*Crc^0.15)*exp( -Acr/(R*Tra) );
Rcr = Kcr*Pra*Hra;
Ctf = 1 / ( 1 + Rtf/Rcr);
dTra = ( Sc*Rrc*(Trg-Tra) + Stf*Rtf*(Ttf-Tra) - Rtf*DHfv - Rtf*Ctf*DHcr ) / (Sc*Hra);

Kor = 1.16e-5*exp( Aor/R*( 1/866.7 - 1/Trg ) );
Kod = 6.34e-9*Rai^2;
Ofg = 21*exp( - (Prg*Hrg/Rai) / ( 1/Kod + 1/(Kor*Crc) ) );
Rcb = Rai*(21-Ofg)/242;
dTrg = ( Sc*Rsc*(Tra-Trg) + Sai*Rai*(Tai-Trg) + Rcb*DHrg ) / (Sc*Hrg);

dCrc = ( Rsc*Csc-Rrc*Crc - 100*Rcb )/Hrg;

dCsc = ( Rrc*Crc-Rsc*Csc + 100*Rcf )/Hra;

dCcat = ( -Rrc*Ccat + 100*Rcc ) /Hra;

sys = 3600*[dTra;dTrg;dCrc;dCsc;dCcat];
% End of mdlDerivatives.
%
%==============================================================
% mdlOutputs
% Return the block outputs.
%==============================================================
%
function sys = mdlOutputs(~,x,~)
%save('xfinal.mat','x')
sys=x;
% End of mdlOutputs.