function plot_triangular(x1,x2,op1,op2)

%function plot_triangular(x1,x2)
%function plot_triangular(x1,x2,line_type) 
%function plot_triangular(x1,x2,line_type,grid_off)
%function plot_triangular(x1,x2,grid_off)
%
%Plot composition trajectory x1,x2 in triangular plot
%x1:A composition
%x2:B composition
%
%Use hold to plot several curves in same plot
%Use grid_off=1 to turn off grid
%Use line_type as in ordinary plot commands ('-.')
%
%Example: plot_triangular(x1,x2,'-r',1)  

%Written by I. Halvosen
%Based on scripts by K. Hilmen


if(nargin>=3)
  if(isstr(op1) & nargin==3)
    grid_off=0;
    line_type=op1;
  else
    grid_ff=op1;
  end  
  if (nargin>3)
    line_type=op1;
    grid_off=op2;
  end
else
  grid_off=0;%default
  line_type=[];
end

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


if(isempty(line_type))
  plot(t2,t1);
else
  plot(t2,t1,line_type);
end

if (ishold)
  return;
end


%Plot triangle
hold on

X = [0 2*sqrt(3)/3];
Y = [0 0];
plot(X,Y,'k') 

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') 

%text(0.57,1.02,'A')
%text(1.18,0.01,'B')
%text(-0.05,0.01,'C')

if(~grid_off)
  grid_triangular;
end
hold off;




function grid_triangular()

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0.1 - 2*sqrt(3)/3*0.9;
gx2_1 = 0.1;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.1;
gx2_2 = 0.1;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0.2 - 2*sqrt(3)/3*0.8;
gx2_1 = 0.2;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.2;
gx2_2 = 0.2;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0.3 - 2*sqrt(3)/3*0.7;
gx2_1 = 0.3;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.3;
gx2_2 = 0.3;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0.4 - 2*sqrt(3)/3*0.6;
gx2_1 = 0.4;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.4;
gx2_2 = 0.4;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0.5 - 2*sqrt(3)/3*0.5;
gx2_1 = 0.5;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.5;
gx2_2 = 0.5;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0.6 - 2*sqrt(3)/3*0.4;
gx2_1 = 0.6;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.6;
gx2_2 = 0.6;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0.7 - 2*sqrt(3)/3*0.3;
gx2_1 = 0.7;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.7;
gx2_2 = 0.7;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0.8 - 2*sqrt(3)/3*0.2;
gx2_1 = 0.8;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.8;
gx2_2 = 0.8;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0.9 - 2*sqrt(3)/3*0.1;
gx2_1 = 0.9;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.9;
gx2_2 = 0.9;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.9;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.1 - 2*sqrt(3)/3*0.9;
gx2_2 = 0.1;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.8;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.2 - 2*sqrt(3)/3*0.8;
gx2_2 = 0.2;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.7;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.3 - 2*sqrt(3)/3*0.7;
gx2_2 = 0.3;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.6;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.4 - 2*sqrt(3)/3*0.6;
gx2_2 = 0.4;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.5;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.5 - 2*sqrt(3)/3*0.5;
gx2_2 = 0.5;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.4;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.6 - 2*sqrt(3)/3*0.4;
gx2_2 = 0.6;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.3;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.7 - 2*sqrt(3)/3*0.3;
gx2_2 = 0.7;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.2;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.8 - 2*sqrt(3)/3*0.2;
gx2_2 = 0.8;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.1;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.9 - 2*sqrt(3)/3*0.1;
gx2_2 = 0.9;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.9;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.9 - 2*sqrt(3)/3*0;
gx2_2 = 0.9;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.8;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.8 - 2*sqrt(3)/3*0;
gx2_2 = 0.8;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.7;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.7 - 2*sqrt(3)/3*0;
gx2_2 = 0.7;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.6;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.6 - 2*sqrt(3)/3*0;
gx2_2 = 0.6;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.5;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.5 - 2*sqrt(3)/3*0;
gx2_2 = 0.5;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.4;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.4 - 2*sqrt(3)/3*0;
gx2_2 = 0.4;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.3;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.3 - 2*sqrt(3)/3*0;
gx2_2 = 0.3;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.2;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.2 - 2*sqrt(3)/3*0;
gx2_2 = 0.2;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')


gx1_1 = 2*sqrt(3)/3 - sqrt(3)/3*0 - 2*sqrt(3)/3*0.1;
gx2_1 = 0;
gx1_2 = 2*sqrt(3)/3 - sqrt(3)/3*0.1 - 2*sqrt(3)/3*0;
gx2_2 = 0.1;
X = [gx1_1 gx1_2];
Y = [gx2_1 gx2_2];
plot(X,Y,':w',X,Y,':w')

