% ------------------------------%
% TERNARY COMPOSITION DIAGRAMS  %
% Stathis Skouras, January 2002 %
% ----------------------------- %
% Number of stages in the rectifying section: 
Nr = 51;
% Total number of 'stages'  
NT = Nr + 2;                                  % We include condenser and reboiler 

% 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.61,0.28,'het.az [un]');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-0.57,0.24,'66.1 ^oC'); h3=get(h,'Children');set(h3(1),'FontSize',10);

text(2*sqrt(3)/3-0.55,0.60,'F');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
%text(0.37,0.77,'x_B');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
%text(2*sqrt(3)/3-0.39,0.71,'x_T');h2=get(h,'Children');set(h2(1),'FontWeight','bold');

%text(2*sqrt(3)/3-0.79,0.61,'still path');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
%text(2*sqrt(3)/3-0.30,0.36,'top');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
%text(2*sqrt(3)/3-0.32,0.31,'vessel');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
%text(2*sqrt(3)/3-1.30,0.85,'Serafimov^,s class 3.1-2');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
%text(2*sqrt(3)/3-1.4,0.85,'------  binodal curve at 25 ^oC');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-1.4,0.75,'.......  still path');h2=get(h,'Children');set(h2(1),'FontWeight','bold');
text(2*sqrt(3)/3-1.4,0.65,'-o-o-  column liquid profile');h2=get(h,'Children');set(h2(1),'FontWeight','bold');

plot(X,Y,'k') 

% --------------------------------------------
% Feed composition
% --------------------------------------------
Xf1=0.6;
Xf2=0.2;
tf1=Xf1;
tf3=1-Xf1-Xf2;
tf2=2*sqrt(3)/3 - sqrt(3)/3*tf1 - 2*sqrt(3)/3*tf3;
plot(tf2,tf1,'ko')

load x
% --------------------------------------------
% Top vessel path
% ---------------------------------------------
% Light component
Xt1=x(:,1);
% Intermediate component
Xt2=x(:,NT+1);
% Creates the 'triangle' compositions
% Light component
tx1=Xt1;
% Heavy component
tx3=1-Xt1-Xt2;
% Intermediate component
tx2=2*sqrt(3)/3 - sqrt(3)/3*tx1 - 2*sqrt(3)/3*tx3;
%plot(tx2,tx1,'r.')


% --------------------------------------------
% Still path
% ---------------------------------------------
% Light component
Xb1=x(:,NT);
% Intermediate component
Xb2=x(:,2*NT);
% Creates the 'triangle' compositions
% Light component
tx1=Xb1;
% Heavy component
tx3=1-Xb1-Xb2;
% Intermediate component
tx2=2*sqrt(3)/3 - sqrt(3)/3*tx1 - 2*sqrt(3)/3*tx3;
plot(tx2,tx1,'k.')


load Xinit

% --------------------------------------------
% Column composition profile
% ---------------------------------------------

% Light component
Xp1=Xinit(1:NT);
% Intermediate component
Xp2=Xinit(NT+1:2*NT);

% Creates the 'triangle' compositions
% Light component
tx1=Xp1;
% Heavy component
tx3=1-Xp1-Xp2;
% Intermediate component
tx2=2*sqrt(3)/3 - sqrt(3)/3*tx1 - 2*sqrt(3)/3*tx3;

plot(tx2,tx1,'k-.o')