function [xdot] = dynModel(t,x,par)

%parametres
F =par.F;
CAin = par.CAin; 
CBin = par.CBin; 
Hrx = par.Hrx;
rho = par.rho;
tau = par.tau;
E1 = par.E1; 
E2 = par.E2;
Cp = par.Cp;
C1 = par.C1;
C2 = par.C2;

R  = par.R;
% States
CA = x(1);
CB = x(2);
T  = x(3);
eint = x(4);
% Input (=setpoint CS)

csSet = par.u.csSet; %actually this is the setpoint for cs now
H       = [-0.7689    0.6394    0.0046];
xnom    = [0.4978    0.5022      426.6998]';
y       = x(1:3)-xnom;
cs      = H*y;
Kc      = par.Kc;
taui	= par.taui;
e       = (cs-csSet);
Ti0 = par.u.Ti0; %initial (nominal input) Process input
Ti =Ti0-Kc*e-(Kc/taui)*eint; % Process input





if t >=par.dist3Time
	E2=par.dist3factor*E2;
end


[CAin CBin]=scenario(t,par);
% kinetics:
k1 = C1*60*exp(-E1/R/T);
k2 = C2*60*exp(-E2/R/T);
r = (k1*CA-k2*CB);

CAdot = 1/tau*(CAin-CA)-r;
CBdot = 1/tau*(CBin-CB)+r;
Tdot  = 1/tau*(Ti-T)+(-Hrx)/(rho*Cp)*r;

xdot = [CAdot;CBdot;Tdot;e];
