% *********************************************************************** %
% ******************* Calculates binodal curve ************************** %
% *********** Based on a code from Tore Haug-Warberg 2003 *************** %
% *********************************************************************** %

% Experimental data for Benzene(1) - Water (2) - 1,4-Dioxane (3)
% Sørensen and Arlt (1980), DECHEMA Data Series, T = 25 C, p = 1 bar.
% Aqueous phase
x1 = [ 0.0405   0.024    0.027   0.029   0.060    0.061    0.103   0.222   0.270   0.313   0.638   1.298    2.088];
x2 = [99.9590  98.888   95.421  93.521  92.043   91.938   87.505  84.475  82.717  82.076  78.207  71.933   67.230];
x3 = 100 - x1 - x2;
% Organic phase
y1 = [ 99.700  94.444  78.786  68.517  63.010  62.062   45.964   36.902  32.146  28.999  23.816  13.850  10.465];
y2 = [  0.290   0.435   0.884   2.645   3.095   3.098    6.571    9.894  11.900  14.204  18.225  30.752  39.086];
y3 = 100 - y1 - y2;

% Critical point (estimated from experimental data)
c=[0.03; 0.54; 0.43];

x  = [x1;x2;x3];         % experimental points for aqueous phase (initial values for the liquid-liquid split)
y  = [y1;y2;y3];         % experimental points for organic phase (initial values for the liquid-liquid split)
k  = ones(size(x));      % same ref.state in both phases

%The NRTL activity model is used for both liquid phases
liq     = @ActivityNRTL;                          % function name

% NRTL parameters from Sørensen and Arlt, 1980
s.alpha = [0    0.2  0.2;                         % non-random par
           0.2  0    0.2;
           0.2  0.2  0];                           

s.tau   = [   0.0    1411.4    129.49; ...        % interaction coefficient parameters [K]
           1320.6      0.0     525.41; ...            
            11.834   49.551       0.0];                 

s.t     = 25 + 273.15;                            % temperature [K]

% Calculate liquid-liquid equilibria
for i=1:max(size(x))
    [x(:,i),y(:,i)] = tpnewton(x(:,i),y(:,i),k(:,i),liq,s,liq,s);
end

% Triangle coordinates for the aqueous phase
tx1=x(1,:)';
tx3=1-x(1,:)'-x(2,:)';
tx2=2*sqrt(3)/3 - sqrt(3)/3*tx1 - 2*sqrt(3)/3*tx3;
plot(tx2,tx1,'bo')

hold on

% Triangle coordinates for the organic phase
ty1=y(1,:)';
ty3=1-y(1,:)'-y(2,:)';
ty2=2*sqrt(3)/3 - sqrt(3)/3*ty1 - 2*sqrt(3)/3*ty3;
plot(ty2,ty1,'bo')

% Show calculated tie-lines (for the experimental data)
plot([tx2';ty2'],[tx1';ty1'],'k-');

% Append the estimated critical end point to the composition arrays.
x(:,end+1) = c;
y(:,end+1) = c;

% Expand composition  arrays to  obtain  more plot points
a = linspace(1,1/99,5);                                     % interpolation vector
x = kron(x(:,1:end-1),a)+kron(x(:,2:end),1-a);              % liquid phase 1
y = kron(y(:,1:end-1),a)+kron(y(:,2:end),1-a);              % liquid phase 2
k  = ones(size(x));

%Calculate several more liquid-liquid equilibria points
for i=1:max(size(x))
   [x(:,i),y(:,i)] = tpnewton(x(:,i),y(:,i),k(:,i),liq,s,liq,s);
end

%--------------------------------------------------------------------------
% 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.095,1.07,'Benzene [s]'); h2=get(h,'Children');set(h2(1),'FontWeight','bold'); 
text(sqrt(3)/3-0.055,1.035,'80.1 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);
text(2*sqrt(3)/3-0.08,-0.04,'Water [sn]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-0.06,-0.08,'100.0 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);
text(-0.14,-0.04,'1,4-Dioxane [sn]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(-0.060,-0.08,'101.3 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);

text(sqrt(3)/3+0.03,0.99,'x_R');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3+0.005,0.008,'x_W');h2=get(h,'Children');set(h2(1),'FontWeight','bold');

text(2*sqrt(3)/3-1.35,0.90,'- - - -   binodal curve at 25 ^oC');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-1.35,0.80,'o o o  experimental points ');h2=get(h,'Children');set(h2(1),'FontWeight','bold','color','b');
text(2*sqrt(3)/3-1.35,0.70,'_____  calculated tie-lines');h2=get(h,'Children');set(h2(1),'FontWeight','bold','color','g');
text(2*sqrt(3)/3-1.35,0.60,'o         critical point');h2=get(h,'Children');set(h2(1),'FontWeight','bold','color','m');

plot(X,Y,'k') 
%--------------------------------------------------------------------------

% Triangle coordinates for the aqueous phase
tx1=x(1,:)';
tx3=1-x(1,:)'-x(2,:)';
tx2=2*sqrt(3)/3 - sqrt(3)/3*tx1 - 2*sqrt(3)/3*tx3;
%plot(tx2,tx1,'ko')

% Triangle coordinates for the organic phase
ty1=y(1,:)';
ty3=1-y(1,:)'-y(2,:)';
ty2=2*sqrt(3)/3 - sqrt(3)/3*ty1 - 2*sqrt(3)/3*ty3;
%plot(ty2,ty1,'ko')

% Show estimated critical point
tc1=c(1);
tc3=1-c(1)-c(2);
tc2=2*sqrt(3)/3 - sqrt(3)/3*tc1 - 2*sqrt(3)/3*tc3;
plot(tc2,tc1,'mo')

% Show calculated tie-lines
plot([tx2';ty2'],[tx1';ty1'],'g');

% Show liquid-liquid binodal curve 
plot(tx2,tx1,'k--')
plot(ty2,ty1,'k--')
