function [Z,phi,H,V]= srksis(x,T,P,phase)
% In: compositon x (vector), Temperature T [K], pressure p [N/m2], phase ['liquid' or 'vapor']
% OUT: Compressibility Z, Fugacity coeff.(vector) phi [-], enthalpy H [J/mol], Molar Volume V [m3/mol]

% Authors: Andreas Linhart and Sigurd Skogestad. NTNU. 2008

%thermodynamic/SRK data
%Component data from: Reid RC, Prausnitz JM, Poling BE. The properties of gases \& liquids (4th edition). New York: McGraw-Hill, Inc., 1987.

NC=4; % propane (1) - isobutane (2) - nbutane (3) - 2methylbutane (4)
compData.Pc=[42.5 36.5 38 33.9]*1.e5;    % [N/m2]
compData.Tc=[369.8 408.2 425.18 460.4];  % [K]
compData.w=[0.153 0.183 0.2008 0.227];   % acentric factor [-]
% Ideal has heat capacity (See Hid below): 
compData.Cp(:, 1) = [-4.224 3.063e-1 -1.586e-4 3.215e-8]';    % propane
compData.Cp(:, 2) = [-1.39, 0.3847, -0.0001846, 2.895E-008]'; % iso-butane
compData.Cp(:, 3) = [9.487 3.313e-1 -1.108e-4 -2.822e-9]';    % n-butane
compData.Cp(:, 4) = [-9.525 5.066e-1 -2.729e-4 5.723e-8 ]';   % 2-methylbutane
compData.Tref =298.15;  % for Ideal gas heat capacity

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

% initialize
liquid=1;vapor=2;
R = 8.314;

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)

Tre=T./Tc;
Pre=P./Pc;
m=0.480+1.574.*w-0.176.*w.^2;
a=(1+m.*(1-Tre.^0.5)).^2;
Ap=0.42747.*a.*Pre./Tre.^2;
Bp=0.08664.*Pre./Tre;

% 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

% Solve cubic equation to find compressibility Z = P*V/(R*T)
Zall=roots([1 -1 A-B-B^2 -A*B]);  
% use real roots only
        Zreal=[];
	for i=1:3 
	  if isreal(Zall(i))==1 
	    Zreal=[Zreal Zall(i)]; 
	  end 
	end
% Select right root depending on phase
if     phase==liquid  Z=min(Zreal);
elseif phase==vapor   Z=max(Zreal);
else   write('error in specifying phase')
end

% Fugacity coefficient
phi=exp((Z-1).*Bp/B-log(Z-B)...
    -A/B*log((Z+B)/Z).*(2.*Ap.^0.5./A^0.5-Bp./B));

% Enthalpy
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

H=Hsrk+Hid;

% Density (more precisely: molar volume)

if phase==liquid % 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


