function DP=bp_wilson(T,xA,xB,NT,TB1,TB2,TB3,V1,V2,V3,A1,A2,A3,B1,B2,B3,C1,C2,C3,A11,A22,A33,A12,A21,A13,A31,A23,A32,R)
% *********************************************************************************************************************
% ******************************** Program written by Stathis Skouras, March 2001 *************************************
% *********************************************************************************************************************

% This subroutine is for bubble point calculation in each stage  

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;
