function plot_data()
% Plots all the paper figures.

disp('Preparing to plot 39 figures, which may take a couple of minutes...')

%% Set text interpreter to 'tex' - default restored on cleanup
tick_interp = get(groot, 'defaultAxesTickLabelInterpreter');
leg_interp = get(groot, 'defaultLegendInterpreter');
restore_interp = onCleanup(@() set(groot, ...
    'defaultAxesTickLabelInterpreter', tick_interp, ...
    'defaultLegendInterpreter', leg_interp));
set(groot, 'defaultAxesTickLabelInterpreter', 'tex', ...
    'defaultLegendInterpreter', 'tex');

%% Fig. 2
plot_color_fig('fig_2.dat',false,true)

load('cmap_fig2.dat','-ascii','cmap_fig2')
colormap(cmap_fig2)
caxis([14 1e3])

xlim([0 2.15])
ylim([0 4.499])
xlabel('Wave vector')
ylabel('Energy (meV)')
c=colorbar(gca,'Location','northoutside');
c.Label.String='Intensity (arb. units)';

xticks([0 0.72 1.08 1.44 2.15])
xticklabels({'\Gamma''','K''_1','M''_2','K''_2','\Gamma'''})

drawnow

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

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 500])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

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

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 500])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

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

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 500])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

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

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 500])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

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

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 500])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

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

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 500])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

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

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 500])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

%% Fig. 3(h)
plot_color_fig('fig_3h.dat',false,true)

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 500])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

%% Fig. 3(i)
plot_color_fig('fig_3i.dat',false,true)

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 126])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

%% Fig. 3(j)
plot_color_fig('fig_3j.dat',false,true)

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 126])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

%% Fig. 3(k)
plot_color_fig('fig_3k.dat',false,true)

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 126])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

%% Fig. 3(l)
plot_color_fig('fig_3l.dat',false,true)

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 126])

c=colorbar;
c.Label.String='Intensity (arb. units)';

axis([0 1 -0.25 1.25])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

%% Fig. 4
data=importdata('fig_4.xy');

figure
h1=plot(data(:,1),data(:,5),'m--'); hold on
h2=plot(data(:,1),data(:,6),'m-');
h3=plot(data(:,1),data(:,7),'g--');
plot(data(:,1),data(:,8),'g-')
h4=plot(data(:,1),data(:,9),'b--');
plot(data(:,1),data(:,10),'b-'), hold off

axis([0 1.44 0 3.7])
xlabel('Wave vector')
ylabel('Energy (meV)')
legend([h1 h2 h3 h4],'\omega_{LSWT}','\omega','\Omega_{K/2}','\Omega_M')

xticks([0 0.72 1.08 1.44])
xticklabels({'\Gamma' 'K_1' 'M_2' 'K_2'})

drawnow

%% Fig. 5(a)
plot_color_fig('fig_5a.dat',false,false)

load('cmap_fig5a.dat','-ascii','cmap_fig5a')
colormap(cmap_fig5a)
caxis([0 8/3])

c=colorbar(gca,'Location','northoutside');
c.Label.String='\omega_{LSWT} (meV)';

axis([-0.417 0.417 -0.65 0.65])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

drawnow

%% Fig. 5(b)
plot_color_fig('fig_5b.dat',false,false)

load('cmap_fig5b.dat','-ascii','cmap_fig5b')
colormap(cmap_fig5b)
caxis([0 1.6473])

c=colorbar(gca,'Location','northoutside');
c.Label.String='\omega (meV)';

axis([-0.417 0.417 -0.65 0.65])
set(gca,'DataAspectRatio',[0.46434 0.80425 0.37868])

title('')
xlabel('(h, h)')
ylabel('(-k, k)')

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

load('cmap_fig2.dat','-ascii','cmap_fig2')
colormap(cmap_fig2)
caxis([4 1e3])

xlabel('Wave vector')
ylabel('Energy (meV)')              
c=colorbar;
c.Label.String='Intensity (arb. units)';

xticks([0 0.72 1.08 1.44 2.15])
xticklabels({'\Gamma''','K''_1','M''_2','K''_2','\Gamma'''})

data=importdata('fig_6.xy');

hold on
h1=plot(data(:,1),data(:,5),'m-','LineWidth',2);
h2=plot(data(:,1),data(:,6),'m--','LineWidth',2);
h3=plot(data(:,1),data(:,7),'m:','LineWidth',2);

data=importdata('fig_6.xye');

plot(data(:,1),data(:,5),'Color','k','Marker','o','MarkerFaceColor','w','LineStyle','none')
xlim([0 2.15])
hold off

legend([h1 h2 h3],{'\omega(k)' '\omega(k-Q)' '\omega(k+Q)'},'Location','SouthEast')

drawnow

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

load('cmap_fig2.dat','-ascii','cmap_fig2')
colormap(cmap_fig2)
caxis([4 1e3])

xlabel('Wave vector')
ylabel('Energy (meV)')              
c=colorbar;
c.Label.String='Intensity (arb. units)';

xticks([0 0.72 1.08 1.44 2.15])
xticklabels({'\Gamma''','K''_1','M''_2','K''_2','\Gamma'''})

data=importdata('fig_6.xy');

hold on
h1=plot(data(:,1),data(:,5),'m-','LineWidth',2);
h2=plot(data(:,1),data(:,6),'m--','LineWidth',2);
h3=plot(data(:,1),data(:,7),'m:','LineWidth',2);

data=importdata('fig_6.xye');

plot(data(:,1),data(:,5),'Color','k','Marker','o','MarkerFaceColor','w','LineStyle','none')
xlim([0 2.15])
hold off

legend([h1 h2 h3],{'\omega(k)' '\omega(k-Q)' '\omega(k+Q)'},'Location','SouthEast')

drawnow

%% Fig. 7(a)
plot_color_fig('fig_7a.dat',false,true)

load('cmap_fig2.dat','-ascii','cmap_fig2')
colormap(cmap_fig2)
caxis([4 1e3])

xlabel('Wave vector')
ylabel('Energy (meV)')
c=colorbar;
c.Label.String='Intensity (arb. units)';

xticks([0 0.36 0.98 1.70 2.13 2.49 3.11])
xticklabels({'K_2','M_2','\Gamma','K_2','K''_2','M''_2','\Gamma'''})

data=importdata('fig_7.xy');

hold on
h1=plot(data(:,1),data(:,5),'m-','LineWidth',2);
h2=plot(data(:,1),data(:,6),'m--','LineWidth',2);
h3=plot(data(:,1),data(:,7),'m:','LineWidth',2);
hold off

xlim([0 3.11])

legend([h1 h2 h3],{'\omega(k)' '\omega(k-Q)' '\omega(k+Q)'},'Location','SouthEast')

drawnow

%% Fig. 7(b)
plot_color_fig('fig_7b.dat',false,false)

load('cmap_fig2.dat','-ascii','cmap_fig2')
colormap(cmap_fig2)
caxis([4 1e3])

xlabel('Wave vector')
ylabel('Energy (meV)')
c=colorbar;
c.Label.String='Intensity (arb. units)';

xticks([0 0.36 0.98 1.70 2.13 2.49 3.11])
xticklabels({'K_2','M_2','\Gamma','K_2','K''_2','M''_2','\Gamma'''})

data=importdata('fig_7.xy');

hold on
h1=plot(data(:,1),data(:,5),'m-','LineWidth',2);
h2=plot(data(:,1),data(:,6),'m--','LineWidth',2);
h3=plot(data(:,1),data(:,7),'m:','LineWidth',2);
hold off

xlim([0 3.11])

legend([h1 h2 h3],{'\omega(k)' '\omega(k-Q)' '\omega(k+Q)'},'Location','SouthEast')

drawnow

%% Fig. 8(a)
data=importdata('fig_8a.xy');

figure, hold on
fill([flip(data(:,1));data(:,1)],[flip(data(:,8));data(:,9)],0.7*ones(1,3),'LineStyle','none')
plot(data(:,1),data(:,5),'m-','LineWidth',2)
plot(data(:,1),data(:,6),'g--','LineWidth',2)
plot(data(:,1),data(:,7),'c:','LineWidth',2)
hold off

box on
axis([0 1.70 0 3])
xlabel('Wave vector')
ylabel('Energy (meV)')
legend('Two-magnon continuum','\omega(k)','\omega(k-Q)','\omega(k+Q)','Location','SouthEast')

xticks([0 0.72 1.08 1.70])
xticklabels({'\Gamma','K','M','\Gamma'})
set(gca,'Layer','Top')

drawnow

%% Fig. 8(b)
data=importdata('fig_8b.xy');

figure, hold on
fill([flip(data(:,1));data(:,1)],[flip(data(:,8));data(:,9)],0.7*ones(1,3),'LineStyle','none')
plot(data(:,1),data(:,5),'m-','LineWidth',2)
plot(data(:,1),data(:,6),'g--','LineWidth',2)
plot(data(:,1),data(:,7),'c:','LineWidth',2)
hold off

box on
axis([0 1.7 0 3])
xlabel('Wave vector')
ylabel('Energy (meV)')
legend('Two-magnon continuum','\omega(k)','\omega(k-Q)','\omega(k+Q)','Location','SouthEast')

xticks([0 0.72 1.08 1.70])
xticklabels({'\Gamma','K','M','\Gamma'})
set(gca,'Layer','Top')

drawnow

%% Fig. 8(c)
data=importdata('fig_8c.xy');

figure, hold on
fill([flip(data(:,1));data(:,1)],[flip(data(:,8));data(:,9)],0.7*ones(1,3),'LineStyle','none')
plot(data(:,1),data(:,5),'m-','LineWidth',2)
plot(data(:,1),data(:,6),'g--','LineWidth',2)
plot(data(:,1),data(:,7),'c:','LineWidth',2)
hold off

box on
axis([0 1.70 0 3])
xlabel('Wave vector')
ylabel('Energy (meV)')
legend('Two-magnon continuum','\omega(k)','\omega(k-Q)','\omega(k+Q)','Location','NorthEast')

xticks([0 0.72 1.08 1.70])
xticklabels({'\Gamma','K','M','\Gamma'})
set(gca,'Layer','Top')

drawnow

%% Fig. 8(d)
data=importdata('fig_8d.xy');

figure, hold on
fill([flip(data(:,1));data(:,1)],[flip(data(:,8));data(:,9)],0.7*ones(1,3),'LineStyle','none')
plot(data(:,1),data(:,5),'m-','LineWidth',2)
plot(data(:,1),data(:,6),'g--','LineWidth',2)
plot(data(:,1),data(:,7),'c:','LineWidth',2)
hold off

box on
axis([0 1.70 0 3])
xlabel('Wave vector')
ylabel('Energy (meV)')
legend('Two-magnon continuum','\omega(k)','\omega(k-Q)','\omega(k+Q)','Location','NorthEast')

xticks([0 0.72 1.08 1.70])
xticklabels({'\Gamma','K','M','\Gamma'})
set(gca,'Layer','Top')

drawnow

%% Fig. 9(a)
plot_color_fig('fig_9a.dat',false,true)

axis([0 1 0 5])
caxis([0 577])
ylabel('Energy (meV)')
xlabel('(h, h, 0)')

load('cmap_fig9.dat','-ascii','cmap_fig9')
colormap(cmap_fig9)

c=colorbar(gca,'Location','northoutside');
c.Label.String='Intensity (arb. units)';

pbaspect([0.98101 1 0.98101])

drawnow

%% Fig. 9(b)
plot_color_fig('fig_9b.dat',false,false)

axis([0 1 0 5])
caxis([0 577])
ylabel('Energy (meV)')
xlabel('(h, h, 0)')

load('cmap_fig9.dat','-ascii','cmap_fig9')
colormap(cmap_fig9)

c=colorbar(gca,'Location','northoutside');
c.Label.String='Intensity (arb. units)';

pbaspect([0.98101 1 0.98101])

drawnow

%% Fig. 10(a)
data=importdata('fig_10a.xye');

figure, hold on
idx=~isnan(data(:,8));
h2m_area=fill(data(idx,7),data(idx,8),[0.7 0.7 0.7],'LineStyle','none'); 
idx=~isnan(data(:,18));
h2m_sim_area=fill(data(idx,17),data(idx,18),[1 0.7 0.7],'LineStyle','none'); 
ebar_pars={'CapSize',0,'LineStyle','none','Color','k'};
hint_loE=errorbar(data(:,1),data(:,2),data(:,3),'Marker','o','MarkerFaceColor','k',ebar_pars{:});
hint_hiE=errorbar(data(:,4),data(:,5),data(:,6),'Marker','d','MarkerFaceColor','k',ebar_pars{:});
hbkg_loE=errorbar(data(:,9),data(:,10),data(:,11),'Marker','o','MarkerFaceColor','w',ebar_pars{:});
hbkg_hiE=errorbar(data(:,12),data(:,13),data(:,14),'Marker','d','MarkerFaceColor','w',ebar_pars{:});
h1m_sim=plot(data(:,15),data(:,16),'r-','LineWidth',2);
h2m_sim=plot(data(:,17),data(:,18),'r--','LineWidth',2);
hold off

legend([hint_loE;hint_hiE;hbkg_loE;hbkg_hiE;h2m_area;h1m_sim;h2m_sim;h2m_sim_area],...
    {'Total intensity [M (-1/2, 3/2, 0)]' '' 'Background (-1/4, 1, 0)' ...
    'Background (-1, 2, 0)' 'Estimated continuum scattering' ...
    'One magnon calculation' 'Two magnon calculation' ''},...
    'Box','off');

box on
axis([0.5 6 -10 225])
xlabel('Energy (meV)')
ylabel('Intensity (arb. units)')

drawnow

%% Fig. 10(b)
data=importdata('fig_10b.xye');

figure, hold on
idx=~isnan(data(:,8));
h2m_area=fill(data(idx,7),data(idx,8),[0.7 0.7 0.7],'LineStyle','none'); 
h2m_sim_area=fill(data(:,20),data(:,21),[1 0.7 0.7],'LineStyle','none'); 
ebar_pars={'CapSize',0,'LineStyle','none','Color','k'};
hdata_loE=errorbar(data(:,1),data(:,2),data(:,3),'Marker','o','MarkerFaceColor','k',ebar_pars{:});
hdata_hiE=errorbar(data(:,4),data(:,5),data(:,6),'Marker','d','MarkerFaceColor','k',ebar_pars{:});
hbkg1=errorbar(data(:,9),data(:,10),data(:,11),'Marker','s','MarkerFaceColor','w',ebar_pars{:});
hbkg2=errorbar(data(:,12),data(:,13),data(:,14),'Marker','d','MarkerFaceColor','w',ebar_pars{:});
hbkg2_loE=errorbar(data(:,15),data(:,16),data(:,17),'Marker','o','MarkerFaceColor','w',ebar_pars{:});
h1m_sim=plot(data(:,18),data(:,19),'r-','LineWidth',2);
h2m_sim=plot(data(:,20),data(:,21),'r--','LineWidth',2);
idx=data(:,7)<1.8;
h2m=plot(data(idx,7),data(idx,8),'k--','LineWidth',2);
hold off

legend([hdata_loE hdata_hiE hbkg1 hbkg2_loE hbkg2 h2m h2m_area h1m_sim h2m_sim h2m_sim_area],...
    {'Total intensity [K''_2 (-1/3, 2/3, 1)]' '' ...
    'Background [M_2 (0, 1/2, 0)]' 'Background [\Gamma (0, 1, 0)]' '' ...
    'Estimated continuum scattering' '' 'One magnon calculation' 'Two magnon calculation' ''},...
    'Box','off','Location','NorthEast')

box on
axis([0.14 5 -20 275])
xlabel('Energy (meV)')
ylabel('Intensity (arb. units)')

drawnow

%% Fig. 11
hCont=plot_color_fig('fig_11.dat',true,false);

cdata=hCont.ZData;
cdata(cdata<0)=NaN;
log_cdata=log10(cdata);
hCont.LevelList=-2+log10(2):1/3:-1/3+log10(2);
hCont.ZData=log_cdata;
hBar=colorbar(gca,'Location','NorthOutside');
ticks=[0.02 0.2 2];
caxis(log10([0.02 2]))
hBar.Ticks=log10(ticks);
hBar.TickLabels=arrayfun(@(x) {num2str(x)},ticks);
hBar.Label.Interpreter='LaTeX';
hBar.Label.String='$\widetilde{D}(\mathbf{k},\omega)$ (meV$^{-1}$Co$^{-1}$)';

load('cmap_fig11.dat','-ascii','cmap_fig11')
colormap(cmap_fig11)

data=importdata('fig_11.xy');

hold on
h2=plot(data(:,1),data(:,5),'g--','LineWidth',2);
h3=plot(data(:,1),data(:,6),'c:','LineWidth',2);
hG=plot(data(:,7),data(:,11),'r.','MarkerSize',2);
hK=plot(data(:,12),data(:,16),'b.','MarkerSize',2);
hold off

legend([h2 h3 hG hK],{'\omega(k-Q)','\omega(k+Q)',...
    '\omega_{2M}(k_1 \sim k_\Gamma)','\omega_{2M}(k_1 \sim k_{K''})'},... 
    'Box','off','Location','South')

xlabel('Wave vector')
ylabel('Energy (meV)')

xlim([0 3.11])
xticks([0 0.36 0.98 1.70 2.13 2.49 3.11])
xticklabels({'K_1','M_2','\Gamma','K_1','K''_1','M''_2','\Gamma'''})

pbaspect([0.98101 1 0.98101])

drawnow

%% Fig. 12(a)
plot_color_fig('fig_12a.dat',false,true)

axis([0.1 0.9 0 3])
xlabel('(h, h, 0)')
ylabel('Energy (meV)')

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 800])

c=colorbar;
c.Label.String='Intensity (arb. units)';

drawnow

%% Fig. 12(b)
plot_color_fig('fig_12b.dat',false,true)

axis([0.1 0.9 0 3])
xlabel('(h, h, 0)')
ylabel('Energy (meV)')

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 800])

c=colorbar;
c.Label.String='Intensity (arb. units)';

drawnow

%% Fig. 12(c)
plot_color_fig('fig_12c.dat',false,true)

axis([0.1 0.9 0 3])
xlabel('(h, h, 0)')
ylabel('Energy (meV)')

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 800])

c=colorbar;
c.Label.String='Intensity (arb. units)';

drawnow

%% Fig. 12(d)
plot_color_fig('fig_12d.dat',false,true)

axis([0.1 0.9 0 3])
xlabel('(h, h, 0)')
ylabel('Energy (meV)')

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 800])

c=colorbar;
c.Label.String='Intensity (arb. units)';

drawnow

%% Fig. 12(e)
plot_color_fig('fig_12e.dat',false,false)

axis([0.1 0.9 0 3])
xlabel('(h, h, 0)')
ylabel('Energy (meV)')

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 800])

c=colorbar;
c.Label.String='Intensity (arb. units)';

drawnow

%% Fig. 12(f)
plot_color_fig('fig_12f.dat',false,false)

axis([0.1 0.9 0 3])
xlabel('(h, h, 0)')
ylabel('Energy (meV)')

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 800])

c=colorbar;
c.Label.String='Intensity (arb. units)';

drawnow

%% Fig. 12(g)
plot_color_fig('fig_12g.dat',false,false)

axis([0.1 0.9 0 3])
xlabel('(h, h, 0)')
ylabel('Energy (meV)')

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 800])

c=colorbar;
c.Label.String='Intensity (arb. units)';

drawnow

%% Fig. 12(h)
plot_color_fig('fig_12h.dat',false,false)

axis([0.1 0.9 0 3])
xlabel('(h, h, 0)')
ylabel('Energy (meV)')

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 800])

c=colorbar;
c.Label.String='Intensity (arb. units)';

drawnow

%% Fig. 13(a)
plot_color_fig('fig_13a.dat',false,false)

axis([0.1 0.9 0 3])
xlabel('(h, h, 0)')
ylabel('Energy (meV)')

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 800])

c=colorbar;
c.Label.String='Intensity (arb. units)';

data=importdata('fig_13a.xy');

hold on
h1=plot(data(:,1),data(:,5),'m-','LineWidth',2);
h2=plot(data(:,1),data(:,6),'m--','LineWidth',2);
h3=plot(data(:,1),data(:,7),'m:','LineWidth',2);
hold off

legend([h1 h2 h3],'\omega(k)','\omega(k-Q)','\omega(k+Q)','Location','South')

drawnow

%% Fig. 13(b)
plot_color_fig('fig_13b.dat',false,false)

axis([0.1 0.9 0 3])
xlabel('(h, h, 0)')
ylabel('Energy (meV)')

load('cmap_fig3.dat','-ascii','cmap_fig3')
colormap(cmap_fig3)
caxis([0 800])

c=colorbar;
c.Label.String='Intensity (arb. units)';

data=importdata('fig_13b.xy');

hold on
h1=plot(data(:,1),data(:,5),'m-','LineWidth',2);
h2=plot(data(:,1),data(:,6),'m--','LineWidth',2);
h3=plot(data(:,1),data(:,7),'m:','LineWidth',2);
hold off

legend([h1 h2 h3],'\omega(k)','\omega(k-Q)','\omega(k+Q)','Location','South')

drawnow

end

%% Subroutines to plot data
function hout=plot_color_fig(fdata,cont,isErr,varargin)
% Plots color data in .dat file, with cont flag specifying if contour plot.
% Additionally evaluates code strings specified as optional arguments.
fdata=importdata(fdata,' ',1);

sz=eval(fdata.textdata{1});

% Reshape to 2D array for plotting
X=reshape(fdata.data(:,1),sz);
Y=reshape(fdata.data(:,end-isErr-1),sz);
C=reshape(fdata.data(:,end-isErr),sz);

[X,Y,C]=cen2edge(X,Y,C);  % convert bin centers to edges for pcolor

figure
if cont  % if contour plot
    [~,h]=contourf(X,Y,C,'LineWidth',1);
else  
    h=pcolor(X,Y,C);
end
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

if nargout>0
    hout=h;
end
end

function [Xedge,Yedge,varargout]=cen2edge(Xcen,Ycen,varargin)
% Converts bin edges to centers
diffX=diff(Xcen,1,2);
Xedge=Xcen(:,1:end-1)-diffX/2;
Xedge(:,end+1)=Xcen(:,end)+diffX(:,end)/2;

diffY=diff(Ycen,1,1);
Yedge=Ycen(1:end-1,:)-diffY/2;
Yedge(end+1,:)=Ycen(end,:)+diffY(end,:)/2;

varargout=cell(size(varargin));
for i=1:numel(varargin)
    array=reshape(varargin{i},size(Xcen));
    varargout{i}=array;
end
end
