function Y=NRTL(xA,xB,T,NT)
%*******************************************************************************
% Program written by Stathis Skouras, April 2003
%******************************************************************************* 
% Program for ternary VLE by using NRTL equation
%*******************************************************************************
% Ternary system
%                (1): Ethanol
%                (2): Benzene
%                (3): Water
%*******************************************************************************
% Pressure is 1 atm = 1.013 bar = 760 mmHg
% Composition vectors xA, xB, are line vectors of NT columns
%*******************************************************************************
% Gas constant
R=1.98721;     % cal/(mol*K)
%*******************************************************************************
% Data that have to be changed for different mixtures

% Pure component constants for Antoine equation (Gmehling et al, 1977-1990)
A1=8.11220; B1=1592.864; C1=226.184;  % (T-range:   20C -  93C, Component 1)
A2=6.87987; B2=1196.760; C2=219.161;  % (T-range:   8C -  80C, Component 1)
A3=8.07131; B3=1730.630; C3=233.426;  % (T-range:   1C - 100C, Component 2)

% Interaction parameters Aij for NRTL 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=775.7625;               % Binary mixture 1-2 (cal/mol) 
A21=1118.9759;               % Binary mixture 1-2 (cal/mol) 
A13=-109.6339;               % Binary mixture 1-3 (cal/mol) 
A31=1332.3134;                % Binary mixture 1-3 (cal/mol) 
A23=1522.2402;               % Binary mixture 2-3 (cal/mol) 
A32=1821.8322;               % Binary mixture 2-3 (cal/mol) 

% Alpha interaction parameter
alpha11=0;
alpha22=0;
alpha33=0;
alpha12=0.6268;                    % Binary mixture 1-2 (cal/mol) 
alpha21=alpha12;                   % Binary mixture 2-1 (cal/mol) 
alpha13=0.3031;                    % Binary mixture 1-3 (cal/mol) 
alpha31=alpha13;                   % Binary mixture 3-1 (cal/mol) 
alpha23=0.3547;                    % Binary mixture 2-3 (cal/mol) 
alpha32=alpha23;                   % Binary mixture 3-2 (cal/mol) 


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

%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)

% Binary interaction parameters 'tafji' for NRTL equation 
taf11=A11./(R*T');            % Pure component 1 
taf22=A22./(R*T');            % Pure component 2 
taf33=A33./(R*T');            % Pure component 3 
taf12=A12./(R*T');            % Binary mixture 1-2 
taf21=A21./(R*T');            % Binary mixture 1-2 
taf13=A13./(R*T');            % Binary mixture 1-3  
taf31=A31./(R*T');            % Binary mixture 1-3  
taf23=A23./(R*T');            % Binary mixture 2-3  
taf32=A32./(R*T');            % Binary mixture 2-3 

% Binary interaction parameters 'Gji' for NRTL equation 
G11=exp(-alpha11*taf11);            % Pure component 1 
G22=exp(-alpha22*taf22);            % Pure component 2 
G33=exp(-alpha33*taf33);            % Pure component 3 
G12=exp(-alpha12*taf12);            % Binary mixture 1-2 
G21=exp(-alpha21*taf21);            % Binary mixture 1-2 
G13=exp(-alpha13*taf13);            % Binary mixture 1-3  
G31=exp(-alpha31*taf31);            % Binary mixture 1-3  
G23=exp(-alpha23*taf23);            % Binary mixture 2-3  
G32=exp(-alpha32*taf32);            % Binary mixture 2-3 

i=1:NT;

% Logarithm of activity coefficients 'LNGAMMA'
LNGAMMA1(i)=(taf11.*G11.*xA(i)+taf21.*G21.*xB(i)+taf31.*G31.*(1-xA(i)-xB(i)))./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i)))...
    + xA(i).*G11./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))).*(taf11 - (xA(i).*taf11.*G11+xB(i).*taf21.*G21+(1-xA(i)-xB(i)).*taf31.*G31)./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))))...
    + xB(i).*G12./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))).*(taf12 - (xA(i).*taf12.*G12+xB(i).*taf22.*G22+(1-xA(i)-xB(i)).*taf32.*G32)./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))))...
    + (1-xA(i)-xB(i)).*G13./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))).*(taf13 - (xA(i).*taf13.*G13+xB(i).*taf23.*G23+(1-xA(i)-xB(i)).*taf33.*G33)./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))));

LNGAMMA2(i)=(taf12.*G12.*xA(i)+taf22.*G22.*xB(i)+taf32.*G32.*(1-xA(i)-xB(i)))./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i)))...
    + xA(i).*G21./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))).*(taf21 - (xA(i).*taf11.*G11+xB(i).*taf21.*G21+(1-xA(i)-xB(i)).*taf31.*G31)./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))))...
    + xB(i).*G22./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))).*(taf22 - (xA(i).*taf12.*G12+xB(i).*taf22.*G22+(1-xA(i)-xB(i)).*taf32.*G32)./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))))...
    + (1-xA(i)-xB(i)).*G23./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))).*(taf23 - (xA(i).*taf13.*G13+xB(i).*taf23.*G23+(1-xA(i)-xB(i)).*taf33.*G33)./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))));
  
LNGAMMA3(i)=(taf13.*G13.*xA(i)+taf23.*G23.*xB(i)+taf33.*G33.*(1-xA(i)-xB(i)))./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i)))...
    + xA(i).*G31./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))).*(taf31 - (xA(i).*taf11.*G11+xB(i).*taf21.*G21+(1-xA(i)-xB(i)).*taf31.*G31)./(G11.*xA(i)+G21.*xB(i)+G31.*(1-xA(i)-xB(i))))...
    + xB(i).*G32./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))).*(taf32 - (xA(i).*taf12.*G12+xB(i).*taf22.*G22+(1-xA(i)-xB(i)).*taf32.*G32)./(G12.*xA(i)+G22.*xB(i)+G32.*(1-xA(i)-xB(i))))...
    + (1-xA(i)-xB(i)).*G33./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))).*(taf33 - (xA(i).*taf13.*G13+xB(i).*taf23.*G23+(1-xA(i)-xB(i)).*taf33.*G33)./(G13.*xA(i)+G23.*xB(i)+G33.*(1-xA(i)-xB(i))));

% Activity coefficients 'GAMMA'
GAMMA1(i)=exp(LNGAMMA1(i));                          % Component 1
GAMMA2(i)=exp(LNGAMMA2(i));                          % Component 2
GAMMA3(i)=exp(LNGAMMA3(i));                          % Component 3

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

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

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