%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function to simulate the CSTR using a 
% pi controller and a self-optimizing variable 
% combination
% Johannes Jäschke
% 2010

function [xdot] = dynModelSOC(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;
tsample = par.tsample;
tsamplep= par.tsamplep;

R  = par.R;
% States
CA = x(1);
CB = x(2);
T  = x(3);
eint = x(4); % added state for the PI controller
% Nominal Input 
Ti0 = par.u.Ti0;

%self-Optimizing Variable Combination:
H	= [-0.7689    0.6394    0.0046]; % from conventional method

xnom 	=[0.4978    0.5022      426.6998]';

y	= x(1:3)-xnom;
cs	= H*y;
Kc   	= par.Kc;
taui 	= par.taui;

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

csSet	= 0;
e	= cs-csSet;
Ti =Ti0-Kc*e-(Kc/taui)*eint;

% Add noise
k=floor(t/(tsample+tsamplep))+1; % find noise value


% load scenario
[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];
