function f=fun(xx)
%Main Program that assembles all equations for A2 P2
%fugacities from Peng-Robinson equation
%Components - propane, n-butane

Tc  =[369.83 425.12]; %K
Pc  =[41.9 37.46]*1.01325; %bar
w   =[0.152291 0.200164];
k(2,2) = 0;
k(1,2) = 8.19e-004; %propane, n-butane binary interaction for PR EOS
k(2,1) = k(1,2);
P = 1.0135;  %bar (14.7 psia)

%Transfer the guesses to composition vectors and T
x3 = xx(1:2);
y2 = xx(3:4);
T2 = xx(5);
x4 = xx(6:7);
y5 = xx(8:9);
T4 = xx(10);
x = [x3 x4];
y = [y2 y5];
T = [T2 T4];
% initialize variables
fl(2,2) = 0;
fv(2,2) = 0;
f(10) = 0;

%Compute fugacity of each species in both phases in both stages
for j = 1:2,    % stage #
    for i=1:2,  % species 
        fl(i,j) = fugL(T(j),P,Tc,Pc,w,k,x(:,j),i);
        fv(i,j) = fugV(T(j),P,Tc,Pc,w,k,y(:,j),i);
        f(2*(j-1)+i) = fl(i,j) - fv(i,j);
    end
end

%Summation constraints
f(5) = sum(x3) - 1;
f(6) = sum(y2) - 1;
f(7) = sum(x4) - 1;
f(8) = sum(y5) - 1;

% Propane material balance
f(9) = 0.6*x3(1) + 0.4*y2(1) - 0.4;
f(10) = 0.6*x4(1) + 0.4*y5(1) - y2(1);

%--------------------------------------------------------
function f=fugL(T,P,Tc,Pc,w,k,x,id)
%Computes fugacity of i (liquid phase) in a mix of comp x
% Tc,T in K
% Pc,P in bar
% x in mole fraction
% k matrix of interaction coefficients
% w accentric factros
% id component for which fugacity is needed
% Usage f=fugL(T,P,Tc,Pc,w,k,x,id)

R   = 83.14 ;   %bar.cm3/mol/K
n=length(Tc);  %Find out number of components

for i=1:n
    m(i)   = 0.37464 + 1.54226*w(i) - 0.26922*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;
for i=1:n
    for j=1:n
        aij(i,j) = (1-k(i,j))*sqrt(a(i)*a(j));
        am = am + x(i)*x(j)*aij(i,j);
    end
end
bm = 0;
for i=1:n
    bm = bm + x(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
zL=min(zz); %need only liquid fugacities

sum=0;
for i=1:n
    sum = sum + x(i)*aij(i,id);
end
tt = 2*sum/am - b(id)/bm;
tempL = b(id)*(zL-1)/bm - log(zL-B) - (A/2/sqrt(2)/B)* ...
    tt * log((zL+2.414*B)/(zL-0.414*B));
f=P*x(id)*exp(tempL);  %Calculates liquid fugacity


%--------------------------------------------------------
function f=fugV(T,P,Tc,Pc,w,k,y,id)
%Computes fugacity of i (vapour phase) in a mix of comp y
% Tc,T in K
% Pc,P in bar
% y in mole fraction
% k matrix of interaction coefficients
% w accentric factros
% id component for which fugacity is needed
% Usage f=fugV(T,P,Tc,Pc,w,k,y,id)

R   = 83.14 ;   %bar.cm3/mol/K
n=length(Tc);  %Find out number of components

for i=1:n
    m(i)   = 0.37464 + 1.54226*w(i) - 0.26922*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;
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;
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
zV=max(zz); %need only vapor fugacities

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