function [phi,H,Vm,A,B]= srkHEXV3(Z,x,T,P)
% In: compositon x (vector), Temperature T [K], pressure p [MPa]
% OUT: Compressibility Z, Fugacity coeff.(vector) phi [-], 
%enthalpy H [MJ/kmol], Molar Volume V [m3/kmol]

% Authors: Andreas Linhart and Sigurd Skogestad. NTNU. 2008
 CORRECTED phi to include kij correctly 
%(authors: Marie Solvik and  Sigurd Skogestad, NTNU, Nov. 2013)
% Adapted to a mixed refirgerent stream (Cristina Zotica, Feb 2017)
% removed 'liquid', 'vapor' as global variables
% added 5th term of Cp in % calculating H;
% changed P in MPa and H in MJ/kmol

%Components N2 C1 C2 C3 C4 
NC=4;  % Nr of components (1)Nitrogen - (2)Methane - (3)Ethane - (4)Propane 

%Component data from: Reid RC, Prausnitz JM, Poling BE. 
%The properties of gases \& liquids (5th edition)

compData.Pc = [33.98  45.99  48.72  42.48 ]*0.1;     % [MPa]
compData.Tc = [126.20 190.56 305.32 369.83];         % [K]
compData.w  = [0.037  0.011  0.099  0.152 ];         % acentric factor [-]

R = 8.314/1000; %MJ/kmol/K

% Ideal has heat capacity (See Hid below): 
compData.Cp(:, 1) = [3.539 -0.261e-3 0.007e-5   0.157e-8 -2.099e-11]' * R;    
compData.Cp(:, 2) = [4.568 -8.975e-3 3.631e-5  -3.407e-8  1.091e-11]' * R;    
compData.Cp(:, 3) = [4.178 -4.427e-3 5.660e-5  -6.651e-8  2.487e-11]' * R;    
compData.Cp(:, 4) = [3.847  5.131e-3 6.011e-5  -7.893e-8  3.079e-11]' * R;    

compData.Tref = 298.15;     % [K], for Ideal gas heat capacity
compData.Pref = 0.1;        % [MPa]   

% SRK binary interaction parameters from Aspen Hysys V9
kinteraction(:,1) = [0 3.11990007758141e-002 3.18990014493465e-002...
                     8.86000022292137e-002];     %Nitrogen
kinteraction(:,2) = [ 3.11990007758141e-002	0	2.24137306213379e-003...
                    6.82878494262695e-003]; %Methane
kinteraction(:,3) = [3.18990014493465e-002	2.24137306213379e-003	0 ...
                    1.25795602798462e-003]; %Ethane
kinteraction(:,4) = [8.86000022292137e-002	6.82878494262695e-003...
                     1.25795602798462e-003	0]; %Propane
 %kinteraction = zeros(4,4);

%% Initialize

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

% 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: 
%Mixing rules
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

%phi
corrphi=zeros(1,NC);
for i = 1:NC
    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));

%end

% 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

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


for i = 1:NC
    H = H+(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) + 1/5 * Cp(5,i) * (T^5 - Tref^5)));
end


% 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
  Vm = ((Z * R * T / P)- c);
else  %vapor
    Vm = Z * R * T / P;
end

