function y = NRTL(Tk); 

% This function returns the activity coefficient of each component in 
% an NC component mixture of liquid composition X and temperature Tk [K] 
% as calculated by the NRTL equation.
 
global X P NC R W A alpha 

% Calculate the matrix of NRTL coefficients TAF_ji
TAF=W/(R*Tk);

% Calculate the matrix of NRTL coefficients G_ji
G=exp(-alpha.*TAF);

% Calculate the log of the activity coefficients
Y = G'*X;
Z = (TAF.*G)'*X./Y;
E = ones(size(X));
LNGAMMA = Z + G.*(TAF-E*Z')*(X./Y);

% Calculate activity coefficients
GAMMA = exp(LNGAMMA);

y=GAMMA;
