% This program calculates distillation lines for specific composition points for ternary mixtures.
% Isobaric bubble point calculations  
% VAPOR-LIQUID EQUILIBRIUM for nonideal multicomponent mixtures 
%                                                      
%      y_i P = x_i Psat_i(T) gamma_i(T,x)              
%                                                      
% 1. Specify P and x              
% 2. Guess T                                                   
% 3. Calculate Psat_i, i=1,...,nc from Antoine equation      
% 4. Calculate gamma_i, i=1,...,nc from Wilson equation 
%    (or another liquid activity coefficient model) 
% 5. Calculate y_i from VLE equation above  
% 6. If sum(y_i) is sufficiently close to unity, then stop.      
%    Otherwise, adjust T and go to step 3.
global X P NC R W A alpha 

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

NC=3;   % number of components
comp1='Ethanol';
comp2='Benzene';
comp3='Water';

P=1;    % pressure in atm

% Antoine coefficients where Psat_i [mmHg] and T [K]
A=[8.11220  1592.864  226.184     
   6.87987  1196.760  219.161
   8.07131  1730.630  233.426];

% NRTL parameters (small parameters, Wij)
% where V [cm3/mol], W_ij [cal/mol] and T [K]
W=[   0        775.7625    -109.6339
  1118.9759     0          1522.2402
  1332.3134   1821.8322       0     ];

% Alpha parameter
alpha=[0        0.6268   0.3031
       0.6268   0        0.3547      
       0.3031   0.3547   0     ];

%---------------------%
% DISTILLATION LINES  %
%---------------------%
X = [0 2*sqrt(3)/3];
Y = [0 0];
plot(X,Y,'k') 
hold on
axis off
axis([0 2*sqrt(3)/3 0 1])
axis equal
X = [0 sqrt(3)/3];
Y = [0 1];
plot(X,Y,'k') 
X = [sqrt(3)/3 2*sqrt(3)/3]; 
Y = [1 0];  
plot(X,Y,'k') 

% Creates the ternary diagram
X = [0 2*sqrt(3)/3];
Y = [0 0];
plot(X,Y,'k') 
hold on
axis off
axis([0 2*sqrt(3)/3 0 1])
axis equal
X = [0 sqrt(3)/3];
Y = [0 1];
plot(X,Y,'k') 
X = [sqrt(3)/3 2*sqrt(3)/3]; 
Y = [1 0];  
h=gca;
text(sqrt(3)/3-0.1,1.07,'Ethanol [sn]'); h2=get(h,'Children');set(h2(1),'FontWeight','bold'); 
text(sqrt(3)/3-0.055,1.035,'78.3 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);
text(2*sqrt(3)/3-0.12,-0.04,'Benzene [sn]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-0.06,-0.08,'80.1 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);
text(-0.1,-0.04,'Water [sn]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(-0.060,-0.08,'100 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);

text(2*sqrt(3)/3-0.25,0.5,'hom.az [s]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-0.20,0.46,'68 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);
text(2*sqrt(3)/3-0.82,0.97,'hom.az [s]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-0.79,0.93,'78.2 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);

text(2*sqrt(3)/3-0.45,-0.04,'het.az [s]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-0.41,-0.08,'69 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);
text(2*sqrt(3)/3-0.60,0.32,'het.az [un]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-0.57,0.28,'66.1 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);

%text(2*sqrt(3)/3-1.30,0.75,'Serafimov^,s class 3.1-2');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
%text(2*sqrt(3)/3-1.40,0.55,'------  distillation boundaries');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
%text(2*sqrt(3)/3-1.40,0.45,'-.-.-.  batch distillation');h2=get(h,'Children');set(h2(1),'FontWeight','bold');%'color','r');
%text(2*sqrt(3)/3-1.30,0.40,'boundary');h2=get(h,'Children');set(h2(1),'FontWeight','bold');%'color','r');

%text(2*sqrt(3)/3-1.40,0.55,'......  distillation lines');h2=get(h,'Children');set(h2(1),'FontWeight','bold','color','r');
%text(2*sqrt(3)/3-0.60,0.69,'*');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
%text(2*sqrt(3)/3-0.62,0.73,'F');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
plot(X,Y,'k') 

% --------------------------------------------
% Sketch the batch distillation boundary
% --------------------------------------------
XL1=[0.2514;0];
XL2=[0.5356;0];
tf1=XL1;
tf3=1-XL1-XL2;
tf2=2*sqrt(3)/3 - sqrt(3)/3*tf1 - 2*sqrt(3)/3*tf3;
%plot(tf2,tf1,'k-.')

% --------------------------------------------
% Ternary azeotrope composition
% --------------------------------------------
Xf1=0.2514;
Xf2=0.5356;
tf1=Xf1;
tf3=1-Xf1-Xf2;
tf2=2*sqrt(3)/3 - sqrt(3)/3*tf1 - 2*sqrt(3)/3*tf3;
plot(tf2,tf1,'ko')

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%% Sketch the 3 distillation boundaries %%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initial temperatures to start the bubble point calculations (intermediate component)
T_K=353.25;

% Prepare xy equilibrium vector 
txy1=ones(2,101)';
txy2=ones(2,101)';

Y=[0.000000001 0.01 0.999999999]'; % initial value (k=34)
for i=1:101           % number of vectors/steps
  X=Y;
  T_K=fzero('sumy',T_K);
  temp_T=T_K-273.15;
  Y=ycomp(T_K);
  xy1=[X(1) Y(1)];
  xy2=[X(2) Y(2)];
  xy3=[X(3) Y(3)];
  txy1(i,:)=xy1;
  txy2(i,:)=2*sqrt(3)/3 - sqrt(3)/3*xy1 - 2*sqrt(3)/3*xy3;
end
plot(txy2,txy1,'k-.')

Y=[0.9995 0.0004 0.0001]';  % initial value  (k=18)
for i=1:101           % number of vectors/steps
  X=Y;
  T_K=fzero('sumy',T_K);
  temp_T=T_K-273.15;
  Y=ycomp(T_K);
  xy1=[X(1) Y(1)];
  xy2=[X(2) Y(2)];
  xy3=[X(3) Y(3)];
  txy1(i,:)=xy1;
  txy2(i,:)=2*sqrt(3)/3 - sqrt(3)/3*xy1 - 2*sqrt(3)/3*xy3;
end
plot(txy2,txy1,'k-.')


Y=[0.92 0.00001 0.07999]';  % initial value  (k=18)
for i=1:101           % number of vectors/steps
  X=Y;
  T_K=fzero('sumy',T_K);
  temp_T=T_K-273.15;
  Y=ycomp(T_K);
  xy1=[X(1) Y(1)];
  xy2=[X(2) Y(2)];
  xy3=[X(3) Y(3)];
  txy1(i,:)=xy1;
  txy2(i,:)=2*sqrt(3)/3 - sqrt(3)/3*xy1 - 2*sqrt(3)/3*xy3;
end
plot(txy2,txy1,'k-.')

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sketching the distillation curves
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y=[0.4 0.001 0.599]';  % initial value  (k=18)
for i=1:101           % number of vectors/steps
  X=Y;
  T_K=fzero('sumy',T_K);
  temp_T=T_K-273.15;
  Y=ycomp(T_K);
  xy1=[X(1) Y(1)];
  xy2=[X(2) Y(2)];
  xy3=[X(3) Y(3)];
  txy1(i,:)=xy1;
  txy2(i,:)=2*sqrt(3)/3 - sqrt(3)/3*xy1 - 2*sqrt(3)/3*xy3;
end
%plot(txy2,txy1,'r-.')


Y=[0.015 0.001 0.984]';  % initial value  (k=18)
for i=1:101           % number of vectors/steps
  X=Y;
  T_K=fzero('sumy',T_K);
  temp_T=T_K-273.15;
  Y=ycomp(T_K);
  xy1=[X(1) Y(1)];
  xy2=[X(2) Y(2)];
  xy3=[X(3) Y(3)];
  txy1(i,:)=xy1;
  txy2(i,:)=2*sqrt(3)/3 - sqrt(3)/3*xy1 - 2*sqrt(3)/3*xy3;
end
%plot(txy2,txy1,'r-.')


Y=[0.005 0.97 0.025]';  % initial value  (k=18)
for i=1:101           % number of vectors/steps
  X=Y;
  T_K=fzero('sumy',T_K);
  temp_T=T_K-273.15;
  Y=ycomp(T_K);
  xy1=[X(1) Y(1)];
  xy2=[X(2) Y(2)];
  xy3=[X(3) Y(3)];
  txy1(i,:)=xy1;
  txy2(i,:)=2*sqrt(3)/3 - sqrt(3)/3*xy1 - 2*sqrt(3)/3*xy3;
end
%plot(txy2,txy1,'r-')