function [qu] = qSensitivities_GL(x_val,z_val,u_val)
import casadi.*

% plant model
par = GasLiftRiser_Param();
par.tSim = 1;

n_w = par.n_w;
R = par.R;
Mw = par.Mw;

%% Modelling

L_w = par.L_w;
H_w = par.H_w;
D_w = par.D_w;

L_bh = par.L_bh;
H_bh = par.H_bh;
D_bh = par.D_bh;

L_a = par.L_a;
H_a = par.H_a;
D_a = par.D_a;

L_r = par.L_r;
H_r = par.H_r;
D_r = par.D_r;

rho_o = par.rho_o;
C_iv = par.C_iv;
C_pc = par.C_pc;
C_pr = par.C_pr;
rho_ro = sum(rho_o)/2;

mu_oil = 1*0.001; % 1cP oil viscosity

A_w = pi.*(D_w/2).^2;
A_bh = pi.*(D_bh/2).^2;
V_a = L_a.*(pi.*(D_a/2).^2 - pi.*(D_w/2).^2);
A_r = pi.*(D_r/2).^2;

% differential states
m_ga = MX.sym('m_ga',n_w); % 1-2
m_gt = MX.sym('m_gt',n_w); % 3-4
m_ot = MX.sym('m_ot',n_w); % 5-6
m_gr = MX.sym('m_gr',1);   % 7
m_or = MX.sym('m_or',1);   % 8

% Algebraic states
p_ai = MX.sym('p_ai',n_w);      % 1-2
p_wh = MX.sym('p_wh',n_w);      % 3-4
p_wi = MX.sym('p_wi',n_w);      % 5-6
p_bh = MX.sym('p_bh',n_w);      % 7-8
rho_ai = MX.sym('rho_ai',n_w);  % 9-10
rho_m = MX.sym('rho_m',n_w);    % 11-12
w_iv = MX.sym('w_iv',n_w);      % 13-14
w_pc = MX.sym('w_pc',n_w);      % 15-16
w_pg = MX.sym('w_pg',n_w);      % 17-18
w_po = MX.sym('w_po',n_w);      % 19-20
w_ro = MX.sym('w_ro',n_w);      % 21-22
w_rg = MX.sym('w_rg',n_w);      % 23-24
p_rh = MX.sym('p_rh',1);        % 25
rho_r =  MX.sym('rho_r',1);     % 26
p_m = MX.sym('p_m',1);          % 27
w_pr = MX.sym('w_pr',1);        % 28
w_to = MX.sym('w_to',1);        % 29
w_tg = MX.sym('w_tg',1);        % 30

% control input
w_gl = MX.sym('w_gl',n_w);

% parameters
p_res = MX.sym('p_res',n_w);
PI = MX.sym('PI',n_w);
GOR = MX.sym('GOR',n_w);
T_a = MX.sym('T_a',n_w);
T_w = MX.sym('T_w',n_w);
T_r = MX.sym('T_r',1);
p_s = MX.sym('p_s',1);

% algebraic equations
f1 = -p_ai.*1e5 + ((R.*T_a./(V_a.*Mw) + 9.81.*H_a./V_a).*m_ga.*1e3) + (Mw./(R.*T_a).*((R.*T_a./(V_a.*Mw) + 9.81.*H_a./V_a).*m_ga.*1e3)).*9.81.*H_a;
f2 = -p_wh.*1e5 + ((R.*T_w./Mw).*(m_gt.*1e3./(L_w.*A_w + L_bh.*A_bh - m_ot.*1e3./rho_o))) - ((m_gt.*1e3+m_ot.*1e3 )./(L_w.*A_w)).*9.81.*H_w/2;
f3 = -p_wi.*1e5 + (p_wh.*1e5 + 9.81./(A_w.*L_w).*max(0,(m_ot.*1e3+m_gt.*1e3-rho_o.*L_bh.*A_bh)).*H_w + 128.*mu_oil.*L_w.*w_pc./(3.14.*D_w.^4.*((m_gt.*1e3 + m_ot.*1e3).*p_wh.*1e5.*Mw.*rho_o)./(m_ot.*1e3.*p_wh.*1e5.*Mw + rho_o.*R.*T_w.*m_gt.*1e3)));
f4 = -p_bh.*1e5 + (p_wi.*1e5 + rho_o.*9.81.*H_bh + 128.*mu_oil.*L_bh.*w_ro./(3.14.*D_bh.^4.*rho_o));
f5 = -rho_ai.*1e2 +(Mw./(R.*T_a).*p_ai.*1e5);
% f6 = -rho_m.*1e2 + (m_gt.*1e3+m_ot.*1e3 )./(L_w.*A_w);
f6 = -rho_m.*1e2 + ((m_gt.*1e3 + m_ot.*1e3).*p_wh.*1e5.*Mw.*rho_o)./(m_ot.*1e3.*p_wh.*1e5.*Mw + rho_o.*R.*T_w.*m_gt.*1e3);
f7 = -w_iv + C_iv.*sqrt(rho_ai.*1e2.*(p_ai.*1e5 - p_wi.*1e5));
f8 = -w_pc + 1.*C_pc.*sqrt(rho_m.*1e2.*(p_wh.*1e5 - p_m.*1e5));
f9 = -w_pg + (m_gt.*1e3./max(1e-3,(m_gt.*1e3+m_ot.*1e3))).*w_pc;
f10 = -w_po + (m_ot.*1e3./max(1e-3,(m_gt.*1e3+m_ot.*1e3))).*w_pc;
f11 = -w_ro + PI.*1e-6.*(p_res.*1e5 - p_bh.*1e5);
f12 = -w_rg.*1e-1 + GOR.*w_ro;  % 23 - 24

f13 = -p_rh.*1e5 + ((R.*T_r./Mw).*(m_gr.*1e3./(L_r.*A_r))) - ((m_gr.*1e3+m_or.*1e3 )./(L_r.*A_r)).*9.81.*H_r/2; %25
% f14 = -rho_r.*1e2 + (m_gr.*1e3+m_or.*1e3 )./(L_r.*A_r); % 26
f14 = -rho_r.*1e2 + ((m_gr.*1e3 + m_or.*1e3).*p_rh.*1e5.*Mw.*rho_ro)./(m_or.*1e3.*p_rh.*1e5.*Mw + rho_ro.*R.*T_r.*m_gr.*1e3);
f15 = -p_m.*1e5 + (p_rh.*1e5 + 9.81./(A_r.*L_r).*(m_or.*1e3+m_gr.*1e3).*H_r + 128.*mu_oil.*L_r.*w_pr./(3.14.*D_r.^4.*((m_gr.*1e3 + m_or.*1e3).*p_rh.*1e5.*Mw.*rho_ro)./(m_or.*1e3.*p_rh.*1e5.*Mw + rho_ro.*R.*T_r.*m_gr.*1e3))) + 0e5;%27
f16 = -w_pr + 1.*C_pr.*sqrt(rho_r.*1e2.*(p_rh.*1e5 - p_s.*1e5)); %28
f17 = -w_to + (m_or.*1e3./(m_gr.*1e3+m_or.*1e3)).*w_pr; % 29
f18 = -w_tg + (m_gr.*1e3./(m_gr.*1e3+m_or.*1e3)).*w_pr; %30


% differential equations
df1 = (w_gl - w_iv).*1e-3;
df2 = (w_iv + w_rg.*1e-1 - w_pg).*1e-3;
df3 = (w_ro - w_po).*1e-3;
df4 = (sum(w_pg) - w_tg).*1e-3;
df5 = (sum(w_po) - w_to).*1e-3;

% Form the DAE system
diff = vertcat(df1,df2,df3,df4,df5);
alg = vertcat(f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18);

% give parameter values
alg = substitute(alg,p_res,par.p_res);
alg = substitute(alg,PI,par.PI);
alg = substitute(alg,p_s,par.p_s);
alg = substitute(alg,T_a,par.T_a);
alg = substitute(alg,T_w,par.T_w);
alg = substitute(alg,T_r,par.T_r);

% concatenate the differential and algebraic states
x_var = vertcat(m_ga,m_gt,m_ot,m_gr,m_or);
z_var = vertcat(p_ai,p_wh,p_wi,p_bh,rho_ai,rho_m,w_iv,w_pc,w_pg,w_po,w_ro,w_rg,p_rh,rho_r,p_m,w_pr,w_to,w_tg);
p_var = vertcat(w_gl,GOR);

qG = sum(w_pg);
L = -(w_to) + 0.15.*sum(w_gl);

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

rf_function = Function('rf_function',{vertcat(x_var,z_var),p_var(1:2),p_var(3:4)},...
                    {vertcat(g{:}),qG},{'x','u','d'},{'g','qG'});
rf  = rootfinder('rf','newton',rf_function);

%% Get sensitivity matrices

qu = rf.jacobian('u','qG');
% Juu = rf.hessian('u','J');

qu = full(qu(vertcat(x_val,z_val),u_val(1:2),u_val(3:4)));
% Juu = full(Juu(x_val,u_val,d_val));
