% *************** Finding Initial States *******************************
function [x0,y0,par] = v1_new_6d_initialize(u,par)

% clc
% clear all
% par = v1_new_6d_parameters();
% z1 = 1;
% z2 = 0.04;
% d1 = 0;
% d2 = 0;
% u = [z1;z2;d1;d2];


%  P1: Pressure in pipeline section
%  P2: Pressure at top of the riser
%  h1: Level of liquid in pipeline

%  x1: Mass of gas in the pipelien (m_G1)
%  x2: Mass of liquid in the pipeline (m_L1)
%  x3: Mass of gas in the riser (m_G2)
%  x4: Mass of liquid in the riser (m_L2)

    
if (u(2)<0.6)
    P_wh = 7019947;
    P_bh = 27854510;
else
    P_wh =6980175;
    P_bh = 27799650;
end

w_r = par.PI*(par.Pr-P_bh);
Alpha_Gwmt = 0.0395;
rho_Gw = P_wh*par.M_G/(par.R*par.Tp);

rho_mix_w = (P_bh - P_wh)/(par.g*par.Lw);
V_Gw = par.Vw*(par.rho_L-rho_mix_w)/(par.rho_L-rho_Gw);
m_Gw_ss = V_Gw*rho_Gw;
m_Lw_ss = (par.Vw - V_Gw)*par.rho_L;

w_Gin=Alpha_Gwmt*w_r;
w_Lin=(1-Alpha_Gwmt)*w_r;

    
if (u(2)<0.6)
    P_p_ss = 7011254;
    P_rt_ss = 5069538;
else
    P_p_ss = 6971484;
    P_rt_ss = 5018719;
end

VLr = par.Ar*par.Lr/2;
m_Lr_ss = VLr*par.rho_L;
VGr = par.Vr - VLr;
rho_Gr = P_rt_ss*par.M_G/(par.R*par.Tr);
m_Gr_ss = VGr*rho_Gr;


%h1_ss = double(ss.h1(1));
rho_Gp_n = P_p_ss*par.M_G/(par.R*par.Tp);
a_new = w_Lin*rho_Gp_n/(w_Lin*rho_Gp_n + w_Gin*par.rho_L);
a_old = a_new - 0.001; % to make sure loop runs at least once
%m_Lp_ss = a_old*par.rho_L*par.Vp;
h_old = par.k_h*a_new*par.hc;
h_new = par.k_h*par.Alpha_Lp_av*par.hc;

opt = optimset('Display','off','TolFun',1e-50,'TolX',1e-50,'MaxIter',1e4,'MaxFunEvals',1e10);
y0(4) = P_p_ss/1e5;
while(abs(a_old - a_new)>1e-4)
    P_p_ss = 1e5*y0(4);
    rho_Gp = P_p_ss*par.M_G/(par.R*par.Tp);
    par.Alpha_Lp_av = w_Lin*rho_Gp/(w_Lin*rho_Gp + w_Gin*par.rho_L);
    m_Lp_ss = par.rho_L*par.Vp*par.Alpha_Lp_av + (h_old - h_new)*par.Ap*par.rho_L*(1-par.Alpha_Lp_av)/sin(par.theta);
    m_Gp_ss = (par.Vp - m_Lp_ss/par.rho_L)*rho_Gp;
    x0 = [m_Gw_ss;m_Lw_ss;m_Gp_ss;m_Lp_ss;m_Gr_ss;m_Lr_ss];
    [x0,fval,exitflag,output,jacobian] = fsolve(@v2_new_6d_model0,x0,opt,u,par);
    y0 = v1_new_6d_model(0,x0,u,'measurements',par);
    m_Gw_ss = x0(1);
    m_Lw_ss = x0(2);
    %m_Lp_ss = x0(4);
    m_Gr_ss = x0(5);
    m_Lr_ss = x0(6);
     rho_Gp = 1e5*y0(4)*par.M_G/(par.R*par.Tp);
     w_Gin = y0(12); w_Lin = y0(13);
     a_old = a_new;
     a_new = w_Lin*rho_Gp/(w_Lin*rho_Gp + w_Gin*par.rho_L);
     h_old = h_new;
     h_new = par.k_h*a_new*par.hc;
     %disp(a_new)
%     disp(y0(12))
%     disp(y0(13))
end
fval
%par.Alpha_Lp_av = a_new;
str1 = ['Pressure at bottom-hole (P_bh):' num2str(y0(1)) ' [bar]' ];
disp(str1)
str2 = ['Pressure at well-head (P_wh):' num2str(y0(2)) ' [bar]' ];
disp(str2)
str3 = ['Inlet mass flow (w_in):' num2str(y0(3)) ' [kg/s]' ];
disp(str3)
str4 = ['Pressure at inlet of pipeline (P_in):' num2str(y0(4)) ' [bar]' ];
disp(str4)
str5 = ['Pressure at riser base (P_rb):' num2str(y0(5)) ' [bar]' ];
disp(str5)
str6 = ['Pressure drop over riser (DP_r):' num2str(y0(6)) ' [bar]' ];
disp(str6)
str7 = ['Pressure at top of riser (P2_t):' num2str(y0(7)) ' [bar]' ];
disp(str7)
str8 = ['Volumetric flow of choke (Q_out):' num2str(y0(8)) ' [L/s]' ];
disp(str8)
str9 = ['Mass flow of choke (w_out):' num2str(y0(9)) ' [kg/s]' ];
disp(str9)
str10 = ['Density at top (Rho_t):' num2str(y0(10)) ' [kg/m3]' ];
disp(str10)
str11 = ['Liquid volume fraction at top (Alpha_t): ' num2str(y0(11)) ' [-]' ];
disp(str11)
disp(y0(14))
end