ControlNet.m

%% Genarate the control net for the spline functions
%% Peng Robinson mixture equation of state
%% 7-point stencil to get the control net
%% Saikat Mukherjee Purdue University November 1 2022

clear all;
close all;

eta = 1e0;
pAmb = 1e5;
theta = 300;
Ru = 8314.32;
Ma = 28.0134;
Mw = 18.0152;
Ra = Ru/Ma;
Rw = Ru/Mw;
theta_c_w = 647.10;
theta_c_a = 126.03;
p_c_w = 22.064e6;
p_c_a = 3.395e6;
omega_a = 0.039;
R = Rw;
bw = 12.8538*p_c_w/Rw/theta_c_w;
a_c_w = 0.457235*(Rw*theta_c_w)^2/p_c_w;
a_c_a = 0.457235*(Ra*theta_c_a)^2/p_c_a;
alpha_w = ( 1.00095 + 0.39222*(1-theta/theta_c_w) - 0.07294*(1-(theta/theta_c_w)^(-1)) + 0.00706*(1-(theta/theta_c_w)^(-2)) )^2;
alpha_a = exp( (0.13280-0.05052*omega_a+0.25948*omega_a^2)*(1-theta/theta_c_a) + 0.81769*log( 1 + ( 0.31355+1.86745*omega_a-0.52604*omega_a^2 )*(1-sqrt(theta/theta_c_a)) ) );
aw = a_c_w*alpha_w;
aa = a_c_a*alpha_a;
delta = -1.70235+0.44338*theta/theta_c_a;

fun   = @(Sol) rootMaxwellPR(Sol,pAmb,theta,R,bw,bw,aw,aa,delta);
Sol0  = [1,840,0.9,0.1];
Sol   = fsolve(fun,Sol0);
rho_v = real(Sol(1));
rho_l = real(Sol(2));
c_v   = real(Sol(3));
c_l   = real(Sol(4));

%% Stencil construction
%% 1. Density
rho1 = rho_v;
rho2 = rho_v + 0.01*(rho_l-rho_v);
rho3 = rho_v + 0.05*(rho_l-rho_v);
rho4 = rho_v + 0.65*(rho_l-rho_v);
rho5 = rho_v + 0.95*(rho_l-rho_v);
rho6 = rho_v + 0.99*(rho_l-rho_v);
rho7 = rho_l;

%% 2. Concentration
c1 = c_v;
c2 = c_v + 0.01*(c_l-c_v);
c3 = c_v + 0.05*(c_l-c_v);
c4 = c_v + 0.65*(c_l-c_v);
c5 = c_v + 0.95*(c_l-c_v);
c6 = c_v + 0.99*(c_l-c_v);
c7 = c_l;

%% 3. Pressure
p11 = ( R*theta.*rho1*bw./(bw-rho1) ...
    + (c1.^2*aa + (1-c1).^2*aw + 2.*c1.*(1-c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho1.^2./(rho1-bw*(1-sqrt(2)))./(rho1-bw*(1+sqrt(2))) );
p12 = ( R*theta.*rho1*bw./(bw-rho1) ...
    + (c2.^2*aa + (1-c2).^2*aw + 2.*c2.*(1-c2)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho1.^2./(rho1-bw*(1-sqrt(2)))./(rho1-bw*(1+sqrt(2))) );
p13 = ( R*theta.*rho1*bw./(bw-rho1) ...
    + (c3.^2*aa + (1-c3).^2*aw + 2.*c3.*(1-c3)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho1.^2./(rho1-bw*(1-sqrt(2)))./(rho1-bw*(1+sqrt(2))) );
p14 = ( R*theta.*rho1*bw./(bw-rho1) ...
    + (c4.^2*aa + (1-c4).^2*aw + 2.*c4.*(1-c4)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho1.^2./(rho1-bw*(1-sqrt(2)))./(rho1-bw*(1+sqrt(2))) );
p15 = ( R*theta.*rho1*bw./(bw-rho1) ...
    + (c5.^2*aa + (1-c5).^2*aw + 2.*c5.*(1-c5)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho1.^2./(rho1-bw*(1-sqrt(2)))./(rho1-bw*(1+sqrt(2))) );
p16 = ( R*theta.*rho1*bw./(bw-rho1) ...
    + (c6.^2*aa + (1-c6).^2*aw + 2.*c6.*(1-c6)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho1.^2./(rho1-bw*(1-sqrt(2)))./(rho1-bw*(1+sqrt(2))) );
p17 = ( R*theta.*rho1*bw./(bw-rho1) ...
    + (c7.^2*aa + (1-c7).^2*aw + 2.*c7.*(1-c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho1.^2./(rho1-bw*(1-sqrt(2)))./(rho1-bw*(1+sqrt(2))) );

p21 = ( R*theta.*rho2*bw./(bw-rho2) ...
    + (c1.^2*aa + (1-c1).^2*aw + 2.*c1.*(1-c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho2.^2./(rho2-bw*(1-sqrt(2)))./(rho2-bw*(1+sqrt(2))) );
p27 = ( R*theta.*rho2*bw./(bw-rho2) ...
    + (c7.^2*aa + (1-c7).^2*aw + 2.*c7.*(1-c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho2.^2./(rho2-bw*(1-sqrt(2)))./(rho2-bw*(1+sqrt(2))) );
p22 = p21 + (c2-c1)*( (2.0*c1*aa - 2.0*(1-c1)*aw + 2.*(1.0-2.0*c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho2.^2./(rho2-bw*(1-sqrt(2)))./(rho2-bw*(1+sqrt(2))) );
p26 = p27 + (c6-c7)*( (2.0*c7*aa - 2.0*(1-c7)*aw + 2.*(1.0-2.0*c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho2.^2./(rho2-bw*(1-sqrt(2)))./(rho2-bw*(1+sqrt(2))) );
p23 = p13 + (rho2-rho1)*( R*theta*bw^2./(bw-rho1).^2 ...
    + (c3.^2.*aa + (1-c3).^2.*aw + 2.*c3.*(1-c3).*sqrt(aa*aw)*(1-delta)).*bw.^2.*...
    ( 2.*rho1.*(rho1-bw*(1-sqrt(2))).*(rho1-bw*(1+sqrt(2))) - rho1.^2.*2.*(rho1-bw) )./(rho1-bw*(1-sqrt(2))).^2./(rho1-bw*(1+sqrt(2))).^2 );
p24 = p14 + (rho2-rho1)*( R*theta*bw^2./(bw-rho1).^2 ...
    + (c4.^2.*aa + (1-c4).^2.*aw + 2.*c4.*(1-c4).*sqrt(aa*aw)*(1-delta)).*bw.^2.*...
    ( 2.*rho1.*(rho1-bw*(1-sqrt(2))).*(rho1-bw*(1+sqrt(2))) - rho1.^2.*2.*(rho1-bw) )./(rho1-bw*(1-sqrt(2))).^2./(rho1-bw*(1+sqrt(2))).^2 );
p25 = p15 + (rho2-rho1)*( R*theta*bw^2./(bw-rho1).^2 ...
    + (c5.^2.*aa + (1-c5).^2.*aw + 2.*c5.*(1-c5).*sqrt(aa*aw)*(1-delta)).*bw.^2.*...
    ( 2.*rho1.*(rho1-bw*(1-sqrt(2))).*(rho1-bw*(1+sqrt(2))) - rho1.^2.*2.*(rho1-bw) )./(rho1-bw*(1-sqrt(2))).^2./(rho1-bw*(1+sqrt(2))).^2 );

p31 = ( R*theta.*rho3*bw./(bw-rho3) ...
    + (c1.^2*aa + (1-c1).^2*aw + 2.*c1.*(1-c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho3.^2./(rho3-bw*(1-sqrt(2)))./(rho3-bw*(1+sqrt(2))) );
p37 = ( R*theta.*rho3*bw./(bw-rho3) ...
    + (c7.^2*aa + (1-c7).^2*aw + 2.*c7.*(1-c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho3.^2./(rho3-bw*(1-sqrt(2)))./(rho3-bw*(1+sqrt(2))) );
p32 = p31 + (c2-c1)*( (2.0*c1*aa - 2.0*(1-c1)*aw + 2.*(1.0-2.0*c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho3.^2./(rho3-bw*(1-sqrt(2)))./(rho3-bw*(1+sqrt(2))) );
p36 = p37 + (c6-c7)*( (2.0*c7*aa - 2.0*(1-c7)*aw + 2.*(1.0-2.0*c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho3.^2./(rho3-bw*(1-sqrt(2)))./(rho3-bw*(1+sqrt(2))) );
p33 = ( R*theta.*rho3*bw./(bw-rho3) ...
    + (c3.^2*aa + (1-c3).^2*aw + 2.*c3.*(1-c3)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho3.^2./(rho3-bw*(1-sqrt(2)))./(rho3-bw*(1+sqrt(2))) )/eta;
p34 = ( R*theta.*rho3*bw./(bw-rho3) ...
    + (c4.^2*aa + (1-c4).^2*aw + 2.*c4.*(1-c4)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho3.^2./(rho3-bw*(1-sqrt(2)))./(rho3-bw*(1+sqrt(2))) )/eta;
p35 = ( R*theta.*rho3*bw./(bw-rho3) ...
    + (c5.^2*aa + (1-c5).^2*aw + 2.*c5.*(1-c5)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho3.^2./(rho3-bw*(1-sqrt(2)))./(rho3-bw*(1+sqrt(2))) )/eta;

p71 = ( R*theta.*rho7*bw./(bw-rho7) ...
    + (c1.^2*aa + (1-c1).^2*aw + 2.*c1.*(1-c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho7.^2./(rho7-bw*(1-sqrt(2)))./(rho7-bw*(1+sqrt(2))) );
p72 = ( R*theta.*rho7*bw./(bw-rho7) ...
    + (c2.^2*aa + (1-c2).^2*aw + 2.*c2.*(1-c2)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho7.^2./(rho7-bw*(1-sqrt(2)))./(rho7-bw*(1+sqrt(2))) );
p73 = ( R*theta.*rho7*bw./(bw-rho7) ...
    + (c3.^2*aa + (1-c3).^2*aw + 2.*c3.*(1-c3)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho7.^2./(rho7-bw*(1-sqrt(2)))./(rho7-bw*(1+sqrt(2))) );
p74 = ( R*theta.*rho7*bw./(bw-rho7) ...
    + (c4.^2*aa + (1-c4).^2*aw + 2.*c4.*(1-c4)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho7.^2./(rho7-bw*(1-sqrt(2)))./(rho7-bw*(1+sqrt(2))) );
p75 = ( R*theta.*rho7*bw./(bw-rho7) ...
    + (c5.^2*aa + (1-c5).^2*aw + 2.*c5.*(1-c5)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho7.^2./(rho7-bw*(1-sqrt(2)))./(rho7-bw*(1+sqrt(2))) );
p76 = ( R*theta.*rho7*bw./(bw-rho7) ...
    + (c6.^2*aa + (1-c6).^2*aw + 2.*c6.*(1-c6)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho7.^2./(rho7-bw*(1-sqrt(2)))./(rho7-bw*(1+sqrt(2))) );
p77 = ( R*theta.*rho7*bw./(bw-rho7) ...
    + (c7.^2*aa + (1-c7).^2*aw + 2.*c7.*(1-c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho7.^2./(rho7-bw*(1-sqrt(2)))./(rho7-bw*(1+sqrt(2))) );

p61 = ( R*theta.*rho6*bw./(bw-rho6) ...
    + (c1.^2*aa + (1-c1).^2*aw + 2.*c1.*(1-c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho6.^2./(rho6-bw*(1-sqrt(2)))./(rho6-bw*(1+sqrt(2))) );
p67 = ( R*theta.*rho6*bw./(bw-rho6) ...
    + (c7.^2*aa + (1-c7).^2*aw + 2.*c7.*(1-c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho6.^2./(rho6-bw*(1-sqrt(2)))./(rho6-bw*(1+sqrt(2))) );
p62 = p61 + (c2-c1)*( (2.0*c1*aa - 2.0*(1-c1)*aw + 2.*(1.0-2.0*c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho6.^2./(rho6-bw*(1-sqrt(2)))./(rho6-bw*(1+sqrt(2))) );
p66 = p67 + (c6-c7)*( (2.0*c7*aa - 2.0*(1-c7)*aw + 2.*(1.0-2.0*c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho6.^2./(rho6-bw*(1-sqrt(2)))./(rho6-bw*(1+sqrt(2))) );
p63 = p73 + (rho6-rho7)*( R*theta*bw^2./(bw-rho7).^2 ...
    + (c3.^2.*aa + (1-c3).^2.*aw + 2.*c3.*(1-c3).*sqrt(aa*aw)*(1-delta)).*bw.^2.*...
    ( 2.*rho7.*(rho7-bw*(1-sqrt(2))).*(rho7-bw*(1+sqrt(2))) - rho7.^2.*2.*(rho7-bw) )./(rho7-bw*(1-sqrt(2))).^2./(rho7-bw*(1+sqrt(2))).^2 );
p64 = p74 + (rho6-rho7)*( R*theta*bw^2./(bw-rho7).^2 ...
    + (c4.^2.*aa + (1-c4).^2.*aw + 2.*c4.*(1-c4).*sqrt(aa*aw)*(1-delta)).*bw.^2.*...
    ( 2.*rho7.*(rho7-bw*(1-sqrt(2))).*(rho7-bw*(1+sqrt(2))) - rho7.^2.*2.*(rho7-bw) )./(rho7-bw*(1-sqrt(2))).^2./(rho7-bw*(1+sqrt(2))).^2 );
p65 = p75 + (rho6-rho7)*( R*theta*bw^2./(bw-rho7).^2 ...
    + (c5.^2.*aa + (1-c5).^2.*aw + 2.*c5.*(1-c5).*sqrt(aa*aw)*(1-delta)).*bw.^2.*...
    ( 2.*rho7.*(rho7-bw*(1-sqrt(2))).*(rho7-bw*(1+sqrt(2))) - rho7.^2.*2.*(rho7-bw) )./(rho7-bw*(1-sqrt(2))).^2./(rho7-bw*(1+sqrt(2))).^2 );

p51 = ( R*theta.*rho5*bw./(bw-rho5) ...
    + (c1.^2*aa + (1-c1).^2*aw + 2.*c1.*(1-c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho5.^2./(rho5-bw*(1-sqrt(2)))./(rho5-bw*(1+sqrt(2))) );
p57 = ( R*theta.*rho5*bw./(bw-rho5) ...
    + (c7.^2*aa + (1-c7).^2*aw + 2.*c7.*(1-c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho5.^2./(rho5-bw*(1-sqrt(2)))./(rho5-bw*(1+sqrt(2))) );
p52 = p51 + (c2-c1)*( (2.0*c1*aa - 2.0*(1-c1)*aw + 2.*(1.0-2.0*c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho5.^2./(rho5-bw*(1-sqrt(2)))./(rho5-bw*(1+sqrt(2))) );
p56 = p57 + (c6-c7)*( (2.0*c7*aa - 2.0*(1-c7)*aw + 2.*(1.0-2.0*c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho5.^2./(rho5-bw*(1-sqrt(2)))./(rho5-bw*(1+sqrt(2))) );
p53 = ( R*theta.*rho5*bw./(bw-rho5) ...
    + (c3.^2*aa + (1-c3).^2*aw + 2.*c3.*(1-c3)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho5.^2./(rho5-bw*(1-sqrt(2)))./(rho5-bw*(1+sqrt(2))) )/eta;
p54 = ( R*theta.*rho5*bw./(bw-rho5) ...
    + (c4.^2*aa + (1-c4).^2*aw + 2.*c4.*(1-c4)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho5.^2./(rho5-bw*(1-sqrt(2)))./(rho5-bw*(1+sqrt(2))) )/eta;
p55 = ( R*theta.*rho5*bw./(bw-rho5) ...
    + (c5.^2*aa + (1-c5).^2*aw + 2.*c5.*(1-c5)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho5.^2./(rho5-bw*(1-sqrt(2)))./(rho5-bw*(1+sqrt(2))) )/eta;

p41 = ( R*theta.*rho4*bw./(bw-rho4) ...
    + (c1.^2*aa + (1-c1).^2*aw + 2.*c1.*(1-c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho4.^2./(rho4-bw*(1-sqrt(2)))./(rho4-bw*(1+sqrt(2))) );
p47 = ( R*theta.*rho4*bw./(bw-rho4) ...
    + (c7.^2*aa + (1-c7).^2*aw + 2.*c7.*(1-c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho4.^2./(rho4-bw*(1-sqrt(2)))./(rho4-bw*(1+sqrt(2))) );
p42 = p41 + (c2-c1)*( (2.0*c1*aa - 2.0*(1-c1)*aw + 2.*(1.0-2.0*c1)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho4.^2./(rho4-bw*(1-sqrt(2)))./(rho4-bw*(1+sqrt(2))) );
p46 = p47 + (c6-c7)*( (2.0*c7*aa - 2.0*(1-c7)*aw + 2.*(1.0-2.0*c7)*sqrt(aa*aw)*(1-delta)).*bw^2.*rho4.^2./(rho4-bw*(1-sqrt(2)))./(rho4-bw*(1+sqrt(2))) );
p43 = p53;
p44 = p54;
p45 = p55;

rho_intermediate = [ rho1,rho1,rho1,rho1,rho1,rho1,rho1,...
    rho2,rho2,rho2,rho2,rho2,rho2,rho2,...
    rho3,rho3,rho3,rho3,rho3,rho3,rho3,...
    rho4,rho4,rho4,rho4,rho4,rho4,rho4,...
    rho5,rho5,rho5,rho5,rho5,rho5,rho5,...
    rho6,rho6,rho6,rho6,rho6,rho6,rho6,...
    rho7,rho7,rho7,rho7,rho7,rho7,rho7];
c_intermediate = [ c1,c2,c3,c4,c5,c6,c7,...
    c1,c2,c3,c4,c5,c6,c7,...
    c1,c2,c3,c4,c5,c6,c7,...
    c1,c2,c3,c4,c5,c6,c7,...
    c1,c2,c3,c4,c5,c6,c7,...
    c1,c2,c3,c4,c5,c6,c7,...
    c1,c2,c3,c4,c5,c6,c7];
p_intermediate = [ p11,p12,p13,p14,p15,p16,p17,...
    p21,p22,p23,p24,p25,p26,p27,...
    p31,p32,p33,p34,p35,p36,p37,...
    p41,p42,p43,p44,p45,p46,p47,...
    p51,p52,p53,p54,p55,p56,p57,...
    p61,p62,p63,p64,p65,p66,p67,...
    p71,p72,p73,p74,p75,p76,p77];

plot3(rho_intermediate,c_intermediate,p_intermediate,'O')

function F = rootMaxwellPR(Sol,pAmb,theta,R,b_v,b_l,aw,aa,delta)
    rho_v = real(Sol(1));
    rho_l = real(Sol(2));
    c_v   = real(Sol(3));
    c_l   = real(Sol(4));
    x_v = c_v; % Simplifying assumption
    x_l = c_l; % Simplifying assumption
    a_v = x_v^2*aa + (1-x_v)^2*aw + 2*x_v*(1-x_v)*sqrt(aa*aw)*(1-delta);
    a_l = x_l^2*aa + (1-x_l)^2*aw + 2*x_l*(1-x_l)*sqrt(aa*aw)*(1-delta);
    da_v = 2*aa*x_v - 2*aw*(1-x_v) + 2*sqrt(aw*aa)*(1-delta)*(1-2*x_v);
    da_l = 2*aa*x_l - 2*aw*(1-x_l) + 2*sqrt(aw*aa)*(1-delta)*(1-2*x_l); 
    psi_vapor  = R*theta*rho_v*log(rho_v/(b_v-rho_v)) - a_v*b_v*rho_v/(2*sqrt(2))*log(abs(rho_v-b_v*(1-sqrt(2)))/abs(rho_v-b_v*(1+sqrt(2)))) ...
        + R*theta*rho_v*( c_v*log(c_v) + (1-c_v)*log(1-c_v) ) ;
    psi_liquid = R*theta*rho_l*log(rho_l/(b_l-rho_l)) - a_l*b_l*rho_l/(2*sqrt(2))*log(abs(rho_l-b_l*(1-sqrt(2)))/abs(rho_l-b_l*(1+sqrt(2)))) ...
        + R*theta*rho_l*( c_l*log(c_l) + (1-c_l)*log(1-c_l) ) ;
    mu_rho_vapor  = R*theta*b_v/(b_v-rho_v) + R*theta*log(rho_v/(b_v-rho_v)) - a_v*b_v/(2*sqrt(2))*log(abs(rho_v-b_v*(1-sqrt(2)))/abs(rho_v-b_v*(1+sqrt(2)))) ...
        + a_v*b_v^2*rho_v/(rho_v-b_v*(1-sqrt(2)))/(rho_v-b_v*(1+sqrt(2))) + R*theta*( c_v*log(c_v)+(1-c_v)*log(1-c_v) );
    mu_rho_liquid = R*theta*b_l/(b_l-rho_l) + R*theta*log(rho_l/(b_l-rho_l)) - a_l*b_l/(2*sqrt(2))*log(abs(rho_l-b_l*(1-sqrt(2)))/abs(rho_l-b_l*(1+sqrt(2)))) ...
        + a_l*b_l^2*rho_l/(rho_l-b_l*(1-sqrt(2)))/(rho_l-b_l*(1+sqrt(2))) + R*theta*( c_l*log(c_l)+(1-c_l)*log(1-c_l) );
    mu_c_vapor  = -b_v*rho_v/(2*sqrt(2))*da_v*log(abs(rho_v-b_v*(1-sqrt(2)))/abs(rho_v-b_v*(1+sqrt(2)))) + R*theta*rho_v*log(c_v/(1-c_v));
    mu_c_liquid = -b_l*rho_l/(2*sqrt(2))*da_l*log(abs(rho_l-b_l*(1-sqrt(2)))/abs(rho_l-b_l*(1+sqrt(2)))) + R*theta*rho_l*log(c_l/(1-c_l));
    F(1) = mu_rho_vapor + (1-c_v)/rho_v*mu_c_vapor - mu_rho_liquid - (1-c_l)/rho_l*mu_c_liquid;
    F(2) = mu_rho_vapor - c_v/rho_v*mu_c_vapor - mu_rho_liquid + c_l/rho_l*mu_c_liquid;
    F(3) = pAmb - rho_l*mu_rho_liquid + psi_liquid;
    F(4) = pAmb - rho_v*mu_rho_vapor + psi_vapor;
end