function [fl,fg]=A1_fugP(T,P,Tc,Pc,w)
%Computes fugacity of pure component
% input Tc,T in deg K
% input Pc,[P1:incP:p2] in bar 
% input w accentric factros

% Usage 

%%%[fl,fg]=A1_fugP(280,[0.4:.05:0.8],369.83,4.245517,0.152291)

%Note that this function takes a range of pressure as input
% and computes fugacities of liquid and vapor - returns them as a vector

R=8.314472E-3; %MPa l/mol/K

n=length(P);  %Find out number of components

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

for i=1:n
  %Proceed with calculation normally
  A   = a*P(i)/(R*T)^2 ;
  B   = b*P(i)/(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)   ;%guess what is happening here?
        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(i)*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(i)*exp(tempG);  %Calculates vapor fugacity
end
