function [F_wL,F_wG,F_i]=friction(alpha_L,rho_G,u_L,u_G,data)
%   Functions for computing frictional forces on the fluids
%
%   Written by Espen Storkaas, Jan 16th 2003

%   Averaging mass data and picking out internal velocities
alpha_L2=mean([alpha_L(data.I_i) alpha_L(data.I_i+1)],2);
rho_G=mean([rho_G(data.I_i) rho_G(data.I_i+1)],2);
u_L=u_L(data.I_i+1);
u_G=u_G(data.I_i+1);
%   Classification of flow regime
I_Lam=+(data.theta<pi/20);
I_Trans=(data.theta>=pi/20)&(alpha_L2>data.Ann_lim)&(alpha_L2<=data.Bub_lim);
I_Trans_rel=I_Trans.*sin((0.5/(data.Bub_lim-data.Ann_lim))*pi*(alpha_L2-data.Ann_lim));
I_Bub=+((data.theta>=pi/20)&(alpha_L2>data.Bub_lim))+I_Trans_rel;
I_Ann=+((data.theta>=pi/20)&(alpha_L2<=data.Ann_lim))+(+I_Trans).*(1-I_Trans_rel);

S_L=pi*(I_Lam.*alpha_L2+abs(I_Lam-1))*data.D;
S_G=pi*(I_Lam-I_Lam.*alpha_L2)*data.D;
S_i=I_Lam.*(alpha_L2.^2-alpha_L2)*(-4*data.D)+...
    I_Ann.*pi*data.D.*sqrt(1-alpha_L2)+...
    I_Bub.*pi.*(1-alpha_L2)*data.D^2/data.Db;

D_hL=4*alpha_L2*data.A./S_L;
D_hG=4*(1-alpha_L2)*data.A./(S_G+S_i);
Re_L=abs(u_L*data.rho_L.*D_hL/data.my_L);
Re_G=abs(u_G.*rho_G.*D_hG/data.my_G);
lambda_L=data.Lparam*(Re_L~=0).*max(64./Re_L,0.0055*(1+(2e4*data.eps./D_hL+1e6./Re_L).^(1/3)));
lambda_G=data.Gparam*(Re_G~=0).*max(64./Re_G,0.0055*(1+(2e4*data.eps./D_hG+1e6./Re_G).^(1/3)));
lambda_i=data.iparam*0.02*(1+75*alpha_L2)/4;
[d_L,d_G]=flowdir(u_L,u_G);%d_L(end)=0;d_G(end)=0;
F_xL=500000*(u_L).*exp(5-500*alpha_L(data.I_i+d_L));
F_xG=500000*(u_G).*exp(5-500*(1-alpha_L(data.I_i+d_G)));

F_wL=S_L/data.A.*lambda_L/8*data.rho_L.*abs(u_L).*u_L+F_xL;
F_wG=S_G/data.A.*lambda_G/8.*rho_G.*abs(u_G).*u_G+F_xG;
F_i=S_i/data.A.*lambda_i/8.*rho_G.*abs(u_G-u_L).*(u_G-u_L);
