function [fl,fg,zL,zG]=A1_fugT(T,P,Tc,Pc,w)
%Computes fugacity of pure component
% Tc,T in deg K
% Pc,P in MPa
% w accentric factors

% Usage

%%%[fl,fg]=A1_fugT([260:5:300],0.7,369.83,4.245517,0.152291)

R=8.314472E-3; %MPa l/mol/K
n=length(T);  %Find out number of components

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

%Proceed with calculation normally
A   = a*P/(R*T(i))^2 ;
B   = b*P/(R*T(i));
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   
zL=min(zz); zG=max(zz);

tempL = (zL-1) - log(zL-B) - (A/2/sqrt(2)/B)* ...
        log((zL+2.414*B)/(zL-0.414*B));
     fl(i)=P*exp(tempL);  %Calculates liquid fugacity
     
tempG = (zG-1) - log(zG-B) - (A/2/sqrt(2)/B)* ...
        log((zG+2.414*B)/(zG-0.414*B));
     fg(i)=P*exp(tempG);  %Calculates liquid fugacity
end


