function par = v1_DiMeglio_parameters()
clc
% ********** Constants and parameters for Storkaas simple model ***********
par.g = 9.81;                   % Gravity (m/s2)
par.R = 8314;                   % Gas constant (J/(K.Kmol))
par.rho_L = 832.2;              % Liquid density (Kg/m3)
par.theta1= pi/180;             % Feed pipe inclination (Rad) = 1 Deg
par.theta2= pi/2;               % Mean inclination of the riser (Rad) = 90 Deg
par.r1 = 0.06;                  % Raduis of pipe (m)
par.r2 = 0.05;                  % Raduis of riser (m)
par.A1 = pi*par.r1^2;           % Cross section area of pipeline (m2)
par.A2 = pi*par.r2^2;           % Cross section area of riser (m2)
par.L1 = 4300;                  % Length of upstream pipe (m)
par.L2 = 300;                   % Height of riser
par.L3 = 100;                   % length of horozontal top section (m)
par.T1 = 335;                   % Pipeline temprature (K)
par.T2 = 298.3;                 % Riser temprature (K)
%par.M_G = 21.43;                % Molecular weight of Gas (kg/kmol)
par.M_G = 23;                % Molecular weight of Gas (kg/kmol)
par.P0 = 50.1e5;                % Pressure after choke valve (Pa)
par.V2 = par.A2*(par.L2+par.L3);% Total volume in riser (m3)
par.Vr = par.A2*par.L2;         % Total volume in riser (m3)
par.V1 = par.A1*par.L1;         % Total volume in upstream (m3)
par.ep = 0.35;                   % Tuning parameter
F_G1 = 0.2715;                     % Average liquid volume fraction in pipe
par.Veb = F_G1*par.V1;      % Volume of gas in pipe
% *************************************************************************
% *************** Finding choke Coefficients ****************************
% ******** Values Baesd on OLGA *************
z0= 0.04;
P1 = 76.59951e5;
P2_t = 57.81898e5;
P2_tt = 58e5;
P2_b = 75.28658e5;
%P2_t_u = 54.99592e5;
y_b_max = 74.57604e5;
y_t = 50.3294e5;
wG_in = 0.36;
wL_in = 8.64;
% **************************************

par.C_c = 1*(wL_in+wG_in)/(z0*sqrt(P2_t-par.P0));


x = wG_in/wL_in;
b = par.rho_L*par.R*par.T2/par.M_G;
par.m_lstill = 0;
% par.m_lstill = (1+P2_tt/(b*x))*(P2_b - P2_tt)/(par.g*sin(par.theta2)/par.A2)...
%      - P2_tt*par.Vr*par.rho_L/(b*x);
% 
% m_L2 = P2_t*(par.V2*par.rho_L - par.m_lstill)/(P2_t + b*x);
% %m_G2 = x*m_L2;
% 
% %P2b =  P2_t + par.g*sin(par.theta2)*(m_L2+par.m_lstill)/par.A2;
% %P2b =  P2_t + (P2_t*par.V2*par.rho_L + b*x*par.m_lstill)/(P2_t+b*x)*par.g*sin(par.theta2)/par.A2;
% 
% m_ep = (y_b_max-y_t)*par.A2/(par.g*sin(par.theta2));
% par.ep = (y_t/(b*x))*(par.rho_L*par.Vr-m_ep)/(m_ep-par.m_lstill);
% par.C_g = 3*(1-par.ep)*wG_in/(P1-P2_b);
par.C_g = 5.3468e-6;
% opt = optimset('Display','off','TolFun',1e-200,'TolX',1e-200,'MaxIter',1e4,'MaxFunEvals',1e10);
% [v,fval] = fsolve(@lambda_fun,0.30*par.V1,opt,par);
% par.F_G1 = v/par.V1;
% par.Veb = v;

end