function F_integrator = CSTR_plant(par)
import casadi.*


cp = 1000;
R = 1.987;
Hr = 5000;
rho = 1;
tau = 1*60;
pCB = 2.009;
pTi = 1.657e-3;

dist = [1;0.0];

% x
CA = MX.sym('CA',1);
CB = MX.sym('CB',1);
T = MX.sym('T',1);

% dist
CAin = MX.sym('CAin',1);
CBin = MX.sym('CBin',1);

% u
Tin = MX.sym('Tin',1);

% p
A1 = MX.sym('A1',1);
A2 = MX.sym('A2',1);

E1 = MX.sym('E1',1);
E2 = MX.sym('E2',1);

du = 1;
% ODE system
dx1 = (CAin - CA)/tau - (A1*exp(-E1/(R*T))*CA - A2*exp(-E2/(R*T))*CB);
dx2 = (CBin - CB)/tau + (A1*exp(-E1/(R*T))*CA - A2*exp(-E2/(R*T))*CB);
dx3 = (Tin - T)/tau + Hr*((A1*exp(-E1/(R*T)))*CA - (A2*exp(-E2/(R*T)))*CB)/(rho*cp);

diff = vertcat(dx1,dx2,dx3);

diff = substitute(diff,A1,par.A1);
diff = substitute(diff,A2,par.A2);
diff = substitute(diff,E1,par.E1);
diff = substitute(diff,E2,par.E2_plant);

x_var = vertcat(CA,CB,T);
p_var = vertcat(Tin,CAin,CBin);
d_var = vertcat(CAin,CBin);
L = -(pCB*CB - (pTi*Tin).^2);

ode = struct('x',x_var,'p',p_var,'ode',diff,'quad',L);  % der(m_tot) = w_in - w_out;
opts = struct('tf',par.Ts);


% create IDAS integrator
F_integrator = integrator('F','cvodes',ode,opts);

