% This program calculates and plots VLE phase diagrams 
% (T-x,y and y,x etc.) for binary- and ternary mixtures. 
%
%     
% 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 V A


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

NC=3;   % number of components
comp1='methanol';
comp2='ethanol';
comp3='1-propanol';

P=1;    % pressure in atm

% Wilson coefficients at 1 atm collected from the DECHEMA series 
% vol.1/2a, p.60,120,336 and Antoine coefficients from p.116.

% Antoine coefficients where Psat_i [mmHg] and T [K]
A=[8.08097  1582.271  239.726     
   8.11220  1592.864  226.184  
   8.37895  1788.020  227.438];

% Wilson parameters (small parameters, Wij)
% where V [cm3/mol], W_ij [cal/mol] and T [K]
W=[  0         135.8113   397.4862
   -132.0576     0        102.6706 
   -265.4397   -40.8745      0     ];

% Molar volumes of pure liquid [cm3/mol] from 
% DECHEMA series vol.1/2a, p.705-706
V=[40.73  58.68 75.14];

% Set up column vectors for composition and temperature
% Ternary mixture of comp1, comp2 and comp3
% Pure components boiling points, Tbp=[64.5 78.4 97.2]
x1=ones(101,101)*diag(0:100)/100;
temp=ones(101,101)*351.55-273.15;   % initialization

for i=0:100
  for j=0:100
    x2(i+1,j+1)=(100-max(100-i,j))/100;
    x3(i+1,j+1)=1-x1(i+1,j+1)-x2(i+1,j+1);
    X=[x1(i+1,j+1) x2(i+1,j+1) 1-x1(i+1,j+1)-x2(i+1,j+1)]';
      if i==0 & j==0
	TK(i+1,j+1)=fzero('sumy',370.35); % Boiling point of 1-propanol
      else
	TK(i+1,j+1)=fzero('sumy',temp(i+1,j+1)+273.15);
      end
    temp(i+1,j+1)=TK(i+1,j+1)-273.15;
    Y=ycomp(TK(i+1,j+1));
    y1(i+1,j+1)=Y(1);
    y2(i+1,j+1)=Y(2);
    y3(i+1,j+1)=Y(3);
  end
end


x1=ones(401,401)*diag(0:400)/400;
temp=ones(401,401)*351.55-273.15;   % initialization

for i=0:400
  for j=0:400
    x2(i+1,j+1)=(400-max(400-i,j))/400;
    x3(i+1,j+1)=1-x1(i+1,j+1)-x2(i+1,j+1);
    X=[x1(i+1,j+1) x2(i+1,j+1) 1-x1(i+1,j+1)-x2(i+1,j+1)]';
      if i==0 & j==0
	TK(i+1,j+1)=fzero('sumy',370.35); % Boiling point of 1-propanol
      else
	TK(i+1,j+1)=fzero('sumy',temp(i+1,j+1)+273.15);
      end
    temp(i+1,j+1)=TK(i+1,j+1)-273.15;
    Y=ycomp(TK(i+1,j+1));
    y1(i+1,j+1)=Y(1);
    y2(i+1,j+1)=Y(2);
    y3(i+1,j+1)=Y(3);
  end
end


for i=0:100
  for j=0:100
     if x1(i+1,j+1)==0
       K_1(i+1,j+1)=1;
     else
       K_1(i+1,j+1)=y1(i+1,j+1)/x1(i+1,j+1);
     end
     if x2(i+1,j+1)==0
       K_2(i+1,j+1)=1;
     else
       K_2(i+1,j+1)=y2(i+1,j+1)/x2(i+1,j+1);
     end
     if x3(i+1,j+1)==0
       K_3(i+1,j+1)=1;
     else
       K_3(i+1,j+1)=y3(i+1,j+1)/x3(i+1,j+1);
     end
       alpha_12(i+1,j+1)=K_1(i+1,j+1)/K_2(i+1,j+1);
       alpha_13(i+1,j+1)=K_1(i+1,j+1)/K_3(i+1,j+1);
       alpha_23(i+1,j+1)=K_2(i+1,j+1)/K_3(i+1,j+1);
  end
end

% Transformation into LIKESIDET TREKANT (60o-60o-60o)
% hvor hoeyden er skalert til 1 (dvs. 100% )
t1=x1;
t2=2*sqrt(3)/3 - sqrt(3)/3*x1 - 2*sqrt(3)/3*x3;


% Plot singular points
for i=0:100
  for j=0:100
     if alpha_12(i+1,j+1)==1
       plot(t2(i+1,j+1),t1(i+1,j+1),'mo')
     end
  end
end
  
for i=0:100
  for j=0:100
     if x1(i+1,j+1)==y1(i+1,j+1) & x2(i+1,j+1)==y2(i+1,j+1)...
	   & x3(i+1,j+1)==y3(i+1,j+1) 
       plot(t2(i+1,j+1),t1(i+1,j+1),'m*')
     end
  end
end
         
%---------------------%
% EQUILIBRIUM VECTORS %
%---------------------%
t1=x1;
t2=2*sqrt(3)/3 - sqrt(3)/3*x1 - 2*sqrt(3)/3*x3;

figure(1)
X = [0 2*sqrt(3)/3];
Y = [0 0];
plot(X,Y,'k') 
hold
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') 

txy1=ones(2,101)';
txy2=ones(2,101)';
T_K=351.55;

Y=[0.01 0.01 0.98]';  % initial value  (k=18)
%Y=[0.01 0.03 0.96]';  % initial value (k=20)
%Y=[0.01 0.06 0.93]';  % initial value (k=21)
%Y=[0.01 0.1 0.89]';   % initial value (k=22)
%Y=[0.01 0.2 0.79]';   % initial value (k=23)
%%Y=[0.01 0.3 0.69]';   % initial value (k=24)
%Y=[0.001 0.2 0.799]'; % initial value (k=28)
%%Y=[0.001 0.1 0.899]'; % initial value (k=26)
%%%Y=[0.0001 0.05 0.9499]'; % initial value (k=26)
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,'ko')

%contour(t2,t1,temp,30,'b:')
%contour(t2,t1,temp,30)

print -depsc distlines_dot_isotherms_c_mep.eps
print -depsc distlines_isotherms_c_mep.eps

print -deps dist_vector_lines_mep.eps   % balck/white
print -depsc c_dist_vector_lines_mep.eps  % colour


%---------------------------------------------------%
% SPLINES cubic interpolation to make smooth curves %
%---------------------------------------------------%
%clear k
for i=1:max(size(txy1(:,1)))
  if txy1(i,1) == 1
    k=max(size(txy1(1:i,1)));
  end
end
s_txy1=(0:.25:100)/100;
s_txy2=spline(txy1(1:k,1)',txy2(1:k,1)',s_txy1);
plot(s_txy2,s_txy1,'k')

print -depsc dot_distlines_isotherms_c_mep.eps  
%print -deps distcurves_mep.eps

mmtext('Methanol 64.5^oC ')
mmtext('Ethanol 78.4^oC ')
mmtext('1-Propanol 97.2^oC ')

%s_tx1=ones(3,101);
%s_tx2=ones(3,101);
%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=xy1;
%  txy2=2*sqrt(3)/3 - sqrt(3)/3*xy1 - 2*sqrt(3)/3*xy3;
%  plot(txy2,txy1,'r-')
%  s_tx1(:,i)=X;
%  s_tx2(:,i)=Y;
%end


%----------------------------------------
% Alpha=unity (univolatility lines) plot
%----------------------------------------
% For a_12, plot i=2:126
% For a_13, plot i=2:

figure(2)
X = [0 2*sqrt(3)/3];
Y = [0 0];
plot(X,Y,'k') 
hold
axis off
axis([0 2*sqrt(3)/3 0 1])
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') 

% univolatility line plotting needs to be "cleaned"
% (not to plot alpha-lines along the binary axis)

[a_12,h]=contour(t2,t1,alpha_12,[1 1],'c-.');

%contour(t2,t1,alpha_13,[1 1],'c-.');
%contour(t2,t1,alpha_23,[1 1],'c-.');

% K-unity lines
contour(t2,t1,K_1,[1 1],'b-');      
contour(t2,t1,K_2,[1 1],'b-');
contour(t2,t1,K_3,[1 1],'b-');

%S=max(size(a_12));
%unity=ones(2,S);

%clear c1
%clear c2
%clear c3
%clear g2_c1
%clear g1_c2

%for i=2:S 
% c1(i-1)=c(1,i); 
% c2(i-1)=c(2,i); 
% c3(i-1)=unity(1,i)-c(1,i)-c(2,i);
%end
  
%g2_c1=c1;
%g1_c2=2*sqrt(3)/3 - sqrt(3)/3*c1 - 2*sqrt(3)/3*c3;
%plot(g1_c2,g2_c1,'m-.')  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3D TEMPERATURE - TRIANGULAR COMPOSITION DIAGRAM
% PLOT LIKESIDET TREKANT (60o-60o-60o)
% hvor hoeyden er skalert til 1 (dvs. 100% )
%
% Temperature on z-axis
% Pure methanol   (1) upper corner
% Pure ethanol    (2) lower right corner
% Pure 1-propanol (3) lower left corner
% Pure components boiling points, Tbp=[64.5 78.4 97.2]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t1=x1;
t2=2*sqrt(3)/3 - sqrt(3)/3*x1 - 2*sqrt(3)/3*x3;

figure(1)
X = [0 2*sqrt(3)/3];
Y = [0 0];
plot(X,Y,'k') 
hold
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') 

%contour(t2,t1,temp,20)

[tt,h]=contour(t2,t1,temp,[71.45 71.45]);
[tt,h]=contour(t2,t1,temp,[88.10 88.10]);

[tt,h]=contour(t2,t1,temp,[67.0 67.0]);
clabel(tt,h,[67.0 67.0]) 

[tt,h]=contour(t2,t1,temp,[84.5 84.5]);  
clabel(tt,h,[84.5 84.5])

print -depsc isotherms_c_mep.eps


g2_x1=x1;
g1_x2=2*sqrt(3)/3 - sqrt(3)/3*x1 - 2*sqrt(3)/3*x3;;

mesh(g1_x2,g2_x1,temp)  

hold
axis off
axis equal
X = [0 2*sqrt(3)/3 0];
Y = [0 0 0];
line(X,Y) 
X = [0 sqrt(3)/3 0];
Y = [0 1 0];
line(X,Y) 
X = [sqrt(3)/3 2*sqrt(3)/3 0]; 
Y = [1 0 0];  
line(X,Y) 

grid_triangular

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TRIANGULAR COMPOSITION DIAGRAM
% PLOT LIKESIDET TREKANT (60o-60o-60o)
% hvor hoeyden er skalert til 1 (dvs. 100% )
% Temperature on z-axis
% Pure methanol   (1) upper corner
% Pure ethanol    (2) lower right corner
% Pure 1-propanol (3) lower left corner
% Pure components boiling points, Tbp=[64.5 78.4 97.2]

g2_x1=x1;
g1_x2=2*sqrt(3)/3 - sqrt(3)/3*x1 - 2*sqrt(3)/3*x3;;

contour(g1_x2,g2_x1,temp) 
hold
axis off
axis equal
X = [0 2*sqrt(3)/3];
Y = [0 0];
line(X,Y) 
X = [0 sqrt(3)/3];
Y = [0 1];
line(X,Y) 
X = [sqrt(3)/3 2*sqrt(3)/3]; 
Y = [1 0];  
line(X,Y) 

grid_triangular

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOTTING
% Ternary T-x1,x2 diagrams
% Use MESH to make 3-D plots of surfaces
% Use CONTOUR to make 2-D contour plots
$
% contour(x1,x2,temp,[58 60 62 64 66 68 70]) % plot lines of constant temp
% contour(x1,x2,temp,20) % 20 lines in contour plot
% [c,h] = contour(x1,x2,temp)

figure(3)
mesh(x2,x1,temp)
%view(0,90)
%view(0,0)
%view(-37.5,160)
view(-45,160)
hold
%contour(x2,x1,temp,50) % plots 50 temp contour lines 
%contour(x2,x1,temp,[80 80]) %plots a single contour at the level 80

X = [1 0];
Y = [0 1];
line(X,Y) 
X = [1 0]; 
Y = [0 0];  
line(X,Y) 
X = [0 0]; 
Y = [0 1]; 
line(X,Y) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 90oC composition triangle
% methanol    (1) upper left  corner
% ethanol     (2) lower right corner
% 1-propanol  (3) lower left  corner
%
% Pure components boiling points, Tbp=[64.5 78.4 97.2]
%
contour(x2,x1,temp,[64 64]);
contour(x2,x1,temp,[65 65]);
contour(x2,x1,temp,[67 67]);
contour(x2,x1,temp,[69 69]);
contour(x2,x1,temp,[73 73]);
contour(x2,x1,temp,[80 80]);
contour(x2,x1,temp,[90 90]);

print -deps isothrems_amw.eps


figure(2)
contour(x2,x1,temp,50)
hold
axis([0 1 0 1 0 100]) 
axis off
X = [1 0];
Y = [0 1];
line(X,Y) 
X = [1 0]; 
Y = [0 0];  
line(X,Y) 
X = [0 0]; 
Y = [0 1]; 
line(X,Y) 

plot(0,0,'+') 
plot(0.1,0,'+')
plot(0.2,0,'+')
plot(0.3,0,'+')
plot(0.4,0,'+') 
plot(0.5,0,'+')
plot(0.6,0,'+') 
plot(0.7,0,'+')
plot(0.8,0,'+') 
plot(0.9,0,'+')
plot(1,0,'+') 
plot(0,0.1,'+') 
plot(0,0.2,'+') 
plot(0,0.3,'+') 
plot(0,0.4,'+') 
plot(0,0.5,'+') 
plot(0,0.6,'+') 
plot(0,0.7,'+') 
plot(0,0.8,'+') 
plot(0,0.9,'+') 
plot(0,1,'+') 

header1=['Isotherms for ',comp1,'/',comp2,'/',comp3,' at ',num2str(P),'atm'];
title(header1)

% mmtext('Methanol')
% mmtext('Ethanol')    
% mmtext('1-Propanol')    
% mmtext('Isotherms for ternary mixture of methanol,... 
% ethanol and water at 1 atm')
% mmtext('Az')
% mmzap



% Binary T-x,y diagrams
% figure(1)
% plot(xone,temp,'r-',yone,temp,'r-')
% header1=['T-x,y diagram for ',comp1,'/',comp2,' at ',num2str(P),'atm'];
% title(header1)
% label1=['Molefraction of ',comp1];
% xlabel(label1)
% ylabel('Temperature [oC]')
% axis('square')
% grid
%
% figure(2)
% plot(xone,yone,'b-',xone,xone,'k-')
% header2=['y-x diagram for ',comp1,'/',comp2,' at ',num2str(P),'atm']
% title(header2)
% label2=['Liquid molefraction of ',comp1];
% xlabel(label2)
% label3=['Vapor molefraction of ',comp1];
% ylabel(label3)
% axis('square')
% grid







