%% SRK EO Package
function [phi,H,V,errceos]= srkseo(x,T,P,Z)
% In: compositon x (vector), Temperature T [K], pressure p [N/m2], Compressibility Z
% OUT: Fugacity coeff.(vector) phi [-], enthalpy H [J/mol], Molar Volume
% V[m3/mol], Error in the ceos errceos. 

%Based on Linhart & Skogestad, 2008

%Components N2 C1 C2 C3 C4 

NC=5;  % Number of components (1)Nitrogen - (2)Methane - (3)Ethane - (4)Propane - (5)Butane 

%Component data from: Reid RC, Prausnitz JM, Poling BE. The properties of gases \& liquids (4th edition). New York: McGraw-Hill, Inc., 1987.
compData.Pc=[33.5 45.4 48.2 41.9 37.5 ]*1.e5;    % [N/m2]
compData.Tc=[126.2 190.6 305.4 369.8 425.2 ];  % [K]
compData.w=[0.040 0.008 0.098 0.152 0.193 ];   % acentric factor [-]
% Ideal has heat capacity (See Hid below): 
compData.Cp(:, 1) = [7.440 -0.324e-2 6.400e-6 -2.790e-9]'; %Nitrogen
compData.Cp(:, 2) = [4.598 1.245e-2 2.86e-6 -2.703e-9]'; % Methane
compData.Cp(:, 3) = [1.292 4.254e-2 -1.657e-5 2.081e-9]'; % Ethane
compData.Cp(:, 4) = [-1.009 7.315e-2 -3.789e-5 7.678e-9]';    % Propane
compData.Cp(:, 5) = [2.266 7.913e-2 -2.647e-5 -0.674e-9]';    % n-butane
compData.Tref = 298.15;  % for Ideal gas heat capacity

kinteraction=zeros(NC,NC);      % Here: SRK binary interaction parameters set to zero

% initialize

R = 8.314; %kJ/kmolK

Pc=compData.Pc;
Tc=compData.Tc;
w=compData.w;
ZRA=0.29056-0.08775*w;
Cp=compData.Cp;
Tref=compData.Tref;

% Calculations for given T, P  and composition (x)

%All equations from Soave,1972
Tre=T./Tc;
Pre=P./Pc;
m=0.480+1.574.*w-0.176.*w.^2; %Eq 15
a=(1+m.*(1-Tre.^0.5)).^2; %Eq 13
Ap=0.42747.*a.*Pre./Tre.^2; %Eq 8
Bp=0.08664.*Pre./Tre; % Eq 9

% Start calculations
% Binary a's:
Ab=zeros(NC,NC);
for i=1:NC
    for j=1:NC
        Ab(i,j)=(Ap(i)*Ap(j))^0.5;
    end
end
% Mixture a and b
A=0;
for i=1:NC
    for j=1:NC
        A=A+x(i)*x(j)*Ab(i,j)*(1-kinteraction(i,j));
    end
end
B=0;
for i=1:NC
    B=B+x(i)*Bp(i);
end

% Compute error for cubic equation
C3=1; C2=-1;  C1=A-B-B^2; C0=-A*B;
errceos = C3*Z^3 + C2*Z^2 + C1*Z + C0;

% Fugacity coefficient
% Correction added Nov. 2013:
for i = 1:NC
   corrphi(i) = 0;
for j = 1:NC
   corrphi(i) = corrphi(i) + x(j) * Ap(j)^0.5 * (1-kinteraction(i,j)); 
end 
end
phi=exp((Z-1).*Bp/B-log(Z-B)...
     -A/B*log((Z+B)/Z).*((corrphi.*2.*Ap.^0.5./A)-Bp./B)); %Eq 21 Soave,1972

% Enthalpy Molar 
dadT=0;
for i = 1:NC
    for j = 1:NC
        dadT = dadT -R / 2 * sqrt(abs(0.42748 / T)) * x(i) * x(j) * (m(j) * sqrt(abs(Ap(i)*(Tc(j)/P*(T^2)/(Pc(j))*(R^2)))) + m(i) * sqrt(abs(Ap(j) * (Tc(i)/P*(T^2)/(Pc(i))*(R^2))))) ;
    end
end

Hsrk=(1/ (B*R*T/P)) * ( A*R^2*T^2/P - T*dadT) * log(Z / (Z + B)) + R * T * (Z - 1) ;

Hid=0;
for i = 1:NC
    Hid = Hid+(x(i)* (Cp(1,i) * (T - Tref) + 1/2 * Cp(2,i) * (T^2 - Tref^2) + 1/3 * Cp(3,i) * (T^3 - Tref^3) + 1/4 * Cp(4,i) * (T^4 - Tref^4)));
end
%Conversion
Hid=Hid*4.184; %From kcal to kJ
H=Hid+Hsrk;
Hsrk;
Hid;


% Density (more precisely: molar volume)

if Z<0.1 % Correct liquid SRK-volume using Peneleoux correction
     c=0;
     for i=1:NC
     c=c+x(i) * (0.40768 * (0.29441 - ZRA(i)) * (R * Tc(i)) / (Pc(i))) ;
     end
   V = ((Z * R * T / P)- c);
else  % vapor
   V = Z * R * T / P;
end