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

%*******************************************************************************
% Program written by Stathis Skouras, 19/3/2000
%******************************************************************************* 
% Program for ternary VLE by using Wilson 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 mixtures

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

% Binary interaction parameters Aij for Wilson 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=135.8113;               % Binary mixture 1-2 (cal/mol) 
A21=-132.0576;                % Binary mixture 1-2 (cal/mol) 
A13=397.4862;                % Binary mixture 1-3 (cal/mol) 
A31=-265.4397;               % Binary mixture 1-3 (cal/mol) 
A23=102.6706;                 % Binary mixture 2-3 (cal/mol) 
A32=-40.8745;                % Binary mixture 2-3 (cal/mol) 

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

% Gas constant
R=1.98721;     % cal/(mol*K)
i=1:NT;

% Binary interaction parameters for Wilson equation 
L11=V1/V1*exp(-A11./(R*T'));            % Pure component 1 
L22=V2/V2*exp(-A22./(R*T'));            % Pure component 2 
L33=V3/V3*exp(-A33./(R*T'));            % Pure component 3 
L12=V2/V1*exp(-A12./(R*T'));            % Binary mixture 1-2 
L21=V1/V2*exp(-A21./(R*T'));            % Binary mixture 1-2 
L13=V3/V1*exp(-A13./(R*T'));            % Binary mixture 1-3  
L31=V1/V3*exp(-A31./(R*T'));            % Binary mixture 1-3  
L23=V3/V2*exp(-A23./(R*T'));            % Binary mixture 2-3  
L32=V2/V3*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)

% Activity coefficients by Wilson equation
G1(i) = exp(1-log(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xA(i).*L11(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L21(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L31(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i)));    % Component 1

G2(i) = exp(1-log(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...
          - xA(i).*L12(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L22(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L32(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i)));    % Component 2

G3(i) = exp(1-log(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i))...
          - xA(i).*L13(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L23(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L33(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(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];
