function [qu] = qSensitivities(par,x_val,u_val,d_val)

% Written by Krishnamoorthy, Oct 2018.
% email - dineshk@ntnu.no

import casadi.*

nA = MX.sym('nA',1);
nB = MX.sym('nB',1);
nC = MX.sym('nC',1);

FA = MX.sym('FA',1);
FB = MX.sym('FB',1);

k1 = MX.sym('k1',1);

CA = nA/(par.V);
CB = nB/(par.V);
CC = nC/(par.V);
qr = (par.H1*k1*CA*CB*(par.V) + par.H2*par.k2*CB^2*(par.V))./1000000;

dx1 = (FA*par.CAin - (FA+FB)*CA - k1*CA*CB*par.V);
dx2 = (FB*par.CBin - (FA+FB)*CB - k1*CA*CB*par.V - 2*par.k2*CB^2*par.V);
dx3 = (-(FA+FB)*CC + k1*CA*CB*par.V);

diff = vertcat(dx1,dx2,dx3);

x_var = vertcat(nA,nB,nC);
p_var = vertcat(FA,FB);
d_var = k1;

g = {};
g = {g{:},vertcat(diff)};

rf_function = Function('rf_function',{x_var,p_var,d_var},...
                    {vertcat(g{:}),qr},{'x','u','d'},{'g','q'});
rf  = rootfinder('rf','newton',rf_function);

%% Get sensitivity matrices

% qu = rf.jacobian('u','q'); % Used in older version of CasADi 3.1.0
qu = rf.factory('qu',{'x','u','d'},{'jac:q:u'});

qu = full(qu(x_val,u_val,d_val));
