function y = wilson(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 Wilson equation.
 
global X P NC R W V A

% Calculate the matrix of Wilson coefficients LAMBDA_ij

for i=1:NC
  for j=1:NC
    if i==j
     LAMBDA(i,j)=1;
    else
     LAMBDA(i,j)=(V(j)/V(i))*exp(-W(i,j)/(R*Tk));
    end
  end
end

for i=1:NC
  for k=1:NC
    RATIO(k)=LAMBDA(k,i)/(X'*LAMBDA(k,:)');
  end
  z(i)=-log(X'*LAMBDA(i,:)')+1-X'*RATIO';
end
  
y=exp(z');
