function Y = NRTL(T); 

% This function uses the NRTL model at given values of temperature (T), pressure (P) and liquid composition (X)
% Returns the vector of vapour compositions (Y).

global X NC P R A W alpha 

% Antoine equation where Psat [mmHg] and temperature T [K]
for i=1:NC
  log10Ps = A(i,1) - A(i,2) / (T - 273.15 + A(i,3));
  Ps(i)   = 10^log10Ps;
end

% Starting NRTL model
% Calculate the matrix of NRTL coefficients TAU_ji
TAU=W/(R*T);

% Calculate the matrix of NRTL coefficients G_ji
G=exp(-alpha.*TAU);
 
% Calculate the log of the activity coefficients
Y = G'*X;
Z = (TAU.*G)'*X./Y;
E = ones(size(X));

LNGAMMA = Z + G.*(TAU-E*Z')*(X./Y);     % logarithm of the activity coefficients
GAMMA = exp(LNGAMMA);                   % activity coefficients
% Finishing NRTL model

% VLE formula, P*Y(i)=X(i)*GAMMA(i)*Ps(i) 
Y=Ps'.*X.*GAMMA./(P*760);                % Vector of vapour compositions