function Y = UNIQUAC2(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 l XU

% Necessary data:
R = 1.98721;    % universal gas constant [cal/molK]
NC = 6;         % number of components
P = 1;          % pressure [atm]

% UNIQUAC Parameters:
%      TG       DG       MG       BD      GL        MetOH
W = [   0     -48.519   213.876  47.076   836.044  3030.060
      59.273     0      -79.159  26.437   757.300  -89.061
     -139.052  114.633    0      219.802  317.795  -113.482
     -45.259   -16.759  -125.764   0      823.302     0
     -266.224  -274.223 -60.390  -262.612    0     411.932
     -704.965  870.689  534.174    0      -267.337    0    ];
 
 % Coefficients for the vapour pressure equations:
%      A            B           C          D          E
A = [505.205  -9.41262*10^4  -59.7879  5.8959*10^-19  6     %TG
     18.1551  -7.16544*10^3  -0.847809 3.18337*10^-20 6     %DG
     145.254  -2.11542*10^4  -16.8740  1.74025*10^-18 6     %MG
     107.534  -1.56037*10^4  -12.2136  1.79996*10^-18 6     %BD
     169.211  -1.68937*10^4  -21.8125  1.14067*10^-5  2     %GL
     59.8373  -6.28289*10^3  -6.37873  4.61746*10^-6  2];   %MetOH
 % Volume parameters of pure components:
%       TG       DG     MG      BD    GL    MetOH
RP = [39.6176 28.0103 16.403 13.5084 4.7957 1.4311];

% Area parameters of pure components:
%      TG     DG    MG    BD    GL   MetOH
Q = [32.094 23.032 13.97 11.11 4.908 1.432];

% Vapour pressure equation where Psat [mmHg] and temperature T [K]
Ps(1,1) = 760*9.869*(10^-3)*exp(A(1,1) + (A(1,2)/T) + (A(1,3)*log(T)) + (A(1,4)*T^A(1,5)));   %TG
Ps(1,2) = 760*9.869*(10^-3)*exp(A(2,1) + (A(2,2)/T) + (A(2,3)*log(T)) + (A(2,4)*T^A(2,5)));   %DG
Ps(1,3) = 760*9.869*(10^-3)*exp(A(3,1) + (A(3,2)/T) + (A(3,3)*log(T)) + (A(3,4)*T^A(3,5)));   %MG
Ps(1,4) = 760*9.869*(10^-3)*exp(A(4,1) + (A(4,2)/T) + (A(4,3)*log(T)) + (A(4,4)*T^A(4,5)));   %BD
Ps(1,5) = 760*9.869*(10^-3)*exp(A(5,1) + (A(5,2)/T) + (A(5,3)*log(T)) + (A(5,4)*T^A(5,5)));   %GL
Ps(1,6) = 760*9.869*(10^-3)*exp(A(6,1) + (A(6,2)/T) + (A(6,3)*log(T)) + (A(6,4)*T^A(6,5)));   %MetOH

% Setting the coordination number
z = 10;

% Calculating the surface area fractions:
F = Q./sum(Q.*XU(l,:));

% Calculating the volume fraction per mixture mole fraction:
PHI = RP./sum(RP.*XU(l,:));

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

% Calculating the combinatorial contribution:
lnGammaC = ones(1,NC);
for i = 1:NC
    lnGammaC(1,i) = (1-PHI(1,i)+log(PHI(1,i)))-((z/2)*Q(1,i)*(1-(PHI(1,i)/F(1,i))+log(PHI(1,i)/F(1,i))));
end


% Calculating the residual part
NOM = ones(1,NC); RATIO = ones(1,NC); lnGammaR = ones(1,NC);

for i = 1:NC
        NOM(1,i) = sum(Q.*XU(l,:).*TAF(i,:));
        DENOM = sum(Q.*XU(l,:));
        RATIO(1,i) = sum((Q.*XU(l,:).*TAF(i,:))./NOM(1,i));
        lnGammaR(1,i) = Q(1,i)*(1-log(NOM(1,i)/DENOM)-RATIO(1,i));
end


% Calculating the activity coefficient:
lnGamma = lnGammaC + lnGammaR;
Gamma = exp(lnGamma);

Y = (Ps.*XU(l,:).*Gamma)./(P*760);

