function Y = WILSON(T); 

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

global X P NC R W V A RP Q

% 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 UNIQUAC model
% Coordination number z
z=10;

% Calculate area fractions THETA
THETA=Q'.*X/(Q*X);

% Calculate volume fractions PHI
PHI=RP'.*X/(RP*X);

% Calculate l
l=z/2*(RP-Q)-(RP-1);

% Calculate the matrix of UNIQUAC coefficients TAF_ij
for i=1:NC
  for j=1:NC
    if i==j
     TAF(j,i)=1;
    else
     TAF(j,i)=exp(-W(j,i)/(R*T));
    end
  end
end

% Calculate combinatorial part of the activity coefficients
lnGammaC=(log(PHI./X) + z/2*Q'.*log(THETA./PHI) + l' - PHI./X*(l*X))';

% Calculate residual part of the activity coefficients
THETAA=THETA';
for i=1:NC
    for j=1:NC
      DENOM(j)=THETAA*TAF(:,j);
      RATIO(j)=TAF(i,j)./DENOM(j);
    end
      lnGammaR(i)=Q(i)*(1-log(THETAA*TAF(:,i))-RATIO*THETA);
end

% Calculate activity coefficients
GAMMA = (exp(lnGammaC + lnGammaR))';
% Finishing UNIQUAC model

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