function [fug]=A1_fug_mix(T,P,y,Tc,Pc,w,k)
%Computes fugacity of pure component
% Tc,T in deg K
% Pc,P in bar
% y - concentration
% w - accentric factors
% k - binary interaction parameters
% Usage fug=A1_fug_mix(T,P,y,Tc,Pc,w,k)

R   = 83.14 ;   %bar.cm3/mol/K
n = length(y) ;%find out the number of components in the mixture

for i=1:n
    m(i)   = 0.37464 + 1.54226*w(i)- 0.26992*w(i)^2  ;
    alp(i) =(1 + m(i)*(1-sqrt(T/Tc(i)) ))^2 ;
    a(i)   = 0.45724 * (R*Tc(i))^2 *alp(i) / Pc(i) ;
    b(i)   = 0.07780*R*Tc(i)/Pc(i) ;
end

%apply mixing rules
am = 0; % a_mixed
for i=1:n
    for j=1:n
        aij(i,j) = (1-k(i,j))*sqrt(a(i)*a(j));
        am = am +y(i)*y(j)*aij(i,j);
    end
end

bm = 0; % b_mixed
for i=1:n
    bm = bm + y(i)*b(i);
end

%Proceed with calculation normally
A   = am*P/(R*T)^2 ;
B   = bm*P/(R*T);
c(1) = 1 ;
c(2) = -(1-B) ;
c(3) = (A-3*B^2-2*B);
c(4) = -(A*B - B^2 - B^3) ;
zz =[];
z = roots(c);
for j=1:3
    if (angle(z(j)) == 0)
        zz = [zz, z(j)];
    end
end
zG=max(zz); %need only vapor fugacities

for id=1:n
    sum = 0;
    for i=1:n
        sum = sum + y(i)*aij(i,id);
    end
    tt = 2*sum/am - b(id)/bm;
    tempG = b(id)*(zG-1)/bm - log(zG-B) - (A/2/sqrt(2)/B)* ...
        tt * log((zG+2.414*B)/(zG-0.414*B));
    fug(id)=P*y(id)*exp(tempG);  %Calculates gas-phase fugacity
end
