function par = v1_new_6d_parameters()
% ********** Constants and parameters for New 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.theta= pi/180;                  % Feed pipe inclination (Rad) = 1 Deg
par.r_w = 0.06;                      % Raduis of well tubing (m)
par.r_p = 0.06;                      % Raduis of pipe (m)
par.r_r = 0.05;                      % Raduis of riser (m)
par.hc = 2*par.r_w/cos(par.theta);   % Critical liquid level (m)
par.Aw = pi*par.r_w^2;               % Cross section area of well tubing (m2)
par.Ap = pi*par.r_p^2;               % Cross section area of pipeline (m2)
par.Ar = pi*par.r_r^2;               % Cross section area of riser (m2)
par.Lw = 3000;                      % Height of oil well (m)
par.Lp = 4300;                      % Length of upstream pipe (m)
par.Lr = 300;                       % Height of riser
par.Lh = 100;                       % length of horozontal top section (m)
par.Tw = 369;                       % Well temprature (K)
par.Tp = 337;                       % Pipeline temprature (K)
%par.T1 = 300;
par.Tr = 298.3;                     % Riser temprature (K)
%par.T2=300;
%par.M_G = 23;                       % Molecular weight of Gas (kg/kmol)
par.M_G = 20;                    % Molecular weight of Gas (kg/kmol)
par.Pr = 320e5;                    % Reservior Pressure
par.P_s = 50.1e5;                    % Pressure after choke valve (Pa)

par.Vw = par.Aw*par.Lw;             % Total volume in upstream (m3)
par.Vp = par.Ap*par.Lp;    % Total volume in riser (m3)
par.Vr = par.Ar*(par.Lr+par.Lh);             % Total volume of well (m3)
par.PI = 2.75e-6;                     % Productivity Index of the well
% *************************************************************************
par.Pnormal = 101325;              % atmospheric conditions for Usg
Tnormal = 293.15;                  % Temprature in Normal conditions
RTnormal =par.R*Tnormal/par.M_G;   % compressibility in Normal Conditions
par.Rog_n = par.Pnormal/RTnormal;  % Gas density in normal conditions
par.visl= 1.4260e-4;
par.visg= 1.39e-5;
par.eps = 2.8e-5;                  % Roughness of pipe
% *************** Finding Orifice Coefficients ****************************
% ******** Steady State Values Baesd on OLGA *************
z1 = 1;
z2 = 0.04;
P_bh = 288.552e5;
w_r = par.PI*max(par.Pr - P_bh,0);
w_Gr = 0;
P_wh = 77.3583e5;
Alpha_Gwhm = 0.0395062;
w_Gwh = Alpha_Gwhm*w_r;
% -------------------------
w_in = 8.64813;
Alpha_Ginm = 0.0395225;
w_Gin = w_in*Alpha_Ginm;
w_Lin = w_in - w_Gin;
P_p = 77.2701e5;
%Alpha_L2_av = 0.717;
P_rt = 58.5413e5;
% ********************************************************
Alpha_Gwmt = w_Gwh/w_r;

w_Gout = w_Gin;
w_Lout = w_Lin;
w_mix = w_Gout + w_Lout;

%par.a = (Alpha_Gwmt-Alpha_Gwmb)/(P_bh - P_wh);
%par.b = Alpha_Gwmb + par.a*P_bh;

rho_Gw = P_wh*par.M_G/(par.R*par.Tp);
rho_mixwt = Alpha_Gwmt*rho_Gw + (1-Alpha_Gwmt)*par.rho_L;
par.K1 = 1*w_mix/(z1*sqrt(rho_mixwt*max(0,P_wh-P_p)));

rho_Gp = P_p*par.M_G/(par.R*par.Tp);
par.Alpha_Lp_av = w_Lin*rho_Gp/(w_Lin*rho_Gp + w_Gin*par.rho_L);

par.w_nom = 12;%w_in;
par.k_a = 1.94;
par.k_h = 0.6;
hss = par.k_h*par.Alpha_Lp_av*par.hc;
h = hss;

rho_mix_p = par.Alpha_Lp_av*par.rho_L + (1-par.Alpha_Lp_av)*rho_Gp;

Uslp= w_Lin/(par.Ap*par.rho_L);
Usgp = w_Lin/(par.Ap*rho_Gp);

Umin = Uslp + Usgp;

vism_p = par.Alpha_Lp_av*par.visl + (1-par.Alpha_Lp_av)*par.visg;
Re_p=rho_mix_p*Umin*(2*par.r_p)/vism_p;

% Fric_p=0.5*Lambda_p*par.Alpha_Lp_av*par.rho_L*Uslp^2*par.Lp/(2*par.r_p);
% temp = -1.8*log10((par.eps/(2*par.r_p)/3.7)^1.11+6.9/Re_p);
% Lambda_p=(1/temp)^2;
Lambda_p=0.0056 + 0.5*Re_p^(-0.32);


Fric_p= 0.5*Lambda_p*par.rho_L*Uslp^2*par.Lp/(2*par.r_p);


Pr_av = (P_p+par.rho_L*par.g*h-Fric_p+P_rt)/2;
rho_Gr_av = Pr_av*par.M_G/(par.R*par.Tr);
Alpha_Lr_av = rho_Gr_av*w_Lin/(rho_Gr_av*w_Lin+par.rho_L*w_Gin);
rho_mixr_av = rho_Gr_av*(1-Alpha_Lr_av)+par.rho_L*Alpha_Lr_av;
rho_Gr = P_rt*par.M_G/(par.R*par.Tr);

Uslr = w_Lin/(par.rho_L*par.Ar);
Usgr = w_Gin/(rho_Gr*par.Ar);
Umr = Uslr+Usgr;

vism_r = Alpha_Lr_av*par.visl + (1-Alpha_Lr_av)*par.visg;
Re_r=rho_mixr_av*Umr*(2*par.r_r)/vism_r;
%Lambda_r=0.0056 + 0.5*Re_r^(-0.32);
%Fric_r=0.5*Lambda_r*Alpha_Lr_av*rho_mixr_av*Umr^2*(par.Lr+par.Lh)/(2*par.r_r);

temp = -1.8*log10((par.eps/(2*par.r_r)/3.7)^1.11+6.9/Re_r);
Lambda2=(1/temp)^2;
%Alpha_Lr_av*
Fric_r= Lambda2*rho_mixr_av*Umr^2*(par.Lr+par.Lh)/(2*par.r_r);


P_rb = P_rt + rho_mixr_av*par.g*par.Lr + Fric_r;


Alpha_Lrt = rho_Gr*w_Lin/(rho_Gr*w_Lin+par.rho_L*w_Gin);

rho_rt = Alpha_Lrt*par.rho_L + (1-Alpha_Lrt)*rho_Gr;

A_g = (h<par.hc)*(par.Ap/par.hc^2)*(par.hc - h)^2;

A_l = par.Ap - A_g;

v_G_lp = w_Gin/(rho_Gp*A_g);

par.K_g = 0.85*v_G_lp/sqrt((P_p-Fric_p-P_rb)/rho_Gr);

par.K_o = 3.6*w_Lin/(A_l*sqrt(par.rho_L*(P_p-Fric_p+par.rho_L*par.g*h-P_rb)));

par.K2 = 1.2*w_mix/(z2*sqrt(rho_rt*(P_rt-par.P_s)));


end