function Y=UNIQUAC(xA,xB,T,NT)

% Program written by Stathis Skouras, 20/3/2000
%******************************************************************************* 
% Program for ternary VLE by using UNIQUAC equation
%*******************************************************************************
% Ternary system
%                (1): Methanol (light)
%                (2): Ethanol (intermediate)
%                (3): 1-Propanol (heavy)
%*******************************************************************************
% Pressure is 1 atm = 1.013 bar = 760 mmHg
% Composition vectors xA, xB, are line vectors of NT columns
%*******************************************************************************
% Data that have to be changed for different mixture

% Pure component constants for Antoine equation (Gmehling et al, 1977-1990)
A1=8.08097; B1=1582.271; C1=239.726;  % (T-range: -15C - 84C, Component 1)
A2=8.11220; B2=1592.864; C2=226.184;  % (T-range: -20C - 93C, Component 2)
A3=8.37895; B3=1788.020; C3=227.438;  % (T-range: -15C - 98C, Component 3)

% Partial molar volumes of pure components (Gmehling et al, 1977-1990)
V1=40.73;     % Component 1 (cm^3/mol)
V2=58.68;     % Component 2 (cm^3/mol)
V3=75.14;     % Component 3 (cm^3/mol)

% Volume parameters of pure components (Gmehling et al, 1977-1990)
R1=1.4311;     % Component 1 (cm^3/mol)
R2=2.1055;     % Component 2 (cm^3/mol)
R3=2.7799;     % Component 3 (cm^3/mol)

% Area parameters of pure components (Gmehling et al, 1977-1990)
Q1=1.4320;     % Component 1 (cm^3/mol)
Q2=1.9720;     % Component 2 (cm^3/mol)
Q3=2.5120;     % Component 3 (cm^3/mol)

% Binary interaction parameters Aij for UNIQUAC equation (Gmehling et al. 1977-1990)
A11=0;                       % Pure component 1 (cal/mol)
A22=0;                       % Pure component 2 (cal/mol)
A33=0;                       % Pure component 3 (cal/mol)
A12=-474.7791;                % Binary mixture 1-2 (cal/mol) 
A21=762.8153;                  % Binary mixture 1-2 (cal/mol) 
A13=-119.5555;                % Binary mixture 1-3 (cal/mol) 
A31=277.2660;                 % Binary mixture 1-3 (cal/mol) 
A23=-35.2815;                % Binary mixture 2-3 (cal/mol) 
A32=71.7593;                 % Binary mixture 2-3 (cal/mol)

% End of data that have to be changed for different mixture
%**********************************************************************************

% Gas constant
R=1.98721;     % cal/mol K

%Coordination number
z=10;

L1=z/2*(R1-Q1)-(R1-1);
L2=z/2*(R2-Q2)-(R2-1);
L3=z/2*(R3-Q3)-(R3-1);


 % Area fraction of pure components
theta1=Q1*xA./(Q1*xA + Q2*xB + (1-xA-xB)*Q3);
theta2=Q2*xB./(Q1*xA + Q2*xB + (1-xA-xB)*Q3);
theta3=Q3*(1-xA-xB)./(Q1*xA + Q2*xB + (1-xA-xB)*Q3);

% Volume fraction of pure components
phi1=R1*xA./(R1*xA + R2.*xB + (1-xA-xB)*R3);
phi2=R2*xB./(R1*xA + R2*xB + (1-xA-xB)*R3);
phi3=R3*(1-xA-xB)./(R1*xA + R2*xB + (1-xA-xB)*R3);

% Binary interaction parameters for UNIQUAC equation 
T11=exp(-A11./(R*T'));            % Pure component 1 
T22=exp(-A22./(R*T'));            % Pure component 2 
T33=exp(-A33./(R*T'));            % Pure component 3 
T12=exp(-A12./(R*T'));            % Binary mixture 1-2 
T21=exp(-A21./(R*T'));            % Binary mixture 1-2 
T13=exp(-A13./(R*T'));            % Binary mixture 1-3  
T31=exp(-A31./(R*T'));            % Binary mixture 1-3  
T23=exp(-A23./(R*T'));            % Binary mixture 2-3  
T32=exp(-A32./(R*T'));            % Binary mixture 2-3 

%Vapor pressures by Antoine equation
P1s=10.^(A1-B1./(T'-273.15 + C1));                % Component 1 (mmHg)
P2s=10.^(A2-B2./(T'-273.15 + C2));                % Component 2 (mmHg)
P3s=10.^(A3-B3./(T'-273.15 + C3));                % Component 3 (mmHg)

i=1:NT;

% Activity coefficients (combinatorial part) by UNIQUAC equation
G1c(i)=log(phi1(i)./xA(i)) + z/2*Q1*log(theta1(i)./phi1(i)) + L1...
      - (phi1(i)./xA(i)).*(xA(i).*L1+xB(i).*L2+(1-xA(i)-xB(i)).*L3);
  
G2c(i)=log(phi2(i)./xB(i)) + z/2*Q2*log(theta2(i)./phi2(i)) + L2...
      - (phi2(i)./xB(i)).*(xA(i).*L1+xB(i).*L2+(1-xA(i)-xB(i)).*L3);
  
G3c(i)=log(phi3(i)./(1-xA(i)-xB(i))) + z/2*Q3*log(theta3(i)./phi3(i)) + L3...
      - (phi3(i)./(1-xA(i)-xB(i))).*(xA(i).*L1+xB(i).*L2+(1-xA(i)-xB(i)).*L3);

% Activity coefficients (residual part) by UNIQUAC equation
G1r(i) = Q1.*(1-log(theta1(i).*T11(i)+theta2(i).*T21(i)+theta3(i).*T31(i))...
          - theta1(i).*T11(i)./(theta1(i).*T11(i)+theta2(i).*T21(i)+theta3(i).*T31(i))...
          - theta2(i).*T12(i)./(theta1(i).*T12(i)+theta2(i).*T22(i)+theta3(i).*T32(i))...    
          - theta3.*T13(i)./(theta1(i).*T13(i)+theta2(i).*T23(i)+theta3(i).*T33(i)));    % Component 1

G2r(i) = Q2.*(1-log(theta1(i).*T12(i)+theta2(i).*T22(i)+theta3(i).*T32(i))...
          - theta1(i).*T21(i)./(theta1(i).*T11(i)+theta2(i).*T21(i)+theta3(i).*T31(i))...
          - theta2(i).*T22(i)./(theta1(i).*T12(i)+theta2(i).*T22(i)+theta3(i).*T32(i))...    
          - theta3.*T23(i)./(theta1(i).*T13(i)+theta2(i).*T23(i)+theta3(i).*T33(i)));    % Component 2

G3r(i) = Q3.*(1-log(theta1(i).*T13(i)+theta2(i).*T23(i)+theta3(i).*T33(i))...
          - theta1(i).*T31(i)./(theta1(i).*T11(i)+theta2(i).*T21(i)+theta3(i).*T31(i))...
          - theta2(i).*T32(i)./(theta1(i).*T12(i)+theta2(i).*T22(i)+theta3(i).*T32(i))...    
          - theta3.*T33(i)./(theta1(i).*T13(i)+theta2(i).*T23(i)+theta3(i).*T33(i)));    % Component 3

% Activity coefficients by UNIQUAC equation
G1(i)=exp(G1c(i)+G1r(i));                        % Component 1
G2(i)=exp(G2c(i)+G2r(i));                        % Component 2
G3(i)=exp(G3c(i)+G3r(i));                        % Component 3

% Total pressure in each stage (bar)
DP=1-(xA(i).*G1(i).*P1s(i) + xB(i).*G2(i).*P2s(i)+(1-xA(i)-xB(i)).*G3(i).*P3s(i))/760;

% Vapor Liquid Equilibria (VLE)
yA(i)=xA(i).*G1(i).*P1s(i)./(xA(i).*G1(i).*P1s(i) + xB(i).*G2(i).*P2s(i) + (1-xA(i)-xB(i)).*G3(i).*P3s(i));            % Component 1
yB(i)=xB(i).*G2(i).*P2s(i)./(xA(i).*G1(i).*P1s(i) + xB(i).*G2(i).*P2s(i) + (1-xA(i)-xB(i)).*G3(i).*P3s(i));            % Component 2
yC(i)=(1-xA(i)-xB(i)).*G3(i).*P3s(i)./(xA(i).*G1(i).*P1s(i) + xB(i).*G2(i).*P2s(i) + (1-xA(i)-xB(i)).*G3(i).*P3s(i));  % Component 3

Y=[yA;yB;T';DP];