%% Fig. 3(a)
plot_color_fig('fig_3a.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([3 8])
axis([0 3 0 20])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (mb sr^{-1} meV^{-1} Cu^{-1})';
c.TickDirection='out';

%% Fig. 3(b)
plot_color_fig('fig_3b.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([3 8])
axis([0 3 0 20])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (mb sr^{-1} meV^{-1} Cu^{-1})';
c.TickDirection='out';

%% Fig. 3(c)
plot_color_fig('fig_3c.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([3 8])
axis([0 3 0 20])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (mb sr^{-1} meV^{-1} Cu^{-1})';
c.TickDirection='out';

%% Fig. 3(d)
plot_color_fig('fig_3d.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([3 8])
axis([0 3 0 20])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (mb sr^{-1} meV^{-1} Cu^{-1})';
c.TickDirection='out';

%% Fig. 3(e)
plot_color_fig('fig_3e.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([3 8])
axis([0 3 0 20])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (mb sr^{-1} meV^{-1} Cu^{-1})';
c.TickDirection='out';

%% Fig. 3(f)
plot_color_fig('fig_3f.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([3 8])
axis([0 3 0 20])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (mb sr^{-1} meV^{-1} Cu^{-1})';
c.TickDirection='out';

%% Fig. 3(g)
plot_color_fig('fig_3g.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([3 8])
axis([0 3 0 20])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (mb sr^{-1} meV^{-1} Cu^{-1})';
c.TickDirection='out';

%% Fig. 3(h)
data=importdata('fig_3h.xye');

figure, hold on
% INS data (scan C)
errorbar(data(:,1),data(:,2),data(:,3),...
    'Color','k','LineStyle','none','Marker','o','MarkerFaceColor','k')
% LSWT (J'=0.35J)
plot(data(:,1),data(:,5),'Color','b','LineWidth',2,'LineStyle','--')
% LSWT (J'=0.55J)
plot(data(:,1),data(:,6),'Color','g','LineWidth',2,'LineStyle',':')
% LSWT (J'=0.45J, scan C)
plot(data(:,1),data(:,4),'Color','r','LineWidth',2,'LineStyle','-')
% INS data (scan D)
errorbar(data(:,7),data(:,8),data(:,9),...
    'Color','k','LineStyle','none','Marker','o','MarkerFaceColor','w')
% LSWT (J'=0.45J, scan D)
plot(data(:,7),data(:,10),'Color','r','LineWidth',2,'LineStyle','-')
hold off

% Instrumental resolution (FWHM)
line([14.7 16.4], [5.8 5.8],'Color','k')
line([14.7 14.7], [5.7 5.9],'Color','k')
line([16.4 16.4], [5.7 5.9],'Color','k')
% QMC broadening (FWHM)
line([13.6 17.5], [5.5 5.5],'Color','k', 'LineStyle', '--')
line([13.6 13.6], [5.4 5.6],'Color','k')
line([17.5 17.5], [5.4 5.6],'Color','k')

box on
title('')
axis([7 18 2.5 9.5])
legend({'INS data (scan C)' 'LSWT (J'' = 0.35J)' 'LSWT (J'' = 0.55J)'...
    'LSWT (J'' = 0.45J)' 'INS data (scan D)'},...
    'FontName','Times','Box','off','Location','NorthWest')
xlabel('Energy (meV)');
ylabel('Intensity (mb sr^{-1} meV^{-1} Cu^{-1})');
set(gca,'TickDir','out')

%% Fig. 3(i)
data=importdata('fig_3i.xye');

figure, hold on
% INS data (scan B)
errorbar(data(:,1),data(:,2),data(:,3),...
    'Color','k','LineStyle','none','Marker','o','MarkerFaceColor','k')
% LSWT (J'=0.45J)
plot(data(:,4),data(:,5),'Color','r','LineWidth',2,'LineStyle','-')
% Magnetic form factor
plot(data(:,4),data(:,6),'k--','LineWidth',2)
% INS data (scan A)
errorbar(data(:,7),data(:,8),data(:,9),...
    'Color','k','LineStyle','none','Marker','o','MarkerFaceColor','w')
% Zero intensity
line([0 10],[0 0],'LineStyle',':','Color','k','LineWidth',2)

box on
title('')
axis([1.3 6 -0.5 4])
legend({'INS data (scan B)' 'LSWT (J'' = 0.45J)' 'Magnetic form factor' 'INS data (scan A)'},...
    'FontName','Times','Box','off','Location','NorthEast')
xlabel(['|Q| (' char(197) '^{-1})']);
ylabel('Intensity (mb sr^{-1} meV^{-1} Cu^{-1})');
set(gca,'TickDir','out')

%% Fig. 4(a)
plot_color_fig('fig_4a.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([0 300])
axis([0 2.5 0 2.3])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (arb. units)';
c.TickDirection='out';

%% Fig. 4(b)
plot_color_fig('fig_4b.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([0 300])
axis([0 2.5 0 2.3])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (arb. units)';
c.TickDirection='out';

%% Fig. 4(c)
plot_color_fig('fig_4c.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([0 300])
axis([0 1.9 0 1.2])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (arb. units)';
c.TickDirection='out';

%% Fig. 4(d)
plot_color_fig('fig_4d.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([0 300])
axis([0 2.5 0 2.3])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (arb. units)';
c.TickDirection='out';

%% Fig. 4(e)
data=importdata('fig_4e.xye');

figure, hold on
% INS data (scan G)
errorbar(data(:,1),data(:,2),data(:,3),...
    'Color','k','LineStyle','none','Marker','o','MarkerFaceColor','k')
% INS data (scan H)
errorbar(data(:,4),data(:,5),data(:,6),...
    'Color','k','LineStyle','none','Marker','s','MarkerFaceColor','w')
% INS data (scan F)
errorbar(data(:,7),data(:,8),data(:,9),...
    'Color','k','LineStyle','none','Marker','o','MarkerFaceColor','w')
% LSWT (3D)
plot(data(:,10),data(:,11),'Color','r','LineWidth',2,'LineStyle','-')
% LSWT (gap)
plot(data(:,10),data(:,12),'Color','b','LineWidth',2,'LineStyle','--')
hold off

box on
axis([0.3 2.3 0 300])
legend({'INS signal (scan G)' 'INS background (scan H)' ...
    'INS background (scan F)' 'LSWT model (3D)' 'LSWT model (gap)'},...
    'FontName','Times','Box','off');
xlabel('Energy (meV)');
ylabel('Intensity (arb. units)');
set(gca,'TickDir','out')

%% Fig. 4(f)
data=importdata('fig_4f.xye');

figure, hold on
% INS data (scan E)
errorbar(data(:,1),data(:,2),data(:,3),'Color','k','LineStyle','none','Marker','o','MarkerFaceColor','k','MarkerSize',8)
% LSWT (3D)
plot(data(:,1),data(:,4),'Color','r','LineWidth',2,'LineStyle','-')
hold off

box on
axis([0.5 2.1 0 300])
legend({'INS signal (scan E)' 'LSWT model (3D)'},'FontName','Times','Box','off')
xlabel(['|Q| (' char(197) '^{-1})']);
ylabel('Intensity (arb. units)');
set(gca,'TickDir','out')

%% Fig. 5(a)
% Plot S(q,w)
plot_color_fig('fig_5a.dat')

load('cmap_Sqw.dat','-ascii')
colormap(cmap_Sqw)
caxis([0 8])
c=colorbar;
c.Label.String='S^{xx}(Q,\omega) (arb. units)';
c.TickDirection='out';

% Plot dispersions
data=importdata('fig_5a.xy');

hold on
h=plot(data(:,1),data(:,5),'m--',data(:,1),data(:,6),'m-'); hold off
set(findobj(gcf,'Type','Line'),'LineWidth',2)

% Ticks and dividing lines at high symmetry points
ticks=[0 0.6538 1.202 1.8558 2.4039];
line([ticks(2) ticks(2)],[0 25],'Color','k','LineStyle','--')
line([ticks(3) ticks(3)],[0 25],'Color','k','LineStyle','--')
line([ticks(4) ticks(4)],[0 25],'Color','k','LineStyle','--')
xticks(ticks)

box on
axis([ticks(1) ticks(end) 0 25])
xticklabels({'(0,1/2,0)' '(0,1/2,1)' '(0,1,1)' '(0,1,0)' '(0,1/2,0)'})
ylabel('Energy (meV)')
xlabel('Wavevector (r.l.u.)')
set(gca,'TickDir','out','Layer', 'Top') 
legend(h,{[char(295) char(969) '^+'] [char(295) char(969) '^-']},'Box','off')

%% Fig. 5(b)
% Plot S(q,w)
plot_color_fig('fig_5b.dat')

load('cmap_Sqw.dat','-ascii')
colormap(cmap_Sqw)
caxis([0 8])
c=colorbar;
c.Label.String='S^{xx}(Q,\omega) (arb. units)';
c.TickDirection='out';

% Ticks and dividing lines at high symmetry points
ticks=[0 0.6538 1.202 1.8558 2.4039];
line([ticks(2) ticks(2)],[0 25],'Color','k','LineStyle','--')
line([ticks(3) ticks(3)],[0 25],'Color','k','LineStyle','--')
line([ticks(4) ticks(4)],[0 25],'Color','k','LineStyle','--')
xticks(ticks)

box on
axis([ticks(1) ticks(end) 0 25])
xticklabels({'(0,1/2,0)' '(0,1/2,1)' '(0,1,1)' '(0,1,0)' '(0,1/2,0)'})
ylabel('Energy (meV)')
xlabel('Wavevector (r.l.u.)')
set(gca,'TickDir','out','Layer', 'Top') 

%% Fig. 6(a)
plot_color_fig('fig_6a.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([0 250])
axis([0 3 0 18])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (arb. units)';
c.TickDirection='out';

%% Fig. 6(b)
plot_color_fig('fig_6b.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([0 250])
axis([0 3 0 18])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (arb. units)';
c.TickDirection='out';

%% Fig. 6(c)
data=importdata('fig_6c.xye');

figure, hold on
% INS data
errorbar(data(:,1),data(:,2),data(:,3),'Color','k','LineStyle','none','Marker','o','MarkerFaceColor','k')
% LSWT
plot(data(:,1),data(:,4),'Color','r','LineWidth',2,'LineStyle','-')
hold off

box on
axis([4 17 56 232])
legend({'INS data' 'LSWT model'},...
    'FontName','Times','Box','off','Location','NorthWest')
xlabel('Energy (meV)');
ylabel('Intensity (arb. units)');
set(gca,'TickDir','out')

%% Fig. 7
data=importdata('fig_7.xy');

figure, hold on
% Two magnon continuum
fill([data(:,1);data(end:-1:1,1)],[data(:,7);data(end:-1:1,8)],0.7*ones(1,3))
% Dispersions
plot(data(:,1),data(:,5),'r-',data(:,1),data(:,6),'b-','LineWidth',2)
hold off

% Ticks and dividing lines at high symmetry points
ticks=[data(1,1) data(101,1) data(201,1) data(301,1) data(401,1)];
line([ticks(2) ticks(2)],[0 35],'Color','k','LineStyle','--')
line([ticks(3) ticks(3)],[0 35],'Color','k','LineStyle','--')
line([ticks(4) ticks(4)],[0 35],'Color','k','LineStyle','--')
xticks(ticks)

box on
xlim([ticks(1) ticks(end)])
xticklabels({'(0,1/2,0)','(0,1/2,1)','(0,1,1)','(0,1,0)','(0,1/2,0)'})
xlabel('Wavevector (r.l.u.)')
ylabel('Energy (meV)')
legend({[char(295) char(969) '^' char(177) '+' char(295) char(969) '^' char(177)] ...
    [char(295) char(969) '^+'] [char(295) char(969) '^-']},...
    'FontName','Times','Box','off','Location','NorthEast')
set(gca,'TickDir','out','Layer','Top')

%% Fig. 8
plot_color_fig('fig_8.dat',true)

load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([3 8])
axis([0 3 0 20])
xlabel(['|Q| (' char(197) '^{-1})'])
ylabel('Energy (meV)')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='Intensity (mb sr^{-1} meV^{-1} Cu^{-1})';
c.TickDirection='out';

%% Fig. 9
data=importdata('fig_9.xy');

figure, hold on
% Experimental data (B||ab)
plot(data(:,1),data(:,2),'r','LineWidth',2)
% Experimental data (B||c)
plot(data(:,3),data(:,4),'b','LineWidth',2)
% MFT calculation
plot(data(:,5),data(:,6),'--g','LineWidth',2)
% QMC calculation
plot(data(:,7),data(:,8),'--k','LineWidth',2)
hold off

box on
ylabel('Magnetization (<S_z> per Cu^{2+})')
xlabel('gB/2 (T)')
legend({'Experimental data (B||ab)' 'Experimental data (B||c)' ...
    'Mean field calculation' 'QMC calculation'},'Location','northwest','Box','off')
xlim([0 52]), ylim([0 0.1])
set(gca,'TickDir','out')

%% Fig. 10
% Plot Chi-squared
plot_color_fig('fig_10.dat',false,'contour(x,y,c,[1 1e3])')

data=importdata('fig_10.xy');
hold on
% Eq.(B2) line
plot(data(:,1),data(:,2),'w--','LineWidth',2)
% J'=Jrung line
plot(0.2:0.1:1.5,0.2:0.1:1.5,'w','LineWidth',2)
% Best fit parameters
plot([0.45 1],[1 0.45],'wo','LineWidth',2)
hold off

axis square
load('cmap_int.dat','-ascii')
colormap(cmap_int)
caxis([0 30])
xlabel('J_{rung}/J_{leg}')
ylabel('J''/J_{leg}')
set(gca,'TickDir','out','Layer', 'Top')
c=colorbar;
c.Label.String='\chi^2';
c.TickDirection='out';
set(findobj(gcf,'Type','Contour'),'LineColor','w','LineWidth',2)

%% Fig. 11(a)
% Plot S(q,w)
plot_color_fig('fig_11a.dat')

load('cmap_Sqw.dat','-ascii')
colormap(cmap_Sqw)
caxis([0 2])
c=colorbar;
c.Label.String='S^{xx}(Q,\omega) (arb. units)';
c.TickDirection='out';

% Ticks and dividing lines at high symmetry points
xmax=max(max(get(findobj(gcf,'Type','Surface'),'XData')));
ticks=[0 xmax/2 xmax];
line([ticks(2) ticks(2)],[0 25],'Color','k','LineStyle','--')
xticks(ticks)

box on
axis([0 xmax 0 4])
xticklabels({'0' '1/2' '1'})
ylabel('Energy/J')
xlabel('k (r.l.u.)')
set(gca,'TickDir','out','Layer', 'Top') 
text(0.1,3.8,'l = 1','FontName',get(gca,'FontName'),'FontSize',get(gca,'FontSize'))

%% Fig. 11(b)
% Plot S(q,w)
plot_color_fig('fig_11b.dat')

load('cmap_Sqw.dat','-ascii')
colormap(cmap_Sqw)
caxis([0 2])
c=colorbar;
c.Label.String='S^{xx}(Q,\omega) (arb. units)';
c.TickDirection='out';

% Ticks and dividing lines at high symmetry points
xmax=max(max(get(findobj(gcf,'Type','Surface'),'XData')));
ticks=[0 xmax/2 xmax];
line([ticks(2) ticks(2)],[0 25],'Color','k','LineStyle','--')
xticks(ticks)

box on
axis([0 xmax 0 4])
xticklabels({'0' '1/2' '1'})
ylabel('Energy/J')
xlabel('k (r.l.u.)')
set(gca,'TickDir','out','Layer', 'Top') 
text(0.1,3.8,'l = 0','FontName',get(gca,'FontName'),'FontSize',get(gca,'FontSize'))

%% Functions to plot data
function plot_color_fig(fdata,ms,varargin)
% Plots color data in .dat file, with ms flag specifying if it mslice data.
% Additionally evaluates code strings specified as optional arguments.
fdata=importdata(fdata,' ',1);
eval(['sz=',fdata.textdata{1},';'])

x=reshape(fdata.data(:,1),sz);
y=reshape(fdata.data(:,2),sz);
c=reshape(fdata.data(:,3),sz);

% Convert bin centers to bin edges for pcolor
if exist('ms', 'var') && ms   
    diffx=diff(x,1,1);
    xx=x(1:end-1,:)+diffx/2;
    x1=x(1,:)-diffx(1,:)/2;
    xn=x(end,:)+diffx(end,:)/2;
    x=[x1;xx;xn];
    x(:,end+1)=x(:,end);
    
    diffy=diff(y,1,2);
    yy=y(:,1:end-1)+diffy/2;
    y1=y(:,1)-diffy(:,1)/2;
    yn=y(:,end)+diffy(:,end)/2;
    y=[y1 yy yn]; 
    y(end+1,:)=y(end,:);
    
    c(end+1,:)=NaN; 
    c(:,end+1)=NaN;
end

figure
pcolor(x,y,c)
shading flat

isNum=~isnan(c);
axis([min(x(isNum)) max(x(isNum)) min(y(isNum)) max(y(isNum))])

% Execute any strings entered as optional parameters
for i=1:length(varargin)
    hold on
    eval(varargin{i})
    hold off
end
end
