function [x_eq,flag,y] = v1_gaslift_initialize(u,par)
% warning('off')
% clear all;
% clc;
% par = v1_gaslift_parameters();
% %% Inputs to model 
% z1= 0.9;           %Valve opening top choke (u1) 
% w_G_a_in=0.5;%2;    %Gas mass inflow to annulus (u2) 
% % Disturbances
% GOR = 0;                     % Mass GOR from reservoir, NB: MASS, increasing gives: lower BH, stable
% P_r = 160;                      % Rservoir pressure [bar]
% P_gs = 150;                     % Gass source pressure [bar]
% P0 = 15;                        % Separator Pressure [bar]
% d = [GOR;P_r;P_gs;P0];         % Disturbance Vector
% 
% u=[z1;w_G_a_in;d]; %Only u(1) and u(2) are MVs

%% ***** Parameter Constants*******
%par=setConstants();
%% **** States **********
syms x1 x2 x3;
% x1: Gas mass in annulus
% x2: Gas mass in riser
% x3: Liquid mass in riser	

%% ***** Inputs ******
u1=u(1);    			%Valve opening top choke

%% **** Disturbances *****
%GOR = u(3);   % Mass GOR from reservior (d1)
P_r = u(4)*1e5;    % Rservior pressure (d2) [bar]
%P_gs = u(5)*1e5;   % Injection gass source pressure (d3) [bar]
P0 = u(6)*1e5;       % Separator Pressure (d4) [bar] 
%w_G_r_in=u(3); 			%Gas mass flow into riser
%w_L_r_in=u(4); 			%Liquid mass flow in to riser
w_G_a_in = 0.8; 			%Gas mass flow into annulus	

%% ********* Pressure in bottom of annulus *********************
P_a_b=(par.R*par.T_a*x1/(par.M_G_a*par.V_a)) + (x1*par.g*par.L_a/par.V_a) ; 

%% *************** Pressure in top of riser **********************%
%Check if volume of liquid is bigger than volume of riser
P_r_t=x2*par.R*par.T_r/(par.M_G_r_t*(par.V_r-(x3/par.rho_L))); 

%% ********** Density of gas in bottom of annulus **********
rho_G_a_b= P_a_b*par.M_G_a/(par.R*par.T_a); 

%% ************* Average mixture density ***************
rho_M_av=(x2+x3)/par.V_r; 

%% *********** Pressure in bottom of riser ******************
Fric = 2.5e5;
P_r_b=P_r_t + (rho_M_av*par.g*par.L_r)+Fric;  
P_bh=P_r_b+par.rho_L*par.g*par.L_bh;

%% Liquid Inflow rate
w_L_r_in=par.PI*(P_r-P_bh);

%% *********  Injection gas flow rate from annulus to riser	*********
w_G_a_out = par.K_a*sqrt(rho_G_a_b*(P_a_b - P_r_b)); 			

%% **********  Mixture flow out of riser ****************************
w_mix_r_out=par.K_r*u1*sqrt(rho_M_av*(P_r_t-P0)); 	        

%% ******** Gas mass flow rate out of riser ************ 
w_G_r_out = x2/(x2+x3)*w_mix_r_out;   

%% ******** Liquid mass flow rate out of riser ************ 
w_L_r_out = x3/(x2+x3)*w_mix_r_out;      							

%% Derivatives
eq1 = simple(w_G_a_in - w_G_a_out);              							
eq2 = simple(w_G_a_out - w_G_r_out);  							
eq3 = simple(w_L_r_in - w_L_r_out);                                        
[x1_bar,x2_bar,x3_bar] = solve(eq1,eq2,eq3, 'MaxDegree', 4);
%[x1_bar,x2_bar,x3_bar] = solve(eq1,eq2,eq3);

for i=1:length(x1_bar)
    if (double(x1_bar(i,1))>0 && double(x2_bar(i,1))>0 && double(x3_bar(i,1))>0)
        x1_hat = double(x1_bar(i,1));
        x2_hat = double(x2_bar(i,1));
        x3_hat = double(x3_bar(i,1));
        x0 = [x1_hat;x2_hat;x3_hat];
        break
    end
end
%disp(x0)
% sub1 = subs(eq1,x1,x0(1));
% sub1 = subs(sub1,x2,x0(2));
% sub1 = subs(sub1,x3,x0(3))
% 
% sub2 = subs(eq2,x1,x0(1));
% sub2 = subs(sub2,x2,x0(2));
% sub2 = subs(sub2,x3,x0(3))
% 
% sub3 = subs(eq3,x2,x0(2));
% sub3 = subs(sub3,x3,x0(3))
%% Options for f-solve
opt=optimset;
opt.MaxIter=2e5;
opt.TolFun=1e-17;
opt.MaxFunEvals=1000000;
opt.TolX=1e-17;
opt.Display='off';
%opt.Algorithm='trust-region-dogleg';

%[x_start,fval,flag,output,jacobian] = fsolve(@model0,x_start,opt,u,par);
[x_eq,fval,flag,output,jacobian] = fsolve(@v1_gaslift_model0,x0,opt,u,par);
if flag<0
    flag=0; %Flag ==0 means that x_eq is not good
else
    flag=1; %Flag ==1 means that x_eq is good
end
% disp('Function value at located EQ-point is :');
% x_eq
y = v1_gaslift_model([],x_eq,u,'measurements',par); %Measurments at EQ point
str1 = ['Injected gas into annulus [kg/s]: ' num2str(y(1))];    disp(str1);
str2 = ['Pressure annulus top [bar]: ' num2str(y(2))];          disp(str2);
str3 = ['Pressure riser top [bar]: ' num2str(y(3))];            disp(str3);
str4 = ['Pressure annulus bottom [bar]: ' num2str(y(4))];       disp(str4);
str5 = ['Pressure riser bottom [bar]: ' num2str(y(5))];         disp(str5);
str6 = ['Choke mass flowrate [kg/s]: ' num2str(y(6))]; disp(str6);
str7 = ['Mixture density top [kg/m^3]: ' num2str(y(7))];        disp(str7);
str8 = ['Liquid volume fraction top [--]: ' num2str(y(8))];     disp(str8);

end
	