% *********************************************************************** %
% ******************* Calculates binodal curve ************************** %
% *********************************************************************** %
% Experimental data for Ethanol(1) - Benzene(2) - Water (3)
% Aqueous (left) phase
x1 = [0.0001  3.817  7.968  12.977  18.134  23.540  24.069  27.892  31.725  35.510  39.382  41.062  41.771];
x2 = [0.0405  0.049  0.052   0.139   0.388   1.001   1.107   1.828   3.371   5.395   9.584  13.309  21.053];
x3 = 100 - x1 - x2;
% Organic (right) phase
y1 = [0.0001   1.010   3.323   5.860   9.121  12.939  13.340  16.090  18.943  22.444  26.216  29.341  33.093];
y2 = [99.700  98.560  95.827  92.059  88.424  83.474  82.694  78.864  74.623  69.829  63.551  58.096  50.300];
y3 = 100 - y1 - y2;
% Critical point (estimated from experimental data)
c=[0.393; 0.367; 0.240];

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

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

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

s.alpha = [0      0.2    0.2;                   % non-random par
           0.2    0      0.2;                   % from LLE data
           0.2    0.2    0];              

s.tau   = [  0.0      87.744   -441.74; ...         % interaction coefficient parameters [K]
           118.04      0.0       986.99; ...         % from LLE data
           376.33   2797.7         0.0];            

s.t     = 25+273.15;                                 % temperature [K]
 
%Calculate liquid-liquid equilibria (for the experimental 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.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.65,'.......  binodal curve at 25 ^oC');h2=get(h,'Children');set(h2(1),'FontWeight','bold');

plot(X,Y,'k')

% --------------------------------------------
% Triangle approximation
% ---------------------------------------------

X1left=[0,0.64];
X2left=[0,0.07];
t1left=X1left;
t3left=1-X1left-X2left;
t2left=2*sqrt(3)/3 - sqrt(3)/3*t1left - 2*sqrt(3)/3*t3left;
%plot(t2left,t1left,'c-')

X1right=[0,0.64];
X2right=[1,0.07];
t1right=X1right;
t3right=1-X1right-X2right;
t2right=2*sqrt(3)/3 - sqrt(3)/3*t1right - 2*sqrt(3)/3*t3right;
%plot(t2right,t1right,'c-')


% Triangle coordinates for the aqueous phase (experimental points)
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 (experimental points)
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 calculated tie-lines (experimental points)
% plot([tx2';ty2'],[tx1';ty1'],'k-');

%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,'co')

%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 2-phase 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,'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 calculated tie-lines
plot([tx2';ty2'],[tx1';ty1'],'k-.');

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