function [W_L_choke,W_G_choke]=Sachdeva_choke(P_1,P_2,rho_G1,alpha_L,z,data)
%   Sachdeva Choke model, implicitt version
%   Inputs:     yc-Critical pressure ratio(from previous timestep/iteration)
%               P_1-Upstream pressure [Pa]
%               P_2-Downstream pressure [Pa]
%               rho_G1- Upstream gas density [kg/m^3]
%               alpha_L-Upstream liquid volume fraction
%   Outputs:    W_L_choke-mass flow of liquid through choke [kg/s]
%               W_G_choke-mass flow of gas through choke [kg/s]
%               yc_dev-error in critical pres. ratio calculation
%
%   Written by Espen Storkaas, Jan. 15th 2003
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1=(1-alpha_L)*rho_G1/(alpha_L*data.rho_L+(1-alpha_L)*rho_G1);
y=P_2/P_1;  %Actual pressure ratio
%y=(y>yc)*y+(y<=yc)*yc; %Corrected for possible critical flow
%rho_m=(x1/(rho_G1*y^(1/data.k))+(1-x1)/data.rho_L)^-1;
rho_m=alpha_L*data.rho_L+(1-alpha_L)*rho_G1;
W_choke=(y<1)*1.4152*z*data.Ac*data.Cd*sqrt(P_1*rho_m^2*...
    ((1-x1)*(1-y)/data.rho_L+x1*data.k*(1-y^((data.k-1)/data.k))/((data.k-1)*rho_G1)));%Rett opp konstanten!!!
W_L_choke=(x1<0.9)*(1-x1)*W_choke;
W_G_choke=(x1>0.001)*x1*W_choke;
% mdot_mix=(y<1)*data.K1*z*sqrt(P2*rho_T^2*...
%     ((1-x1)*(1-y)/data.rho_L+x1*data.k*(1-y^((data.k-1)/data.k))/((data.k-1)*rho_G2)));

% n=1+x1*(data.C_pG-data.C_vG)/(x1*data.C_vG+(1-x1)*data.C_L);
% N=data.k/(data.k-1)+(1-x1)*rho_G1*(1-yc)/(x1*data.rho_L);
% D=data.k/(data.k-1)+n/2+n*(1-x1)*rho_G1*yc^(1/data.k)/(x1*data.rho_L)+...
%     n/2*((1-x1)*rho_G1*yc^(1/data.k)/(x1*data.rho_L))^2;
% yc_dev=(N/D)^(data.k/(data.k-1))-yc;
% [y yc]
