%********************************************************************************************************************
%********************** Program written by Stathis Skouras, January 2005 **********************************************
%*********************** Based on previous programs from Katrine Hilmen ************************************************
%********************************************************************************************************************

% This program does an isobaric bubble point calculation to calculate equilibrium temperature and vapour compositions
% Given: X, P (liquid compositions, pressure) -> Calculated: Y, T (vapour compositions, temperature)
% The program uses the UNIQUAC activity coefficient model
% The program can handle any number of components NC
% The example here is for the quaternary mixture Methanol/Ethanol/1-Propanol/1-Butanol 

global X P NC R W V A RP Q

R=1.98721;   % universal gas constant [cal/molK]

NC=4;   % number of components

% Mixture information (lighter -> heavier)
comp1='Methanol';    
comp2='Ethanol';     
comp3='1-Propanol';  
comp4='1-Butanol';

P=1;    % pressure in atm

% *******************************************************************************
% Data that have to be changed for different mixtures (Dortmund Data Bank, Gmehling et al, 1977-1990)

% Antoine coefficients
A=[8.08097  1582.271  239.726
   8.11220  1592.864  226.184 
   8.37895  1788.020  227.438
   7.36366  1305.198  173.427];

% UNIQUAC parameters
W=[   0        -181.2860      107.1719      1088.8272
   247.3775        0          -35.2815      -157.4798
    -2.0418      71.7539        0            449.0943
  -432.2175     656.5578     -338.6701         0     ];

% Molar volumes of pure components
V=[40.73  58.68  75.14  91.97]; 

% Volume parameters of pure components 
RP=[1.4311  2.1055  2.7799  3.4543];

% Area parameters of pure components 
Q=[1.4320  1.9720  2.5120  3.0520];

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

% Liquid compositions (lighter -> heavier)
X1=0.25;       
X2=0.25;
X3=0.25;

% Vector of liquid compositions
X=[X1;X2;X3;1-X1-X2-X3]; 

% Initial temperature estimation (let's say 75C) (T is the iterating variable)
T = 75 +  273.15;   % [K]

% Function ''fzero('sumY',value) returns the T value (equilibrium temperature) that gives 1-sumY=0
T=fzero('sumY_UNIQUAC',T)        % Equilibrium temperature [K]

% Read vapour compositions at equilibrium
Y=UNIQUAC(T)                % Equilibrium vector of vapour compositions

% *************************************************************************
% % Here we can compare the model with experimental data from Gmehling et al. (1977-1990)
% *************************************************************************

% Experimental data
% Texp = ...;   %C;
% Y1exp = ...;
% Y2exp = ...;
% Y3exp=...; 

% Comparison
%DIFFT = Texp - T + 273.15
%DIFFY1=Y1exp-Y(1)
%DIFFY2=Y2exp-Y(2)
%DIFFY3=Y3exp-Y(3)