% *************** Finding Initial States *******************************
function [x0,y0,par] = v1_DiMeglio_initialize(z,wG_in,wL_in,par)

% Simplified model for riser slugging
% By: Florent Di Meglio
% Mines Paristech

% x1: Mass of gas in the elongated bubble (m_G1)
% x2: Mass of gas in the riser (m_G2)
% x3: Mass of liquid in the riser (m_L2)


syms m_G1 m_G2 m_L2


P1 = m_G1*par.R*par.T1/(par.M_G*par.Veb);

V_G2 = par.V2 - (m_L2+par.m_lstill)/par.rho_L;
P2_t = m_G2*par.R*par.T2/(par.M_G*V_G2);

P2_b = P2_t + par.g*sin(par.theta2)*(m_L2+par.m_lstill)/par.A2;

eq1 = (1 - par.ep)*wG_in - par.C_g*(P1-P2_b);


eq2 = wG_in - (m_G2/(m_L2+m_G2))*par.C_c*z*sqrt(P2_t-par.P0);

eq3 = wL_in - (m_L2/(m_L2+m_G2))*par.C_c*z*sqrt(P2_t-par.P0);

ss = solve(eq1,eq2,eq3);
%%
m_G1_ss = double(ss.m_G1(1));
m_G2_ss = double(ss.m_G2(1));
m_L2_ss = double(ss.m_L2(1));

x0 = [m_G1_ss;m_G2_ss;m_L2_ss];

% u = [z;wG_in;wL_in];
% opt = optimset('Display','off','TolFun',1e-100,'TolX',1e-100,'MaxIter',1e10,'MaxFunEvals',1e10);
% x0 = fsolve(@v1_DiMeglio_model,x0,opt,u,par);
P10 = x0(1)*par.R*par.T1/(par.M_G*par.Veb);
V_G2 = par.V2 - (x0(3)+par.m_lstill)/par.rho_L;
P2_T0 = x0(2)*par.R*par.T2/(par.M_G*V_G2);
y0 = [P10;P2_T0];
disp(y0/1e5)
end