% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Linerized model of final version New model. Function version.
%
% x1=mg1
% x2=ml1
% x3=mg2
% x4=ml2
% u1=z
% u2=w_G_in
% u3=w_L_in
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Written by:  Knut Åge Meland & Esmaeil Jahanshahi
% Date:         27.10.2010
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
function [A,B,C,D]= v2_new_6d_linmod(x0,y0,u,par)

%%
z1 = u(1);
z2 = u(2);
d1 = u(3);
d2 = u(4);

x1 = x0(1);
x2 = x0(2);
x3 = x0(3);
x4 = x0(4);
x5 = x0(5);
x6 = x0(6);

Pp_nom=y0(4)*10^5;
w_Gin_ss = y0(12);
w_Lin_ss = y0(13);

%%
g=par.g;
R=par.R;
rho_L=par.rho_L;
theta=par.theta; 
r_w=par.r_w;
r1=par.r_p;
r2=par.r_r; 
hc=par.hc;
Aw=par.Aw;
Ap=par.Ap; 
Ar=par.Ar;
Lw=par.Lw;
Lp=par.Lp; 
Lr=par.Lr;
Lh=par.Lh;
Tp=par.Tp;
Tr=par.Tr;
M_G=par.M_G;
P0=par.P_s;
Pr=par.Pr;
Vw=par.Vw;
Vr=par.Vr; 
Vp=par.Vp;

K1 =par.K1;
K2 =par.K2;
K_g =par.K_g;
K_L=par.K_o; 
k_a=par.k_a;

PI = par.PI;
w_nom = par.w_nom;

K_h=par.k_h; %correction factor to finetune model
my=par.visl; %friction factor 



%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculating differensials:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
V_Gw = Vw - x2/rho_L;

dV_Gw__dx1 = 0;
dV_Gw__dx2 = -1/rho_L;
dV_Gw__dx3 = 0;
dV_Gw__dx4 = 0;
dV_Gw__dx5 = 0;
dV_Gw__dx6 = 0;
dV_Gw__dz1 = 0;
dV_Gw__dz2 = 0;
%%
rho_Gw = x1/V_Gw;  % Density of gas at well-head

drho_Gw__dx1 = 1/V_Gw;
drho_Gw__dx2 = -x1*dV_Gw__dx2/V_Gw^2;
drho_Gw__dx3 = 0;
drho_Gw__dx4 = 0;
drho_Gw__dx5 = 0;
drho_Gw__dx6 = 0;
drho_Gw__dz1 = 0;
drho_Gw__dz2 = 0;

%%
P_wh = rho_Gw*R*Tp/M_G; % Pressure at well-head

dP_wh__dx1 = R*Tp*drho_Gw__dx1 / M_G;
dP_wh__dx2 = R*Tp*drho_Gw__dx2 / M_G;
dP_wh__dx3 = R*Tp*drho_Gw__dx3 / M_G;
dP_wh__dx4 = R*Tp*drho_Gw__dx4 / M_G;
dP_wh__dx5 = 0;
dP_wh__dx6 = 0;
dP_wh__dz1 = 0;
dP_wh__dz2 = 0;
dP_wh__dd1 = 0;
dP_wh__dd2 = 0;
%%
rho_mix_w = (x1+x2)/Vw;

drho_mix_w__dx1 = 1/Vw;
drho_mix_w__dx2 = 1/Vw;
drho_mix_w__dx3 = 0; 
drho_mix_w__dx4 = 0;
drho_mix_w__dx5 = 0;
drho_mix_w__dx6 = 0;
drho_mix_w__dz1 = 0; 
drho_mix_w__dz2 = 0; 
drho_mix_w__dd1 = 0;
drho_mix_w__dd2 = 0;
%%
Alpha_Lw_av = x2/(rho_L*Vw);

dAlpha_Lw_av__dx1 = 0;
dAlpha_Lw_av__dx2 = 1/(rho_L*Vw);
dAlpha_Lw_av__dx3 = 0; 
dAlpha_Lw_av__dx4 = 0;
dAlpha_Lw_av__dx5 = 0;
dAlpha_Lw_av__dx6 = 0;
dAlpha_Lw_av__dz1 = 0; 
dAlpha_Lw_av__dz2 = 0; 
dAlpha_Lw_av__dd1 = 0;
dAlpha_Lw_av__dd2 = 0;
%%
Uslw = w_nom/(Aw*rho_L);
Re_w=Uslw*(2*r_w)/my;
Lambda_w=0.0056 + 0.5*Re_w^(-0.32);
Fric_w=0.5*Alpha_Lw_av*Lambda_w*rho_L*Uslw^2*Lw/(2*r_w);
%%
P_bh = P_wh + rho_mix_w*g*Lw + Fric_w;

dP_bh__dx1 = dP_wh__dx1 + drho_mix_w__dx1*g*Lw;
dP_bh__dx2 = dP_wh__dx2 + drho_mix_w__dx2*g*Lw;
dP_bh__dx3 = dP_wh__dx3 + drho_mix_w__dx3*g*Lw; 
dP_bh__dx4 = dP_wh__dx4 + drho_mix_w__dx4*g*Lw;
dP_bh__dx5 = dP_wh__dx5 + drho_mix_w__dx5*g*Lw;
dP_bh__dx6 = dP_wh__dx6 + drho_mix_w__dx6*g*Lw;
dP_bh__dz1 = dP_wh__dz1 + drho_mix_w__dz1*g*Lw; 
dP_bh__dz2 = dP_wh__dz2 + drho_mix_w__dz2*g*Lw; 
dP_bh__dd1 = 0;
dP_bh__dd2 = 0;
%%
Alpha_Lwt =  k_a*Alpha_Lw_av - 1;

dAlpha_Lwt__dx1 = k_a*dAlpha_Lw_av__dx1;
dAlpha_Lwt__dx2 = k_a*dAlpha_Lw_av__dx2;
dAlpha_Lwt__dx3 = k_a*dAlpha_Lw_av__dx3; 
dAlpha_Lwt__dx4 = k_a*dAlpha_Lw_av__dx4;
dAlpha_Lwt__dx5 = k_a*dAlpha_Lw_av__dx5;
dAlpha_Lwt__dx6 = k_a*dAlpha_Lw_av__dx6;
dAlpha_Lwt__dz1 = k_a*dAlpha_Lw_av__dz1; 
dAlpha_Lwt__dz2 = k_a*dAlpha_Lw_av__dz2; 
dAlpha_Lwt__dd1 = k_a*dAlpha_Lw_av__dd1;
dAlpha_Lwt__dd2 = k_a*dAlpha_Lw_av__dd2;

if(Alpha_Lwt>Alpha_Lw_av)
    dAlpha_Lwt__dx1 = dAlpha_Lw_av__dx1;
    dAlpha_Lwt__dx2 = dAlpha_Lw_av__dx2;
    dAlpha_Lwt__dx3 = dAlpha_Lw_av__dx3; 
    dAlpha_Lwt__dx4 = dAlpha_Lw_av__dx4;
    dAlpha_Lwt__dx5 = dAlpha_Lw_av__dx5;
    dAlpha_Lwt__dx6 = dAlpha_Lw_av__dx6;
    dAlpha_Lwt__dz1 = dAlpha_Lw_av__dz1; 
    dAlpha_Lwt__dz2 = dAlpha_Lw_av__dz2; 
    dAlpha_Lwt__dd1 = dAlpha_Lw_av__dd1;
    dAlpha_Lwt__dd2 = dAlpha_Lw_av__dd2;

    
elseif(Alpha_Lwt<0)
    dAlpha_Lwt__dx1 = 0;
    dAlpha_Lwt__dx2 = 0;
    dAlpha_Lwt__dx3 = 0; 
    dAlpha_Lwt__dx4 = 0;
    dAlpha_Lwt__dx5 = 0;
    dAlpha_Lwt__dx6 = 0;
    dAlpha_Lwt__dz1 = 0; 
    dAlpha_Lwt__dz2 = 0; 
    dAlpha_Lwt__dd1 = 0;
    dAlpha_Lwt__dd2 = 0;
end
%%
rho_mixwt = (1-Alpha_Lwt)*rho_Gw + Alpha_Lwt*par.rho_L;

drho_mixwt__dx1 = - dAlpha_Lwt__dx1*rho_Gw + (1 - Alpha_Lwt)*drho_Gw__dx1 + dAlpha_Lwt__dx1*par.rho_L;
drho_mixwt__dx2 = - dAlpha_Lwt__dx2*rho_Gw + (1 - Alpha_Lwt)*drho_Gw__dx2 + dAlpha_Lwt__dx2*par.rho_L;
drho_mixwt__dx3 = - dAlpha_Lwt__dx3*rho_Gw + (1 - Alpha_Lwt)*drho_Gw__dx3 + dAlpha_Lwt__dx3*par.rho_L;
drho_mixwt__dx4 = - dAlpha_Lwt__dx4*rho_Gw + (1 - Alpha_Lwt)*drho_Gw__dx4 + dAlpha_Lwt__dx4*par.rho_L;
drho_mixwt__dx5 = - dAlpha_Lwt__dx5*rho_Gw + (1 - Alpha_Lwt)*drho_Gw__dx5 + dAlpha_Lwt__dx5*par.rho_L;
drho_mixwt__dx6 = - dAlpha_Lwt__dx6*rho_Gw + (1 - Alpha_Lwt)*drho_Gw__dx6 + dAlpha_Lwt__dx6*par.rho_L;
drho_mixwt__dz1 = - dAlpha_Lwt__dz1*rho_Gw + (1 - Alpha_Lwt)*drho_Gw__dz1 + dAlpha_Lwt__dz1*par.rho_L;
drho_mixwt__dz2 = - dAlpha_Lwt__dz2*rho_Gw + (1 - Alpha_Lwt)*drho_Gw__dz2 + dAlpha_Lwt__dz2*par.rho_L;
drho_mixwt__dd1 = 0;
drho_mixwt__dd2 = 0;
%%
Alpha_Gwmt = rho_Gw*(1-Alpha_Lwt)/(rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt);

dAlpha_Gwmt__dx1 = ((drho_Gw__dx1*(1-Alpha_Lwt) - dAlpha_Lwt__dx1*rho_Gw)*(rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)...
    - (drho_Gw__dx1*(1-Alpha_Lwt) - dAlpha_Lwt__dx1*rho_Gw + rho_L*dAlpha_Lwt__dx1)*rho_Gw*(1-Alpha_Lwt))/...
    ((rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)^2);

dAlpha_Gwmt__dx2 = ((drho_Gw__dx2*(1-Alpha_Lwt) - dAlpha_Lwt__dx2*rho_Gw)*(rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)...
    - (drho_Gw__dx2*(1-Alpha_Lwt) - dAlpha_Lwt__dx2*rho_Gw + rho_L*dAlpha_Lwt__dx2)*rho_Gw*(1-Alpha_Lwt))/...
    ((rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)^2);

dAlpha_Gwmt__dx3 = ((drho_Gw__dx3*(1-Alpha_Lwt) - dAlpha_Lwt__dx3*rho_Gw)*(rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)...
    - (drho_Gw__dx3*(1-Alpha_Lwt) - dAlpha_Lwt__dx3*rho_Gw + rho_L*dAlpha_Lwt__dx3)*rho_Gw*(1-Alpha_Lwt))/...
    ((rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)^2);

dAlpha_Gwmt__dx4 = ((drho_Gw__dx4*(1-Alpha_Lwt) - dAlpha_Lwt__dx1*rho_Gw)*(rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)...
    - (drho_Gw__dx4*(1-Alpha_Lwt) - dAlpha_Lwt__dx4*rho_Gw + rho_L*dAlpha_Lwt__dx4)*rho_Gw*(1-Alpha_Lwt))/...
    ((rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)^2);

dAlpha_Gwmt__dx5 = ((drho_Gw__dx5*(1-Alpha_Lwt) - dAlpha_Lwt__dx5*rho_Gw)*(rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)...
    - (drho_Gw__dx5*(1-Alpha_Lwt) - dAlpha_Lwt__dx5*rho_Gw + rho_L*dAlpha_Lwt__dx5)*rho_Gw*(1-Alpha_Lwt))/...
    ((rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)^2);

dAlpha_Gwmt__dx6 = ((drho_Gw__dx6*(1-Alpha_Lwt) - dAlpha_Lwt__dx6*rho_Gw)*(rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)...
    - (drho_Gw__dx6*(1-Alpha_Lwt) - dAlpha_Lwt__dx6*rho_Gw + rho_L*dAlpha_Lwt__dx6)*rho_Gw*(1-Alpha_Lwt))/...
    ((rho_Gw*(1-Alpha_Lwt)+rho_L*Alpha_Lwt)^2);

dAlpha_Gwmt__dz1 = 0;

dAlpha_Gwmt__dz2 = 0;
dAlpha_Gwmt__dd1 = 0;
dAlpha_Gwmt__dd2 = 0;
%%
if (Pr - P_bh)<0 
    disp('The flowdirection is wrong. Check parameters 1');
    return
end
w_r = PI*max(Pr - P_bh,0);


dw_r__dx1 = -PI*dP_bh__dx1;
dw_r__dx2 = -PI*dP_bh__dx2;
dw_r__dx3 = -PI*dP_bh__dx3;
dw_r__dx4 = -PI*dP_bh__dx4;
dw_r__dx5 = -PI*dP_bh__dx5;
dw_r__dx6 = -PI*dP_bh__dx6;
dw_r__dz1 = -0;
dw_r__dz2 = -0;
dw_r__dd1 = -0;
dw_r__dd2 = -0;
%%
rho_Gp_av=Pp_nom*M_G/(R*Tp);

%%
a_Lp_av=(rho_Gp_av*w_Lin_ss)/(rho_Gp_av*w_Lin_ss + rho_L*w_Gin_ss);
%a_Lp_av = par.Alpha_Lp_av;

da_Lp_av__dx1 = 0;
da_Lp_av__dx2 = 0;
da_Lp_av__dx3 = 0;
da_Lp_av__dx4 = 0;
da_Lp_av__dx5 = 0;
da_Lp_av__dx6 = 0;
da_Lp_av__dz1 = 0;
da_Lp_av__dz2 = 0;
da_Lp_av__dd1 = 0; %-(rho_L*rho_Gp_av*w_Lin_ss)/(rho_Gp_av*w_Lin_ss +rho_L*w_Gin_ss)^2;
da_Lp_av__dd2 = 0; %(rho_Gp_av*rho_L*w_Gin_ss)/(rho_Gp_av*w_Lin_ss +rho_L*w_Gin_ss)^2;
%%
h1_av = K_h*hc*a_Lp_av;

dh1_av__dx3 = 0;
dh1_av__dx4 = 0;
dh1_av__dx5 = 0;
dh1_av__dx6 = 0;
dh1_av__dz2 = 0;
dh1_av__dd1 = K_h*hc*da_Lp_av__dd1;
dh1_av__dd2 = K_h*hc*da_Lp_av__dd2;
%%
h1 = h1_av + sin(theta)*(x4 - Vp*rho_L*a_Lp_av)/(Ap*(1-a_Lp_av)*rho_L);

dh1__dx1 = 0;
dh1__dx2 = 0;
dh1__dx3 = 0;
dh1__dx4 = sin(theta)/(Ap*(1-a_Lp_av)*par.rho_L);
dh1__dx5 = 0;
dh1__dx6 = 0;
dh1__dz1 = 0;
dh1__dz2 = 0;
dh1__dd1 = dh1_av__dd1 ...
    +sin(theta)*x4*da_Lp_av__dd1/(Ap*rho_L*(1-a_Lp_av)^2)...
    -(Vp*sin(theta)/Ap)*da_Lp_av__dd1/(1-a_Lp_av)^2;
dh1__dd2 = dh1_av__dd2 ...
    +sin(theta)*x4*da_Lp_av__dd2/(Ap*rho_L*(1-a_Lp_av)^2)...
    -(Vp*sin(theta)/Ap)*da_Lp_av__dd2/(1-a_Lp_av)^2;

%%
V_Gp = Vp -x4/rho_L;

dV_Gp__dx1 = 0;
dV_Gp__dx2 = 0;
dV_Gp__dx3 = 0;
dV_Gp__dx4 = -1/rho_L;
dV_Gp__dx5 = 0;
dV_Gp__dx6 = 0;
dV_Gp__dz2 = 0;
dV_Gp__dd1 = 0;
dV_Gp__dd2 = 0;

%%
rho_Gp = x3/V_Gp;

drho_Gp__dx1 = 0;
drho_Gp__dx2 = 0;
drho_Gp__dx3 = 1/V_Gp;
drho_Gp__dx4 = -x3*dV_Gp__dx4/V_Gp^2;
drho_Gp__dx5 = 0;
drho_Gp__dx6 = 0;
drho_Gp__dz1 = 0;
drho_Gp__dz2 = 0;
drho_Gp__dd1 = 0;
drho_Gp__dd2 = 0;

%%
P_p = rho_Gp*R*Tp/M_G;

dP_p__dx1 = 0;
dP_p__dx2 = 0;
dP_p__dx3 = R*Tp*drho_Gp__dx3 / M_G;
dP_p__dx4 = R*Tp*drho_Gp__dx4 / M_G;
dP_p__dx5 = 0;
dP_p__dx6 = 0;
dP_p__dz1 = 0;
dP_p__dz2 = 0;
dP_p__dd1 = 0;
dP_p__dd2 = 0;

%%
% Flow rate form well to pipeline
if (P_wh-P_p)<0 
    disp('The flowdirection is wrong. Check parameters 2');
    return
end

w_mix_wh = K1*z1*sqrt(rho_mixwt*max(0,P_wh-P_p));

dw_mix_wh__dx1 = K1*z1*(drho_mixwt__dx1*(P_wh-P_p)+rho_mixwt*(dP_wh__dx1-dP_p__dx1))/(2*sqrt(rho_mixwt*(P_wh-P_p)));
dw_mix_wh__dx2 = K1*z1*(drho_mixwt__dx2*(P_wh-P_p)+rho_mixwt*(dP_wh__dx2-dP_p__dx2))/(2*sqrt(rho_mixwt*(P_wh-P_p)));
dw_mix_wh__dx3 = K1*z1*(drho_mixwt__dx3*(P_wh-P_p)+rho_mixwt*(dP_wh__dx3-dP_p__dx3))/(2*sqrt(rho_mixwt*(P_wh-P_p)));
dw_mix_wh__dx4 = K1*z1*(drho_mixwt__dx4*(P_wh-P_p)+rho_mixwt*(dP_wh__dx4-dP_p__dx4))/(2*sqrt(rho_mixwt*(P_wh-P_p)));
dw_mix_wh__dx5 = K1*z1*(drho_mixwt__dx5*(P_wh-P_p)+rho_mixwt*(dP_wh__dx5-dP_p__dx5))/(2*sqrt(rho_mixwt*(P_wh-P_p)));
dw_mix_wh__dx6 = K1*z1*(drho_mixwt__dx6*(P_wh-P_p)+rho_mixwt*(dP_wh__dx6-dP_p__dx6))/(2*sqrt(rho_mixwt*(P_wh-P_p)));
dw_mix_wh__dz1 = K1*sqrt(rho_mixwt*(P_wh-P_p));
dw_mix_wh__dz2 = 0;
dw_mix_wh__dd1 = 0;
dw_mix_wh__dd2 = 0;

%%
w_Gwh = w_mix_wh*Alpha_Gwmt;

dw_Gwh__dx1 = dw_mix_wh__dx1*Alpha_Gwmt + dAlpha_Gwmt__dx1*w_mix_wh;
dw_Gwh__dx2 = dw_mix_wh__dx2*Alpha_Gwmt + dAlpha_Gwmt__dx2*w_mix_wh;
dw_Gwh__dx3 = dw_mix_wh__dx3*Alpha_Gwmt + dAlpha_Gwmt__dx3*w_mix_wh;
dw_Gwh__dx4 = dw_mix_wh__dx4*Alpha_Gwmt + dAlpha_Gwmt__dx4*w_mix_wh;
dw_Gwh__dx5 = dw_mix_wh__dx5*Alpha_Gwmt + dAlpha_Gwmt__dx5*w_mix_wh;
dw_Gwh__dx6 = dw_mix_wh__dx6*Alpha_Gwmt + dAlpha_Gwmt__dx6*w_mix_wh;
dw_Gwh__dz1 = dw_mix_wh__dz1*Alpha_Gwmt + dAlpha_Gwmt__dz1*w_mix_wh;
dw_Gwh__dz2 = 0;
dw_Gwh__dd1 = 0;
dw_Gwh__dd2 = 0;
%%
w_Gin = w_Gwh + d1;

dw_Gin__dx1 = dw_Gwh__dx1;
dw_Gin__dx2 = dw_Gwh__dx2;
dw_Gin__dx3 = dw_Gwh__dx3;
dw_Gin__dx4 = dw_Gwh__dx4;
dw_Gin__dx5 = dw_Gwh__dx5;
dw_Gin__dx6 = dw_Gwh__dx6;
dw_Gin__dz1 = dw_Gwh__dz1;
dw_Gin__dz2 = dw_Gwh__dz2;
dw_Gin__dd1 = 1;
dw_Gin__dd2 = 0;

%%
w_Lwh = w_mix_wh*(1-Alpha_Gwmt);
dw_Lwh__dx1 = dw_mix_wh__dx1*(1-Alpha_Gwmt) - dAlpha_Gwmt__dx1*w_mix_wh;
dw_Lwh__dx2 = dw_mix_wh__dx2*(1-Alpha_Gwmt) - dAlpha_Gwmt__dx2*w_mix_wh;
dw_Lwh__dx3 = dw_mix_wh__dx3*(1-Alpha_Gwmt) - dAlpha_Gwmt__dx3*w_mix_wh;
dw_Lwh__dx4 = dw_mix_wh__dx4*(1-Alpha_Gwmt) - dAlpha_Gwmt__dx4*w_mix_wh;
dw_Lwh__dx5 = dw_mix_wh__dx5*(1-Alpha_Gwmt) - dAlpha_Gwmt__dx5*w_mix_wh;
dw_Lwh__dx6 = dw_mix_wh__dx6*(1-Alpha_Gwmt) - dAlpha_Gwmt__dx6*w_mix_wh;
dw_Lwh__dz1 = dw_mix_wh__dz1*(1-Alpha_Gwmt) - dAlpha_Gwmt__dz1*w_mix_wh;
dw_Lwh__dz2 = 0;
dw_Lwh__dd1 = 0;
dw_Lwh__dd2 = 0;

%%
w_Lin = w_Lwh + d2;
dw_Lin__dx1 = dw_Lwh__dx1;
dw_Lin__dx2 = dw_Lwh__dx2;
dw_Lin__dx3 = dw_Lwh__dx3;
dw_Lin__dx4 = dw_Lwh__dx4;
dw_Lin__dx5 = dw_Lwh__dx5;
dw_Lin__dx6 = dw_Lwh__dx6;
dw_Lin__dz1 = dw_Lwh__dz1;
dw_Lin__dz2 = dw_Lwh__dz2;
dw_Lin__dd1 = 0;
dw_Lin__dd2 = 1;
%%
U_slin_av = w_Lin/(pi*r1^2*rho_L);

dU_slin_av__dx1 = dw_Lin__dx1/(pi*r1^2*rho_L);
dU_slin_av__dx2 = dw_Lin__dx2/(pi*r1^2*rho_L);
dU_slin_av__dx3 = dw_Lin__dx3/(pi*r1^2*rho_L);
dU_slin_av__dx4 = dw_Lin__dx4/(pi*r1^2*rho_L);
dU_slin_av__dx5 = dw_Lin__dx5/(pi*r1^2*rho_L);
dU_slin_av__dx6 = dw_Lin__dx6/(pi*r1^2*rho_L);
dU_slin_av__dz1 = dw_Lin__dz1/(pi*r1^2*rho_L);
dU_slin_av__dz2 = dw_Lin__dz2/(pi*r1^2*rho_L);
dU_slin_av__dd1 = dw_Lin__dd1/(pi*r1^2*rho_L);
dU_slin_av__dd2 = 1 / (pi*r1^2*rho_L );

%%
Re_p = 2*rho_L*U_slin_av*r1 / my ;

dRe_p__dx1 = 2*rho_L*dU_slin_av__dx1*r1 / my ;
dRe_p__dx2 = 2*rho_L*dU_slin_av__dx2*r1 / my ;
dRe_p__dx3 = 2*rho_L*dU_slin_av__dx3*r1 / my ;
dRe_p__dx4 = 2*rho_L*dU_slin_av__dx4*r1 / my ;
dRe_p__dx5 = 2*rho_L*dU_slin_av__dx5*r1 / my ;
dRe_p__dx6 = 2*rho_L*dU_slin_av__dx6*r1 / my ;
dRe_p__dz1 = 2*rho_L*r1*dU_slin_av__dz1 / my ;
dRe_p__dz2 = 2*rho_L*r1*dU_slin_av__dz2 / my ;
dRe_p__dd1 = 2*rho_L*r1* dU_slin_av__dd1 / my;
dRe_p__dd2 = 2*rho_L*r1* dU_slin_av__dd2 / my;

%%
Lambda_p = 0.0056 +0.5*Re_p^(-0.32);

dLambda_p__dx1 = -0.5*0.32*Re_p^(-1.32)*dRe_p__dx1;
dLambda_p__dx2 = -0.5*0.32*Re_p^(-1.32)*dRe_p__dx2;
dLambda_p__dx3 = -0.5*0.32*Re_p^(-1.32)*dRe_p__dx3;
dLambda_p__dx4 = -0.5*0.32*Re_p^(-1.32)*dRe_p__dx4;
dLambda_p__dx5 = -0.5*0.32*Re_p^(-1.32)*dRe_p__dx5;
dLambda_p__dx6 = -0.5*0.32*Re_p^(-1.32)*dRe_p__dx6;
dLambda_p__dz1 = -0.5*0.32*Re_p^(-1.32)*dRe_p__dz1;
dLambda_p__dz2 = -0.5*0.32*Re_p^(-1.32)*dRe_p__dz2;
dLambda_p__dd1 = -0.5*0.32*Re_p^(-1.32)*dRe_p__dd1;
dLambda_p__dd2 = -0.5*0.32*Re_p^(-1.32)*dRe_p__dd2;

%%
DP_fp = a_Lp_av*Lambda_p*rho_L*U_slin_av^2*Lp/(4*r1);

dDP_fp__dx1 = (a_Lp_av*rho_L*U_slin_av^2*Lp*dLambda_p__dx1 +2*a_Lp_av*rho_L*Lambda_p*U_slin_av*Lp*dU_slin_av__dx1 ) / (4*r1);
dDP_fp__dx2 = (a_Lp_av*rho_L*U_slin_av^2*Lp*dLambda_p__dx2 +2*a_Lp_av*rho_L*Lambda_p*U_slin_av*Lp*dU_slin_av__dx2 ) / (4*r1);
dDP_fp__dx3 = (a_Lp_av*rho_L*U_slin_av^2*Lp*dLambda_p__dx3 +2*a_Lp_av*rho_L*Lambda_p*U_slin_av*Lp*dU_slin_av__dx3 ) / (4*r1);
dDP_fp__dx4 = (a_Lp_av*rho_L*U_slin_av^2*Lp*dLambda_p__dx4 +2*a_Lp_av*rho_L*Lambda_p*U_slin_av*Lp*dU_slin_av__dx4 ) / (4*r1);
dDP_fp__dx5 = (a_Lp_av*rho_L*U_slin_av^2*Lp*dLambda_p__dx5 +2*a_Lp_av*rho_L*Lambda_p*U_slin_av*Lp*dU_slin_av__dx5 ) / (4*r1);
dDP_fp__dx6 = (a_Lp_av*rho_L*U_slin_av^2*Lp*dLambda_p__dx6 +2*a_Lp_av*rho_L*Lambda_p*U_slin_av*Lp*dU_slin_av__dx6 ) / (4*r1);
dDP_fp__dz1 = ( Lambda_p*rho_L*U_slin_av^2*Lp*da_Lp_av__dz1  ...
    +a_Lp_av*rho_L*U_slin_av^2*Lp*dLambda_p__dz1 ...
    +2*a_Lp_av*rho_L*Lambda_p*U_slin_av*Lp*dU_slin_av__dz1 ) / (4*r1);
dDP_fp__dz2 = 0;
dDP_fp__dd1 = Lambda_p*rho_L*U_slin_av^2*Lp*da_Lp_av__dd1 / (4*r1);
dDP_fp__dd2 = ( Lambda_p*rho_L*U_slin_av^2*Lp*da_Lp_av__dd2  ...
    +a_Lp_av*rho_L*U_slin_av^2*Lp*dLambda_p__dd2 ...
    +2*a_Lp_av*rho_L*Lambda_p*U_slin_av*Lp*dU_slin_av__dd2 ) / (4*r1);

%%
V_Gr = Vr - x6 / rho_L ; 

dV_Gr__dx1 = 0;
dV_Gr__dx2 = 0;
dV_Gr__dx3 = 0;
dV_Gr__dx4 = 0;
dV_Gr__dx5 = 0;
dV_Gr__dx6 = -1/rho_L ;
dV_Gr__dz1 = 0;
dV_Gr__dz2 = 0;
dV_Gr__dd1 = 0;
dV_Gr__dd2 = 0;

%%
rho_Gr = x5 / V_Gr ; 

drho_Gr__dx1 = 0;
drho_Gr__dx2 = 0;
drho_Gr__dx3 = 0;
drho_Gr__dx4 = 0;
drho_Gr__dx5 = 1/V_Gr;
drho_Gr__dx6 = -x5*dV_Gr__dx6 / V_Gr^2;
drho_Gr__dz1 = 0;
drho_Gr__dz2 = 0;
drho_Gr__dd1 = 0;
drho_Gr__dd2 = 0;

%%
Prt = rho_Gr*R*Tr / M_G; 

dPrt__dx1 = 0;
dPrt__dx2 = 0;
dPrt__dx3 = 0;
dPrt__dx4 = 0;
dPrt__dx5 = R*Tr*drho_Gr__dx5/M_G;
dPrt__dx6 = R*Tr*drho_Gr__dx6/M_G;
dPrt__dz1 = 0;
dPrt__dz2 = 0;
dPrt__dd1 = 0;
dPrt__dd2 = 0;

%%
Psqrt = sqrt((Prt-P0)*1e-5);
dPsqrt__dx1 = 0;
dPsqrt__dx2 = 0;
dPsqrt__dx3 = 0;
dPsqrt__dx4 = 0;
dPsqrt__dx5 = 1e-5*dPrt__dx5/(2*Psqrt);
dPsqrt__dx6 = 1e-5*dPrt__dx6/(2*Psqrt);
dPsqrt__dz1 = 0;
dPsqrt__dz2 = 0;
dPsqrt__dd1 = 0;
dPsqrt__dd2 = 0;
%%
Alpha_Lr_av = x6/(rho_L*Vr);

dAlpha_Lr__dx1 = 0; 
dAlpha_Lr__dx2 = 0;
dAlpha_Lr__dx3 = 0; 
dAlpha_Lr__dx4 = 0;
dAlpha_Lr__dx5 = 0;
dAlpha_Lr__dx6 = 1/(rho_L*Vr);
dAlpha_Lr__dz1 = 0;
dAlpha_Lr__dz2 = 0;
dAlpha_Lr__dd1 = 0;
dAlpha_Lr__dd2 = 0;

%%
rho_m_av = (x5+x6)/Vr;

drho_m_av__dx1 = 0; 
drho_m_av__dx2 = 0;
drho_m_av__dx3 = 0; 
drho_m_av__dx4 = 0;
drho_m_av__dx5 = 1/Vr;
drho_m_av__dx6 = 1/Vr;
drho_m_av__dz1 = 0; 
drho_m_av__dz2 = 0; 
drho_m_av__dd1 = 0; 
drho_m_av__dd2 = 0; 

%%
U_sl2_av = w_Lin/(rho_L*pi*r2^2);
U_sg2_av = w_Gin /(rho_Gr*pi*r2^2);

%%
U_m_av = U_sl2_av +U_sg2_av;

dU_m_av__dx1 = dw_Lin__dx1/(rho_L*pi*r2^2) + (dw_Gin__dx1*(rho_Gr*pi*r2^2)-w_Gin*drho_Gr__dx1*pi*r2^2)/((rho_Gr*pi*r2^2)^2);
dU_m_av__dx2 = dw_Lin__dx2/(rho_L*pi*r2^2) + (dw_Gin__dx2*(rho_Gr*pi*r2^2)-w_Gin*drho_Gr__dx2*pi*r2^2)/((rho_Gr*pi*r2^2)^2);
dU_m_av__dx3 = dw_Lin__dx3/(rho_L*pi*r2^2) + (dw_Gin__dx3*(rho_Gr*pi*r2^2)-w_Gin*drho_Gr__dx3*pi*r2^2)/((rho_Gr*pi*r2^2)^2);
dU_m_av__dx4 = dw_Lin__dx4/(rho_L*pi*r2^2) + (dw_Gin__dx4*(rho_Gr*pi*r2^2)-w_Gin*drho_Gr__dx4*pi*r2^2)/((rho_Gr*pi*r2^2)^2);
dU_m_av__dx5 = dw_Lin__dx5/(rho_L*pi*r2^2) + (dw_Gin__dx5*(rho_Gr*pi*r2^2)-w_Gin*drho_Gr__dx5*pi*r2^2)/((rho_Gr*pi*r2^2)^2);
dU_m_av__dx6 = dw_Lin__dx6/(rho_L*pi*r2^2) + (dw_Gin__dx6*(rho_Gr*pi*r2^2)-w_Gin*drho_Gr__dx6*pi*r2^2)/((rho_Gr*pi*r2^2)^2);
dU_m_av__dz1 = dw_Lin__dz1/(rho_L*pi*r2^2) + (dw_Gin__dz1*(rho_Gr*pi*r2^2)-w_Gin*drho_Gr__dz1*pi*r2^2)/((rho_Gr*pi*r2^2)^2);
dU_m_av__dz2 = 0;
dU_m_av__dd1 = dw_Lin__dd1/ (rho_L*pi*r2^2) + dw_Gin__dd1/ (rho_Gr*pi*r2^2); 
dU_m_av__dd2 = dw_Lin__dd2/ (rho_L*pi*r2^2) + dw_Gin__dd2/ (rho_Gr*pi*r2^2); 

%%
Re_r= 2*rho_m_av *U_m_av* r2 / my;

dRe_r__dx1 = 2*rho_m_av * r2* dU_m_av__dx1 / my;
dRe_r__dx2 = 2*rho_m_av * r2* dU_m_av__dx2 / my;
dRe_r__dx3 = 2*rho_m_av * r2* dU_m_av__dx3 / my;
dRe_r__dx4 = 2*rho_m_av * r2* dU_m_av__dx4 / my;
dRe_r__dx5 = 2*rho_m_av * r2* dU_m_av__dx5 / my  ...
    + 2 *U_m_av* r2* drho_m_av__dx5 / my;
dRe_r__dx6 = 2*rho_m_av * r2* dU_m_av__dx6 / my  ...
    + 2 *U_m_av* r2* drho_m_av__dx6 / my;
dRe_r__dz1 = 2*rho_m_av * r2* dU_m_av__dz1 / my;
dRe_r__dz2 = 0;
dRe_r__dd1 = 2*rho_m_av * r2* dU_m_av__dd1 / my;
dRe_r__dd2 = 2*rho_m_av * r2* dU_m_av__dd2 / my;

%%
Lambda_r = 0.0056 + 0.5*Re_r^(-0.32);

dLambda_r__dx1 = -0.5*0.32*Re_r^(-1.32)*dRe_r__dx1;
dLambda_r__dx2 = -0.5*0.32*Re_r^(-1.32)*dRe_r__dx2;
dLambda_r__dx3 = -0.5*0.32*Re_r^(-1.32)*dRe_r__dx3;
dLambda_r__dx4 = -0.5*0.32*Re_r^(-1.32)*dRe_r__dx4;
dLambda_r__dx5 = -0.5*0.32*Re_r^(-1.32)*dRe_r__dx5;
dLambda_r__dx6 = -0.5*0.32*Re_r^(-1.32)*dRe_r__dx6;
dLambda_r__dz1 = -0.5*0.32*Re_r^(-1.32)*dRe_r__dz1;
dLambda_r__dz2 = 0;
dLambda_r__dd1 = -0.5*0.32*Re_r^(-1.32)*dRe_r__dd1;
dLambda_r__dd2 = -0.5*0.32*Re_r^(-1.32)*dRe_r__dd2;

%% 
DP_fr = Alpha_Lr_av * Lambda_r *rho_m_av *U_m_av^2 * (Lr + Lh) / (4*r2);

dDP_fr__dx1 = ( 2*Alpha_Lr_av * Lambda_r *rho_m_av *U_m_av * (Lr + Lh)*dU_m_av__dx1 ...
    + Alpha_Lr_av * Lambda_r *U_m_av^2 * (Lr + Lh) * drho_m_av__dx1 ...
    +Alpha_Lr_av *rho_m_av *U_m_av^2 * (Lr + Lh) *dLambda_r__dx1 ...
    +Lambda_r *rho_m_av *U_m_av^2 * (Lr + Lh)*dAlpha_Lr__dx1 )/ (4*r2);

dDP_fr__dx2 = ( 2*Alpha_Lr_av * Lambda_r *rho_m_av *U_m_av * (Lr + Lh)*dU_m_av__dx2 ...
    + Alpha_Lr_av * Lambda_r *U_m_av^2 * (Lr + Lh) * drho_m_av__dx2 ...
    +Alpha_Lr_av *rho_m_av *U_m_av^2 * (Lr + Lh) *dLambda_r__dx2 ...
    +Lambda_r *rho_m_av *U_m_av^2 * (Lr + Lh)*dAlpha_Lr__dx2 )/ (4*r2);

dDP_fr__dx3 = ( 2*Alpha_Lr_av * Lambda_r *rho_m_av *U_m_av * (Lr + Lh)*dU_m_av__dx3 ...
    + Alpha_Lr_av * Lambda_r *U_m_av^2 * (Lr + Lh) * drho_m_av__dx3 ...
    +Alpha_Lr_av *rho_m_av *U_m_av^2 * (Lr + Lh) *dLambda_r__dx3 ...
    +Lambda_r *rho_m_av *U_m_av^2 * (Lr + Lh)*dAlpha_Lr__dx3 )/ (4*r2);


dDP_fr__dx4 = ( 2*Alpha_Lr_av * Lambda_r *rho_m_av *U_m_av * (Lr + Lh)*dU_m_av__dx4 ...
    + Alpha_Lr_av * Lambda_r *U_m_av^2 * (Lr + Lh) * drho_m_av__dx4 ...
    +Alpha_Lr_av *rho_m_av *U_m_av^2 * (Lr + Lh) *dLambda_r__dx4 ...
    +Lambda_r *rho_m_av *U_m_av^2 * (Lr + Lh)*dAlpha_Lr__dx4 )/ (4*r2);

dDP_fr__dx5 = ( 2*Alpha_Lr_av * Lambda_r *rho_m_av *U_m_av * (Lr + Lh)*dU_m_av__dx5 ...
    + Alpha_Lr_av * Lambda_r *U_m_av^2 * (Lr + Lh) * drho_m_av__dx5 ...
    +Alpha_Lr_av *rho_m_av *U_m_av^2 * (Lr + Lh) *dLambda_r__dx5 ...
    +Lambda_r *rho_m_av *U_m_av^2 * (Lr + Lh)*dAlpha_Lr__dx5 )/ (4*r2);

dDP_fr__dx6 = ( 2*Alpha_Lr_av * Lambda_r *rho_m_av *U_m_av * (Lr + Lh)*dU_m_av__dx6 ...
    + Alpha_Lr_av * Lambda_r *U_m_av^2 * (Lr + Lh) * drho_m_av__dx6 ...
    +Alpha_Lr_av *rho_m_av *U_m_av^2 * (Lr + Lh) *dLambda_r__dx6 ...
    +Lambda_r *rho_m_av *U_m_av^2 * (Lr + Lh)*dAlpha_Lr__dx6 )/ (4*r2);

dDP_fr__dz1 = ( 2*Alpha_Lr_av * Lambda_r *rho_m_av *U_m_av * (Lr + Lh)*dU_m_av__dz1 ...
    +Alpha_Lr_av *rho_m_av *U_m_av^2 * (Lr + Lh) *dLambda_r__dz1 ) / (4*r2);

dDP_fr__dz2 = 0;

dDP_fr__dd1 = ( 2*Alpha_Lr_av * Lambda_r *rho_m_av *U_m_av * (Lr + Lh)*dU_m_av__dd1 ...
    +Alpha_Lr_av *rho_m_av *U_m_av^2 * (Lr + Lh) *dLambda_r__dd1 ) / (4*r2);

dDP_fr__dd2 = ( 2*Alpha_Lr_av * Lambda_r *rho_m_av *U_m_av * (Lr + Lh)*dU_m_av__dd2 ...
    +Alpha_Lr_av *rho_m_av *U_m_av^2 * (Lr + Lh) *dLambda_r__dd2 ) / (4*r2);



%%
DP_r = rho_m_av*g*Lr + DP_fr;
dDP_r__dx1 = g*Lr*drho_m_av__dx1 + dDP_fr__dx1;
dDP_r__dx2 = g*Lr*drho_m_av__dx2 + dDP_fr__dx2;
dDP_r__dx3 = g*Lr*drho_m_av__dx3 + dDP_fr__dx3;
dDP_r__dx4 = g*Lr*drho_m_av__dx4 + dDP_fr__dx4;
dDP_r__dx5 = g*Lr*drho_m_av__dx5 + dDP_fr__dx5;
dDP_r__dx6 = g*Lr*drho_m_av__dx6 + dDP_fr__dx6;
dDP_r__dz1 = g*Lr*drho_m_av__dz1 + dDP_fr__dz1;
dDP_r__dz2 = g*Lr*drho_m_av__dz2 + dDP_fr__dz2;
dDP_r__dd1 = g*Lr*drho_m_av__dd1 + dDP_fr__dd1;
dDP_r__dd2 = g*Lr*drho_m_av__dd2 + dDP_fr__dd2;

%%
Prb = Prt + DP_r;
dPrb__dx1 = dPrt__dx1 + dDP_r__dx1;
dPrb__dx2 = dPrt__dx2 + dDP_r__dx2;
dPrb__dx3 = dPrt__dx3 + dDP_r__dx3;
dPrb__dx4 = dPrt__dx4 + dDP_r__dx4;
dPrb__dx5 = dPrt__dx5 + dDP_r__dx5;
dPrb__dx6 = dPrt__dx6 + dDP_r__dx6;
dPrb__dz1 = dPrt__dz1 + dDP_r__dz1;
dPrb__dz2 = dPrt__dz2 + dDP_r__dz2;
dPrb__dd1 = dPrt__dd1 + dDP_r__dd1;
dPrb__dd2 = dPrt__dd2 + dDP_r__dd2;

%%
DP_G = P_p - Prb;

dDP_G__dx1 = dP_p__dx1 - dPrb__dx1;
dDP_G__dx2 = dP_p__dx2 - dPrb__dx2; 
dDP_G__dx3 = dP_p__dx3 - dPrb__dx3;
dDP_G__dx4 = dP_p__dx4 - dPrb__dx4; 
dDP_G__dx5 = dP_p__dx5 - dPrb__dx5;
dDP_G__dx6 = dP_p__dx6 - dPrb__dx6;
dDP_G__dz1 = dP_p__dz1 - dPrb__dz1;
dDP_G__dz2 = 0;
dDP_G__dd1 = dP_p__dd1 - dPrb__dd1;
dDP_G__dd2 = dP_p__dd2 - dPrb__dd2;

%%
DP_L = P_p -DP_fp +rho_L*g*h1 - Prt -rho_m_av*g *Lr - DP_fr;

dDP_L__dx1 = dP_p__dx1 - dDP_fp__dx1 + rho_L*g*dh1__dx1 - dPrt__dx1 - g* Lr * drho_m_av__dx1 - dDP_fr__dx1;
dDP_L__dx2 = dP_p__dx2 - dDP_fp__dx2 + rho_L*g*dh1__dx2 - dPrt__dx2 - g* Lr * drho_m_av__dx2 - dDP_fr__dx2;
dDP_L__dx3 = dP_p__dx3 - dDP_fp__dx3 + rho_L*g*dh1__dx3 - dPrt__dx3 - g* Lr * drho_m_av__dx3 - dDP_fr__dx3;
dDP_L__dx4 = dP_p__dx4 - dDP_fp__dx4 + rho_L*g*dh1__dx4 - dPrt__dx4 - g* Lr * drho_m_av__dx4 - dDP_fr__dx4;
dDP_L__dx5 = dP_p__dx5 - dDP_fp__dx5 + rho_L*g*dh1__dx5 - dPrt__dx5 - g* Lr * drho_m_av__dx5 - dDP_fr__dx5;
dDP_L__dx6 = dP_p__dx6 - dDP_fp__dx6 + rho_L*g*dh1__dx6 - dPrt__dx6 - g* Lr * drho_m_av__dx6 - dDP_fr__dx6;
dDP_L__dz1 = dP_p__dz1 - dDP_fp__dx1 + rho_L*g*dh1__dz1 - dPrt__dz1 - g* Lr * drho_m_av__dz1 - dDP_fr__dz1;
dDP_L__dz2 = 0;
dDP_L__dd1 = -dDP_fp__dd1 + rho_L*g*dh1__dd1 - dDP_fr__dd1;
dDP_L__dd2 = -dDP_fp__dd2 + rho_L*g*dh1__dd2 - dDP_fr__dd2;

%%
if h1<hc 
    A_G = (Ap/hc^2)*(hc - h1)^2; 

    dA_G__dx1 = (-2*Ap/hc^2)*(hc - h1)*dh1__dx1;
    dA_G__dx2 = (-2*Ap/hc^2)*(hc - h1)*dh1__dx2;
    dA_G__dx3 = (-2*Ap/hc^2)*(hc - h1)*dh1__dx3;
    dA_G__dx4 = (-2*Ap/hc^2)*(hc - h1)*dh1__dx4;
    dA_G__dx5 = 0;
    dA_G__dx6 = 0;
    dA_G__dz1 = (-2*Ap/hc^2)*(hc - h1)*dh1__dz1;
    dA_G__dz2 = 0;
    dA_G__dd1 = (-2*Ap/hc^2)*(hc - h1)*dh1__dd1;
    dA_G__dd2 = (-2*Ap/hc^2)*(hc - h1)*dh1__dd2;
    
else
    A_G=0;
    dA_G__dx1 = 0;
    dA_G__dx2 = 0;
    dA_G__dx3 = 0;
    dA_G__dx4 = 0;
    dA_G__dx5 = 0;
    dA_G__dx6 = 0;
    dA_G__dz2 = 0;
    dA_G__dd1 = 0;
    dA_G__dd2 = 0;
    
end

%%
A_L = Ap - A_G;

dA_L__dx1 = -dA_G__dx1;
dA_L__dx2 = -dA_G__dx2;
dA_L__dx3 = -dA_G__dx3;
dA_L__dx4 = -dA_G__dx4;
dA_L__dx5 = 0;
dA_L__dx6 = 0;
dA_L__dz1 = 0;
dA_L__dz2 = 0;
dA_L__dd1 = -dA_G__dd1;
dA_L__dd2 = -dA_G__dd2;

%%
if h1<hc
    if DP_G<0
        disp('The flowdirection is wrong. Check parameters');
        return
    end
    w_G_lp = K_g*A_G*sqrt(rho_Gp*DP_G); 

    dw_G_lp__dx1 = K_g*dA_G__dx1*sqrt(rho_Gp*DP_G)...
                + K_g*A_G*(DP_G*drho_Gp__dx1 +rho_Gp*dDP_G__dx1)/(2*sqrt(rho_Gp*DP_G));
    dw_G_lp__dx2 = K_g*dA_G__dx2*sqrt(rho_Gp*DP_G)...
                + K_g*A_G*(DP_G*drho_Gp__dx2 +rho_Gp*dDP_G__dx2)/(2*sqrt(rho_Gp*DP_G));
    dw_G_lp__dx3 = K_g*dA_G__dx3*sqrt(rho_Gp*DP_G)...
                + K_g*A_G*(DP_G*drho_Gp__dx3 +rho_Gp*dDP_G__dx3)/(2*sqrt(rho_Gp*DP_G));
    dw_G_lp__dx4 = K_g*dA_G__dx4*sqrt(rho_Gp*DP_G)...
                + K_g*A_G*(DP_G*drho_Gp__dx4 +rho_Gp*dDP_G__dx4)/(2*sqrt(rho_Gp*DP_G));
    dw_G_lp__dx5 = K_g*dA_G__dx5*sqrt(rho_Gp*DP_G)...
                + K_g*A_G*(DP_G*drho_Gp__dx5 +rho_Gp*dDP_G__dx5)/(2*sqrt(rho_Gp*DP_G));
    dw_G_lp__dx6 = K_g*dA_G__dx6*sqrt(rho_Gp*DP_G)...
                + K_g*A_G*(DP_G*drho_Gp__dx6 +rho_Gp*dDP_G__dx6)/(2*sqrt(rho_Gp*DP_G));
    dw_G_lp__dz1 = K_g*dA_G__dz1*sqrt(rho_Gp*DP_G)...
                + K_g*A_G*(DP_G*drho_Gp__dz1 +rho_Gp*dDP_G__dz1)/(2*sqrt(rho_Gp*DP_G));
    dw_G_lp__dz2 = 0;
    dw_G_lp__dd1 = K_g*A_G*rho_Gp*dDP_G__dd1/(2*sqrt(rho_Gp*DP_G))...
        +K_g*dA_G__dd1*sqrt(rho_Gp*DP_G);
    dw_G_lp__dd2 = K_g*A_G*rho_Gp*dDP_G__dd2/(2*sqrt(rho_Gp*DP_G))...
        +K_g*dA_G__dd2*sqrt(rho_Gp*DP_G);
else
    w_G_lp = 0;

    dw_G_lp__dx1 = 0;
    dw_G_lp__dx2 = 0;
    dw_G_lp__dx3 = 0;
    dw_G_lp__dx4 = 0;
    dw_G_lp__dx5 = 0;
    dw_G_lp__dx6 = 0;
    dw_G_lp__dz1 = 0;
    dw_G_lp__dz2 = 0;
    dw_G_lp__dd1 = 0;
    dw_G_lp__dd2 = 0;
end

%%
if DP_L<0
    disp('The flowdirection is wrong. Check parameters 2');
    return
end
    
w_L_lp = K_L*A_L*sqrt(rho_L*DP_L);

dw_L_lp__dx1 = K_L*A_L*rho_L*dDP_L__dx1/(2*sqrt(rho_L*DP_L)) ...
    +K_L*dA_L__dx1*sqrt(rho_L*DP_L);
dw_L_lp__dx2 = K_L*A_L*rho_L*dDP_L__dx2/(2*sqrt(rho_L*DP_L)) ...
    +K_L*dA_L__dx2*sqrt(rho_L*DP_L);
dw_L_lp__dx3 = K_L*A_L*rho_L*dDP_L__dx3/(2*sqrt(rho_L*DP_L)) ...
    +K_L*dA_L__dx3*sqrt(rho_L*DP_L);
dw_L_lp__dx4 = K_L*A_L*rho_L*dDP_L__dx4/(2*sqrt(rho_L*DP_L)) ...
    +K_L*dA_L__dx4*sqrt(rho_L*DP_L);
dw_L_lp__dx5 = K_L*A_L*rho_L*dDP_L__dx5/(2*sqrt(rho_L*DP_L)) ...
    +K_L*dA_L__dx5*sqrt(rho_L*DP_L);
dw_L_lp__dx6 = K_L*A_L*rho_L*dDP_L__dx6/(2*sqrt(rho_L*DP_L)) ...
    +K_L*dA_L__dx6*sqrt(rho_L*DP_L);
dw_L_lp__dz1 = K_L*A_L*rho_L*dDP_L__dz1/(2*sqrt(rho_L*DP_L)) ...
    +K_L*dA_L__dz1*sqrt(rho_L*DP_L);
dw_L_lp__dz2 = 0;
dw_L_lp__dd1 = K_L*A_L*rho_L*dDP_L__dd1/(2*sqrt(rho_L*DP_L)) ...
    +K_L*dA_L__dd1*sqrt(rho_L*DP_L);
dw_L_lp__dd2 = K_L*A_L*rho_L*dDP_L__dd2/(2*sqrt(rho_L*DP_L)) ...
    +K_L*dA_L__dd2*sqrt(rho_L*DP_L);

%%
Alpha_Lt = 2*x6 / (Vr*rho_L) - A_L/Ap ;

dAlpha_Lt__dx1 = -dA_L__dx1/Ap;
dAlpha_Lt__dx2 = -dA_L__dx2/Ap;
dAlpha_Lt__dx3 = -dA_L__dx3/Ap;
dAlpha_Lt__dx4 = -dA_L__dx4/Ap;
dAlpha_Lt__dx5 = 0;
dAlpha_Lt__dx6 = 2/(Vr*rho_L);
dAlpha_Lt__dz1 = -dA_L__dz1/Ap;
dAlpha_Lt__dz2 = 0;
dAlpha_Lt__dd1 = -dA_L__dd1/Ap;
dAlpha_Lt__dd2 = -dA_L__dd2/Ap;

%%
rho_t=Alpha_Lt*rho_L +(1-Alpha_Lt)*rho_Gr;

drho_t__dx1 = (rho_L-rho_Gr)*dAlpha_Lt__dx1 + (1-Alpha_Lt)*drho_Gr__dx1;
drho_t__dx2 = (rho_L-rho_Gr)*dAlpha_Lt__dx2 + (1-Alpha_Lt)*drho_Gr__dx2;
drho_t__dx3 = (rho_L-rho_Gr)*dAlpha_Lt__dx3 + (1-Alpha_Lt)*drho_Gr__dx3;
drho_t__dx4 = (rho_L-rho_Gr)*dAlpha_Lt__dx4 + (1-Alpha_Lt)*drho_Gr__dx4;
drho_t__dx5 = (rho_L-rho_Gr)*dAlpha_Lt__dx5 + (1-Alpha_Lt)*drho_Gr__dx5;
drho_t__dx6 = (rho_L-rho_Gr)*dAlpha_Lt__dx6 + (1-Alpha_Lt)*drho_Gr__dx6;
drho_t__dz1 = (rho_L-rho_Gr)*dAlpha_Lt__dz1 + (1-Alpha_Lt)*drho_Gr__dz1;
drho_t__dz2 = 0;
drho_t__dd1 = (rho_L-rho_Gr)*dAlpha_Lt__dd1 + (1-Alpha_Lt)*drho_Gr__dd1;
drho_t__dd2 = (rho_L-rho_Gr)*dAlpha_Lt__dd2 + (1-Alpha_Lt)*drho_Gr__dd2;

%%
Alpha_Lmt = Alpha_Lt*rho_L/(Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr);

dAlpha_Lmt__dx1 = (rho_L*dAlpha_Lt__dx1*(Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)...
    - rho_L*Alpha_Lt*dAlpha_Lt__dx1*(rho_L-rho_Gr)) ...
    / (Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)^2;
dAlpha_Lmt__dx2 = (rho_L*dAlpha_Lt__dx2*(Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)...
    - rho_L*Alpha_Lt*dAlpha_Lt__dx2*(rho_L-rho_Gr)) ...
    / (Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)^2;
dAlpha_Lmt__dx3 = (rho_L*dAlpha_Lt__dx3*(Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)...
    - rho_L*Alpha_Lt*dAlpha_Lt__dx3*(rho_L-rho_Gr)) ...
    / (Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)^2;
dAlpha_Lmt__dx4 = (rho_L*dAlpha_Lt__dx4*(Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)...
    - rho_L*Alpha_Lt*dAlpha_Lt__dx4*(rho_L-rho_Gr)) ...
    / (Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)^2 ;
dAlpha_Lmt__dx5 = -rho_L*Alpha_Lt*(1-Alpha_Lt)*drho_Gr__dx5 /(Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)^2;
dAlpha_Lmt__dx6 = (rho_L*dAlpha_Lt__dx6*(Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)...
    - rho_L*Alpha_Lt*dAlpha_Lt__dx6*(rho_L-rho_Gr)) ...
    / (Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)^2 ...
    -rho_L*Alpha_Lt*(1-Alpha_Lt)*drho_Gr__dx6 ...
    /(Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)^2;
dAlpha_Lmt__dz1 = (rho_L*dAlpha_Lt__dz1*(Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)...
    - rho_L*Alpha_Lt*dAlpha_Lt__dz1*(rho_L-rho_Gr)) ...
    / (Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)^2 ;
dAlpha__Lmt__dz2 = 0;
dAlpha_Lmt__dd1 = (rho_L*dAlpha_Lt__dd1*(Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)...
    - rho_L*Alpha_Lt*dAlpha_Lt__dd1*(rho_L-rho_Gr)) ...
    / (Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)^2 ;
dAlpha_Lmt__dd2 = (rho_L*dAlpha_Lt__dd2*(Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)...
    - rho_L*Alpha_Lt*dAlpha_Lt__dd2*(rho_L-rho_Gr)) ...
    / (Alpha_Lt*rho_L + (1-Alpha_Lt)*rho_Gr)^2 ;

%%
if Prt-P0 <0 
    disp('The flowdirection is wrong. Check parameters 3');
    return
end
w_mix_out = K2*z2*sqrt(rho_t*(Prt-P0)); 

dw_mix__dx1 = K2*z2*(drho_t__dx1*(Prt-P0))/(2*sqrt(rho_t*(Prt-P0)));
dw_mix__dx2 = K2*z2*(drho_t__dx2*(Prt-P0))/(2*sqrt(rho_t*(Prt-P0)));
dw_mix__dx3 = K2*z2*(drho_t__dx3*(Prt-P0))/(2*sqrt(rho_t*(Prt-P0)));
dw_mix__dx4 = K2*z2*(drho_t__dx4*(Prt-P0))/(2*sqrt(rho_t*(Prt-P0)));
dw_mix__dx5 = K2*z2*(drho_t__dx5*(Prt-P0)+rho_t*dPrt__dx5)/(2*sqrt(rho_t*(Prt-P0)));
dw_mix__dx6 = K2*z2*(drho_t__dx6*(Prt-P0)+rho_t*dPrt__dx6)/(2*sqrt(rho_t*(Prt-P0)));
dw_mix__dz1 = 0;
dw_mix__dz2 = K2*sqrt(rho_t*(Prt-P0));
dw_mix__dd1 = K2*z2*(drho_t__dd1*(Prt-P0))/(2*sqrt(rho_t*(Prt-P0)));
dw_mix__dd2 = K2*z2*(drho_t__dd2*(Prt-P0))/(2*sqrt(rho_t*(Prt-P0)));
%%
Q_out = K2*z2*sqrt((Prt-P0)/rho_t); 

dQ__dx1 = K2*z2*(-drho_t__dx1*(Prt-P0))/(2*rho_t*sqrt(rho_t*(Prt-P0)));
dQ__dx2 = K2*z2*(-drho_t__dx2*(Prt-P0))/(2*rho_t*sqrt(rho_t*(Prt-P0)));
dQ__dx3 = K2*z2*(-drho_t__dx3*(Prt-P0))/(2*rho_t*sqrt(rho_t*(Prt-P0)));
dQ__dx4 = K2*z2*(-drho_t__dx4*(Prt-P0))/(2*rho_t*sqrt(rho_t*(Prt-P0)));
dQ__dx5 = K2*z2*(-drho_t__dx5*(Prt-P0)+rho_t*dPrt__dx5)/(2*rho_t*sqrt(rho_t*(Prt-P0)));
dQ__dx6 = K2*z2*(-drho_t__dx6*(Prt-P0)+rho_t*dPrt__dx6)/(2*rho_t*sqrt(rho_t*(Prt-P0)));
dQ__dz1 = K2*z2*(-drho_t__dz1*(Prt-P0))/(2*rho_t*sqrt(rho_t*(Prt-P0)));
dQ__dz2 = K2*sqrt((Prt-P0)/rho_t);
dQ__dd1 = K2*z2*(-drho_t__dd1*(Prt-P0))/(2*rho_t*sqrt(rho_t*(Prt-P0)));
dQ__dd2 = K2*z2*(-drho_t__dd2*(Prt-P0))/(2*rho_t*sqrt(rho_t*(Prt-P0)));
%%
w_L_out = Alpha_Lmt*w_mix_out;

dw_L_out__dx1 = dAlpha_Lmt__dx1*w_mix_out + dw_mix__dx1*Alpha_Lmt;
dw_L_out__dx2 = dAlpha_Lmt__dx2*w_mix_out + dw_mix__dx2*Alpha_Lmt;
dw_L_out__dx3 = dAlpha_Lmt__dx3*w_mix_out + dw_mix__dx3*Alpha_Lmt;
dw_L_out__dx4 = dAlpha_Lmt__dx4*w_mix_out + dw_mix__dx4*Alpha_Lmt;
dw_L_out__dx5 = dAlpha_Lmt__dx5*w_mix_out + dw_mix__dx5*Alpha_Lmt;
dw_L_out__dx6 = dAlpha_Lmt__dx6*w_mix_out + dw_mix__dx6*Alpha_Lmt;
dw_L_out__dz1 = dAlpha_Lmt__dz1*w_mix_out + dw_mix__dz1*Alpha_Lmt;
dw_L_out__dz2 = dw_mix__dz2*Alpha_Lmt;
dw_L_out__dd1 = dAlpha_Lmt__dd1*w_mix_out + dw_mix__dd1*Alpha_Lmt;
dw_L_out__dd2 = dAlpha_Lmt__dd2*w_mix_out + dw_mix__dd2*Alpha_Lmt;
%%
w_G_out = (1 - Alpha_Lmt)*w_mix_out;

dw_G_out__dx1 = -dAlpha_Lmt__dx1*w_mix_out + dw_mix__dx1*(1-Alpha_Lmt);
dw_G_out__dx2 = -dAlpha_Lmt__dx2*w_mix_out + dw_mix__dx2*(1-Alpha_Lmt);
dw_G_out__dx3 = -dAlpha_Lmt__dx3*w_mix_out + dw_mix__dx3*(1-Alpha_Lmt);
dw_G_out__dx4 = -dAlpha_Lmt__dx4*w_mix_out + dw_mix__dx4*(1-Alpha_Lmt);
dw_G_out__dx5 = -dAlpha_Lmt__dx5*w_mix_out + dw_mix__dx5*(1-Alpha_Lmt);
dw_G_out__dx6 = -dAlpha_Lmt__dx6*w_mix_out + dw_mix__dx6*(1-Alpha_Lmt);
dw_G_out__dz1 = -dAlpha_Lmt__dz1*w_mix_out + dw_mix__dz1*(1-Alpha_Lmt);
dw_G_out__dz2 =  dw_mix__dz2*(1-Alpha_Lmt);
dw_G_out__dd1 = -dAlpha_Lmt__dd1*w_mix_out + dw_mix__dd1*(1-Alpha_Lmt);
dw_G_out__dd2 = -dAlpha_Lmt__dd2*w_mix_out + dw_mix__dd2*(1-Alpha_Lmt);
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Finished calculating the differensials
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% State equation:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
dx1 = Alpha_Gwmt*w_r - w_Gwh;
dx2 = (1-Alpha_Gwmt)*w_r - w_Lwh;
dx3 = w_Gin - w_G_lp;
dx4 = w_Lin - w_L_lp;
dx5 = w_G_lp - w_G_out; 
dx6 = w_L_lp -w_L_out;

% jackobian:
%% For A:
dx1__dx1 = dAlpha_Gwmt__dx1*w_r + dw_r__dx1*Alpha_Gwmt - dw_Gwh__dx1;
dx1__dx2 = dAlpha_Gwmt__dx2*w_r + dw_r__dx2*Alpha_Gwmt - dw_Gwh__dx2;
dx1__dx3 = dAlpha_Gwmt__dx3*w_r + dw_r__dx3*Alpha_Gwmt - dw_Gwh__dx3;
dx1__dx4 = dAlpha_Gwmt__dx4*w_r + dw_r__dx4*Alpha_Gwmt - dw_Gwh__dx4;
dx1__dx5 = dAlpha_Gwmt__dx5*w_r + dw_r__dx5*Alpha_Gwmt - dw_Gwh__dx5;
dx1__dx6 = dAlpha_Gwmt__dx6*w_r + dw_r__dx6*Alpha_Gwmt - dw_Gwh__dx6;

dx2__dx1 = -dAlpha_Gwmt__dx1*w_r + dw_r__dx1*(1 - Alpha_Gwmt) - dw_Lwh__dx1;
dx2__dx2 = -dAlpha_Gwmt__dx2*w_r + dw_r__dx2*(1 - Alpha_Gwmt) - dw_Lwh__dx2;
dx2__dx3 = -dAlpha_Gwmt__dx3*w_r + dw_r__dx3*(1 - Alpha_Gwmt) - dw_Lwh__dx3;
dx2__dx4 = -dAlpha_Gwmt__dx4*w_r + dw_r__dx4*(1 - Alpha_Gwmt) - dw_Lwh__dx4;
dx2__dx5 = -dAlpha_Gwmt__dx5*w_r + dw_r__dx5*(1 - Alpha_Gwmt) - dw_Lwh__dx5;
dx2__dx6 = -dAlpha_Gwmt__dx6*w_r + dw_r__dx6*(1 - Alpha_Gwmt) - dw_Lwh__dx6;

dx3__dx1 = dw_Gin__dx1-dw_G_lp__dx1;
dx3__dx2 = dw_Gin__dx2-dw_G_lp__dx2;
dx3__dx3 = dw_Gin__dx3-dw_G_lp__dx3;
dx3__dx4 = dw_Gin__dx4-dw_G_lp__dx4;
dx3__dx5 = dw_Gin__dx5-dw_G_lp__dx5;
dx3__dx6 = dw_Gin__dx6-dw_G_lp__dx6;

dx4__dx1 = dw_Lin__dx1-dw_L_lp__dx1;
dx4__dx2 = dw_Lin__dx2-dw_L_lp__dx2;
dx4__dx3 = dw_Lin__dx3-dw_L_lp__dx3;
dx4__dx4 = dw_Lin__dx4-dw_L_lp__dx4;
dx4__dx5 = dw_Lin__dx5-dw_L_lp__dx5;
dx4__dx6 = dw_Lin__dx6-dw_L_lp__dx6;

dx5__dx1 = dw_G_lp__dx1 - dw_G_out__dx1;
dx5__dx2 = dw_G_lp__dx2 - dw_G_out__dx2;
dx5__dx3 = dw_G_lp__dx3 - dw_G_out__dx3;
dx5__dx4 = dw_G_lp__dx4 - dw_G_out__dx4;
dx5__dx5 = dw_G_lp__dx5 - dw_G_out__dx5;
dx5__dx6 = dw_G_lp__dx6 - dw_G_out__dx6;

dx6__dx1 = dw_L_lp__dx1 - dw_L_out__dx1;
dx6__dx2 = dw_L_lp__dx2 - dw_L_out__dx2;
dx6__dx3 = dw_L_lp__dx3 - dw_L_out__dx3;
dx6__dx4 = dw_L_lp__dx4 - dw_L_out__dx4;
dx6__dx5 = dw_L_lp__dx5 - dw_L_out__dx5;
dx6__dx6 = dw_L_lp__dx6 - dw_L_out__dx6;
%% For B:
dx1__dz1 = dAlpha_Gwmt__dz1*w_r + dw_r__dz1*Alpha_Gwmt - dw_Gwh__dz1;
dx1__dz2 = dAlpha_Gwmt__dz2*w_r + dw_r__dz2*Alpha_Gwmt - dw_Gwh__dz2;
dx1__dd1 = dAlpha_Gwmt__dd1*w_r + dw_r__dd1*Alpha_Gwmt - dw_Gwh__dd1;
dx1__dd2 = dAlpha_Gwmt__dd2*w_r + dw_r__dd2*Alpha_Gwmt - dw_Gwh__dd2;

dx2__dz1 = -dAlpha_Gwmt__dz1*w_r + dw_r__dz1*(1 - Alpha_Gwmt) - dw_Lwh__dz1;
dx2__dz2 = -dAlpha_Gwmt__dz2*w_r + dw_r__dz2*(1 - Alpha_Gwmt) - dw_Lwh__dz2;
dx2__dd1 = -dAlpha_Gwmt__dd1*w_r + dw_r__dd1*(1 - Alpha_Gwmt) - dw_Lwh__dd1;
dx2__dd2 = -dAlpha_Gwmt__dd2*w_r + dw_r__dd2*(1 - Alpha_Gwmt) - dw_Lwh__dd2;

dx3__dz1 = dw_Gin__dz1-dw_G_lp__dz1;
dx3__dz2 = dw_Gin__dz2-dw_G_lp__dz2;
dx3__dd1 = dw_Gin__dd1-dw_G_lp__dd1;
dx3__dd2 = dw_Gin__dd2-dw_G_lp__dd2;

dx4__dz1 = dw_Lin__dz1-dw_L_lp__dz1;
dx4__dz2 = dw_Lin__dz2-dw_L_lp__dz2;
dx4__dd1 = dw_Lin__dd1-dw_L_lp__dd1;
dx4__dd2 = dw_Lin__dd2-dw_L_lp__dd2;

dx5__dz1 = dw_G_lp__dz1 - dw_G_out__dz1;
dx5__dz2 = dw_G_lp__dz2 - dw_G_out__dz2;
dx5__dd1 = dw_G_lp__dd1 - dw_G_out__dd1;
dx5__dd2 = dw_G_lp__dd2 - dw_G_out__dd2;

dx6__dz1 = dw_L_lp__dz1 - dw_L_out__dz1;
dx6__dz2 = dw_L_lp__dz2 - dw_L_out__dz2;
dx6__dd1 = dw_L_lp__dd1 - dw_L_out__dd1;
dx6__dd2 = dw_L_lp__dd2 - dw_L_out__dd2;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Connecting jackobians to form A, B, C and D:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = [ dx1__dx1 dx1__dx2 dx1__dx3 dx1__dx4 dx1__dx5 dx1__dx6 ;
      dx2__dx1 dx2__dx2 dx2__dx3 dx2__dx4 dx2__dx5 dx2__dx6 ; 
      dx3__dx1 dx3__dx2 dx3__dx3 dx3__dx4 dx3__dx5 dx3__dx6 ; 
      dx4__dx1 dx4__dx2 dx4__dx3 dx4__dx4 dx4__dx5 dx4__dx6 ;
      dx5__dx1 dx5__dx2 dx5__dx3 dx5__dx4 dx5__dx5 dx5__dx6 ; 
      dx6__dx1 dx6__dx2 dx6__dx3 dx6__dx4 dx6__dx5 dx6__dx6 ];

B = [ dx1__dz1 dx1__dz2 dx1__dd1 dx1__dd2 ;
      dx2__dz1 dx2__dz2 dx2__dd1 dx2__dd2 ;
      dx3__dz1 dx3__dz2 dx3__dd1 dx3__dd2 ;
      dx4__dz1 dx4__dz2 dx4__dd1 dx4__dd2 ;
      dx5__dz1 dx5__dz2 dx5__dd1 dx5__dd2 ;
      dx6__dz1 dx6__dz2 dx6__dd1 dx6__dd2 ];

C = zeros(11,6);
C(1,:) = 1e-5.*[ dP_bh__dx1 dP_bh__dx2 dP_bh__dx3 dP_bh__dx4 dP_bh__dx5 dP_bh__dx6 ]; %Pa-> bar
C(2,:) = 1e-5.*[ dP_wh__dx1 dP_wh__dx2 dP_wh__dx3 dP_wh__dx4 dP_wh__dx5 dP_wh__dx6 ]; %Pa-> bar
C(3,:) = [ dw_mix_wh__dx1 dw_mix_wh__dx2 dw_mix_wh__dx3 dw_mix_wh__dx4 dw_mix_wh__dx5 dw_mix_wh__dx6 ]; %[kg/s]    
C(4,:) = 1e-5.*[ dP_p__dx1 dP_p__dx2 dP_p__dx3 dP_p__dx4 dP_p__dx5 dP_p__dx6 ]; %Pa-> bar
C(5,:) = 1e-5.*[ dPrb__dx1 dPrb__dx2 dPrb__dx3 dPrb__dx4 dPrb__dx5 dPrb__dx6 ]; %Pa-> bar
C(6,:) = 1e-5.*[ dDP_r__dx1 dDP_r__dx2 dDP_r__dx3 dDP_r__dx4 dDP_r__dx5 dDP_r__dx6 ]; %Pa-> bar
C(7,:) = 1e-5.*[ dPrt__dx1 dPrt__dx2 dPrt__dx3 dPrt__dx4 dPrt__dx5 dPrt__dx6 ]; %Pa-> bar
C(8,:) =  1000*[ dQ__dx1 dQ__dx2 dQ__dx3 dQ__dx4 dQ__dx5 dQ__dx6 ]; %[L/s]
C(9,:) =  [ dw_mix__dx1 dw_mix__dx2 dw_mix__dx3 dw_mix__dx4 dw_mix__dx5 dw_mix__dx6 ]; %[kg/s]            
C(10,:) = [drho_t__dx1 drho_t__dx2 drho_t__dx3 drho_t__dx4 drho_t__dx5 drho_t__dx6];
C(11,:) = [dAlpha_Lt__dx1 dAlpha_Lt__dx2 dAlpha_Lt__dx3 dAlpha_Lt__dx4 dAlpha_Lt__dx5 dAlpha_Lt__dx6];

D = zeros(11,4);
D(1,:) = 1e-5.*[ dP_bh__dz1 dP_bh__dz2 dP_bh__dd1 dP_bh__dd2 ]; %Pa-> bar
D(2,:) = 1e-5.*[ dP_wh__dz1 dP_wh__dz2 dP_wh__dd1 dP_wh__dd2 ]; %Pa-> bar
D(3,:) = [ dw_mix_wh__dz1 dw_mix_wh__dz2 dw_mix_wh__dd1 dw_mix_wh__dd2 ]; %[kg/s]    
D(4,:) = 1e-5.*[ dP_p__dz1 dP_p__dz2 dP_p__dd1 dP_p__dd2 ]; %Pa-> bar
D(5,:) = 1e-5.*[ dPrb__dz1 dPrb__dz2 dPrb__dd1 dPrb__dd2 ]; %Pa-> bar
D(6,:) = 1e-5.*[ dDP_r__dz1 dDP_r__dz2 dDP_r__dd1 dDP_r__dd2 ]; %Pa-> bar
D(7,:) = 1e-5.*[ dPrt__dz1 dPrt__dz2 dPrt__dd1 dPrt__dd2 ]; %Pa-> bar
D(8,:) =  1000*[ dQ__dz1 dQ__dz2 dQ__dd1 dQ__dd2 ]; %[L/s]
D(9,:) =  [ dw_mix__dz1 dw_mix__dz2 dw_mix__dd1 dw_mix__dd2 ]; %[kg/s]            
D(10,:) = [ drho_t__dz1 drho_t__dz2 drho_t__dd1 drho_t__dd2 ];
D(11,:) = [ dAlpha_Lt__dz1 dAlpha_Lt__dz2 dAlpha_Lt__dd1 dAlpha_Lt__dd2 ];
