%% This Matlab script reproduces the figures / subfigures showing data in 
%% the paper Low energy phonons in single crystal ZrW2O8, by R. A. Ewings et al

%% Data files. Note paths should be the locations on disk of the .sqw files
proj2.u=[1,1,0]; proj2.v=[0,0,1]; proj2.uoffset=[0,0,0,0]; proj2.type='rrr';
sqw_file100='/mnt/ceph/instrument/LET/RBNumber/RB1690453/ZWO_ILL_Mar21/ZrWO_100K.sqw';
sqw_file400='/mnt/ceph/instrument/LET/RBNumber/RB1690453/ZWO_ILL_Mar21/ZrWO_400K.sqw';

proj.u=[1,0,0]; proj.v=[0,1,0]; proj.uoffset=[0,0,0,0]; proj.type='rrr';
sqw_file2='/mnt/ceph/instrument/LET/RBNumber/RB1690453/ZWO_ILL_Mar21/ZrW2O8_2K.sqw';
sqw_file300='/mnt/ceph/instrument/LET/RBNumber/RB1690453/ZWO_ILL_Mar21/ZrW2O8_295K.sqw';


%% Figure 2
bp1=[0 0 4]; bp2=[-2,-2,2]; bp3=[3,-3,0];

%Used this for bp1
rlp7=[0,0,0; 0.5,0,0; 0.5,0.5,0; 0,0.5,0; 0,0,0; 0,0,0.5; 0,0.5,0.5; 0,0.5,0; 0.5,0.5,-0.5; 0,0,0];

%For bp2
rlp8=[0,0,0; -0.5,0,0; -0.5,-0.5,0; 0,-0.5,0; 0,0,0; 0,0,0.5; 0,-0.5,0.5; 0,-0.5,0; -0.5,-0.5,0.5; 0,0,0];

%for bp3
rlp9=[0,0,0; 0.5,0,0; 0.5,0.5,0; 0,0.5,0; 0,0,0; 0,0,0.5; 0,0.5,0.5; 0,0.5,0; 0.5,0.5,-0.5; 0,0,0; -0.5,0.5,0];

rnewquick=bp1+rlp7;
w_100_1=spaghetti_plot_sqw(rnewquick,sqw_file100,'qbin',0.03,'labels',...
    {'G','X','M1','Y','G','Z','M2','Y','R','G'},...
        'clim',[0,0.0007],'ebin',[1,0.09,6],'qwidth',0.06,'noplot');

spaghetti_plot(bose(d2d(w_100_1)-7e-5,100));
lz 0 0.0002; title('');
keep_figure;

w_400_1=spaghetti_plot_sqw(rnewquick,sqw_file400,'qbin',0.03,'labels',...
    {'G','X','M1','Y','G','Z','M2','Y','R','G'},...
        'clim',[0,0.0007],'ebin',[1,0.09,6],'qwidth',0.06,'noplot');

spaghetti_plot(bose(d2d(w_400_1)-2e-4,400));
lz 0 0.0002; title('');
keep_figure;


%% Figure 3

rnewquick2=bp2+rlp8;
w_100_2=spaghetti_plot_sqw(rnewquick2,sqw_file100,'qbin',0.03,'labels',...
    {'G','X','M1','Y','G','Z','M2','Y','R','G'},...
        'clim',[0,0.0007],'ebin',[1,0.09,6],'qwidth',0.06,'noplot');

spaghetti_plot(bose(d2d(w_100_2)-7e-5,100));
lz 0 0.0002; title('');
keep_figure;

w_400_2=spaghetti_plot_sqw(rnewquick2,sqw_file400,'qbin',0.03,'labels',...
    {'G','X','M1','Y','G','Z','M2','Y','R','G'},...
        'clim',[0,0.0007],'ebin',[1,0.09,6],'qwidth',0.06);
spaghetti_plot(bose(d2d(w_400_2)-2e-4,400));
lz 0 0.0002; title('');
keep_figure;

%% Figure 4
rnewquick3=bp3+rlp9;
w_2_3=spaghetti_plot_sqw(rnewquick3,sqw_file2,'qbin',0.03,'labels',...
    {'G','X','M1','Y','G','Z','M2','Y','R','G','M3'},...
        'clim',[0,0.0007],'ebin',[1,0.15,7],'qwidth',0.1,'noplot');
spaghetti_plot(bose(d2d(w_2_3)-7e-7,2));
lz 0 3.5e-6; title('');
keep_figure;

w_300_3=spaghetti_plot_sqw(rnewquick3,sqw_file300,'qbin',0.03,'labels',...
    {'G','X','M1','Y','G','Z','M2','Y','R','G','M3'},...
        'clim',[0,0.0007],'ebin',[1,0.15,7],'qwidth',0.1,'noplot');
spaghetti_plot(bose(d2d(w_300_3)-3.6e-6,300));
lz 0 3.5e-6; title('');
keep_figure;

%% Figure 5

%Midpoint of Y-R trajectory 2K and 300K
proj3.u=[1,0,0]; proj3.v=[0,1,0]; proj3.type='rrr';

ctest=cut_sqw(sqw_file2,proj3,[3.15,3.35],[-2.6,-2.4],[-0.35,-0.15],[0,0.07,7],'-nopix');
ctest2=cut_sqw(sqw_file300,proj3,[3.15,3.35],[-2.6,-2.4],[-0.35,-0.15],[0,0.07,7],'-nopix');
acolor black
amark 8
plot(bose(ctest,2));
acolor red
pp((bose(ctest2,300))*2);
lx 2 6


%% Figure 6

%Midpoint of Z-M2 trajectory at 100K and 400K
crapola1=cut(bose(d2d(w_100_1(6))-7e-5,100),[0.1,0.4],[]);
crapola2=cut(bose(d2d(w_400_1(6))-2e-4,400),[0.1,0.4],[]);

amark(8)
acolor blue; plot(crapola1); acolor red; pp(crapola2*2);
lx 2 5
ly -1e-5 2.2e-4
set(gca,'FontSize',14);
xlabel('Energy (meV)')
ylabel('Intensity (arb. units)');
title('')
grid on
legend('100 K','400 K','Location','NorthWest')
legend('boxoff')
keep_figure;


%% Figure 7

%4 cuts at high symmetry points

% G point
proj.u=[1,0,0]; proj.v=[0,1,0]; proj.type='rrr'; proj.uoffset=[3,-3,0,0];
c(1)=cut_sqw(sqw_file2,proj,[-0.04,0.04],[-0.04,0.04],[-0.04,0.04],[2,0.124,7],'-nopix');
c(2)=cut_sqw(sqw_file300,proj,[-0.04,0.04],[-0.04,0.04],[-0.04,0.04],[2,0.124,7],'-nopix');
acolor blue
%plot(c(1))
plot(bose(c(1)-1e-6,2))
acolor red
%pp(c(2))
pp(2*bose(c(2)-5e-6,295))
lx 2 7
keep_figure;
title('');
set(gca,'FontSize',14);
ylabel('Intensity (arb. units)');
xlabel('Energy (meV)');
grid on
ll=legend('2 K','300 K','Location','NorthWest');
set(ll,'FontSize',14);
legend('boxoff');
text(2.5,4e-6,'\bf Q \rm = (3, -3, 0)','FontSize',14);

% Z point
proj.u=[1,1,0]; proj.v=[0,0,1]; proj.type='rrr'; proj.uoffset=[-2,-2,2.5,0];
cz(1)=cut_sqw(sqw_file100,proj,[-0.04,0.04],[-0.04,0.04],[-0.04,0.04],[1,0.12,6],'-nopix');
cz(2)=cut_sqw(sqw_file400,proj,[-0.04,0.04],[-0.04,0.04],[-0.04,0.04],[1,0.12,6],'-nopix');
acolor blue
%plot(cx(1))
plot(bose(cz(1)-2e-5,100))
acolor red
%pp(cx(2))
pp(2*bose(cz(2)-2e-4,400))
keep_figure;
title('');
set(gca,'FontSize',14);
ylabel('Intensity (arb. units)');
xlabel('Energy (meV)');
grid on
ll=legend('100 K','400 K','Location','NorthWest');
set(ll,'FontSize',14);
legend('boxoff');
text(1.1,1e-4,'\bf Q \rm = (-2, -2, 2.5) = Z','FontSize',14);

%M2 point
proj.u=[1,1,0]; proj.v=[0,0,1]; proj.type='rrr'; proj.uoffset=[-2,-2.5,2.5,0];
cm(3)=cut_sqw(sqw_file100,proj,[-0.04,0.04],[-0.04,0.04],[-0.04,0.04],[1,0.08,6],'-nopix');
cm(4)=cut_sqw(sqw_file400,proj,[-0.04,0.04],[-0.04,0.04],[-0.04,0.04],[1,0.08,6],'-nopix');
acolor blue
%plot(cx(1))
plot(bose(cm(3)-2e-5,100))
acolor red
%pp(cx(2))
pp(2*bose(cm(4)-2e-4,400))
%
keep_figure;
title('');
set(gca,'FontSize',14);
ylabel('Intensity (arb. units)');
xlabel('Energy (meV)');
grid on
ll=legend('100 K','400 K','Location','NorthWest');
set(ll,'FontSize',14);
legend('boxoff');
text(1.1,1.7e-4,'\bf Q \rm = (-2, -2.5, 2.5) = M2','FontSize',14);

%R point
proj.u=[1,1,0]; proj.v=[0,0,1]; proj.type='rrr'; proj.uoffset=[-2.5,-2.5,2.5,0];
cr(1)=cut_sqw(sqw_file100,proj,[-0.04,0.04],[-0.04,0.04],[-0.04,0.04],[1,0.08,6],'-nopix');
cr(2)=cut_sqw(sqw_file400,proj,[-0.04,0.04],[-0.04,0.04],[-0.04,0.04],[1,0.08,6],'-nopix');
acolor blue
amark 8
%plot(cx(1))
plot(bose(cr(1)-2e-5,100))
acolor red
%pp(cx(2))
pp(2*bose(cr(2)-2e-4,400))
%
keep_figure;
title('');
set(gca,'FontSize',14);
ylabel('Intensity (arb. units)');
xlabel('Energy (meV)');
grid on
ll=legend('100 K','400 K','Location','NorthWest');
set(ll,'FontSize',14);
legend('boxoff');
text(1.3,2.7e-4,'\bf Q \rm = (-2.5, -2.5, 2.5) = R','FontSize',14);