clear all; close all;

cvec =[ 0,0,0.85; 0,0.5,0; 0.9,0.5,0; 0.9,0,0; 0.2 0.2 0.2];
cvec2=[ 0,0.7,1;  0,0.8,0; 1,0.7,0;   1.0,0.4,0; 0.2 0.2 0.2];

figure2=figure('units','centimeters','position',[5,4,18,14],'color','white','DefaultTextInterpreter','LaTex','DefaultLineLineWidth',1.5);
set(gca,'LineWidth',1.5,'FontSize',16,'Layer','top','TickLabelInterpreter','LaTex'); box(gca,'on'); hold(gca,'all');

filelist=[011532]; x=[0.3616 ]; dksc = [ -0.0274  ];
filelist=[011533]; x=[x, 0.73  ]; dksc = [dksc, -0.0562  ];
filelist=[011534]; x=[x, 1.10  ]; dksc = [dksc, -0.0866   ];
filelist=[011535]; x=[x, 1.49  ]; dksc = [dksc, -0.1174   ];
filelist=[011536]; x=[x, 1.81 ]; dksc = [dksc, -0.1544  ];
filelist=[011537]; x=[x, 2.16   ]; dksc = [dksc, -0.1909   ];
filelist=[011538]; x=[x, 2.54   ]; dksc = [dksc, -0.2270   ];
filelist=[011539]; x=[x, 2.99 ]; dksc = [dksc,  -0.2608];
filelist=[011540]; x=[x, 3.4 ]; dksc = [dksc, -0.2985];
filelist=[011541]; x=[x, 3.88 ]; dksc = [dksc, -0.3342  ];
filelist=[011542]; x=[x, 4.74 ]; dksc = [dksc, -0.4158   ];
filelist=[011543]; x=[x, 5.15]; dksc = [dksc, -0.4598  ];
filelist=[011544]; x=[x, 5.56  ]; dksc = [dksc, -0.5051  ];
filelist=[011545]; x=[x, 5.98 ]; dksc = [dksc, -0.5514   ];
filelist=[011546]; x=[x, 6.41 ]; dksc = [dksc, -0.5987   ];
filelist=[011547]; x=[x, 3.51 -0.36 ]; dksc = [dksc, 0.3152    0.0257  ];
filelist=[011548]; x=[x, 5.5 -0.71 ]; dksc = [dksc, 0.5086    0.0503   ];
filelist=[011549]; x=[x, 6.71 -0.95 ]; dksc = [dksc, 0.6494    0.0830   ];
filelist=[011550]; x=[x, 7.6 -1.2 ]; dksc = [dksc, 0.7704    0.1134  ];
filelist=[012012]; x=[x,7.75 -1.25 ]; dksc = [dksc,0.7791    0.1099];
filelist=[012013]; x=[x,8.08 -1.37 -2.01]; dksc = [dksc,0.8327    0.1252    0.0627];
filelist=[012014]; x=[x,8.35 -1.53 -2.2]; dksc = [dksc,0.8827    0.1362    0.0693];
filelist=[012015]; x=[x,8.55 -1.678 -2.36]; dksc = [dksc,0.9286    0.1480    0.0783];
filelist=[012016]; x=[x,8.70 -1.826 -2.58]; dksc = [dksc,0.9715    0.1593    0.0804];
filelist=[012017]; x=[x,8.81 -1.97 -2.74]; dksc = [dksc,1.0121    0.1706    0.0882];
filelist=[012018]; x=[x,8.91 -2.11 -3.03]; dksc = [dksc,1.0527    0.1818    0.0805];
filelist=[012019]; x=[x,8.92 -2.22 -3.17]; dksc = [dksc,1.0864    0.1959    0.0891];
filelist=[012020]; x=[x,8.92 -2.37 -3.23]; dksc = [dksc, 1.1201    0.2052    0.1068];
filelist=[012021]; x=[x,7.32 -1.086 -1.576]; dksc = [dksc, 0.7195    0.0984    0.0527];
filelist=[012022]; x=[x,6.8 -0.95 -1.41]; dksc = [dksc,0.6538    0.0838    0.0418];
filelist=[012023]; x=[x,6.2 -0.79 -1.19]; dksc = [dksc,0.5844    0.0711    0.0354];
filelist=[012024]; x=[x,5.51 -0.67]; dksc = [dksc,0.5066    0.0544  ];
filelist=[012025]; x=[x,4.32]; dksc = [dksc,-0.3765  ];
filelist=[012026]; x=[x, 4.51 ]; dksc = [dksc, -0.3984 ];

%folding back into first bz along GX
bz=2.182/2;
dks=dksc;
dks = mod(dks,2*bz); 
dks(dks>bz)=abs(dks(dks>bz)-2*bz);
plot(abs(dks),abs(x),'Marker','o','MarkerSize',7,'LineStyle','None','Color',cvec(1,:),'MarkerFaceColor',cvec2(1,:));

dk=linspace(0,bz,101); cm2meV=1/8.065;
%phonon dispersion at room temperature, adpated from to Kostov et al.
A1=70*cm2meV; A3=-1.3*cm2meV; A5=+0.7*cm2meV; 
Eph=A1*sin(dk/(2*bz/pi)) + A3*sin(3*dk/(2*bz/pi)) + A5*sin(5*dk/(2*bz/pi));
p(1)=plot(dk,Eph,'Color',cvec(4,:));
A1=12.2; A3=0; A5=0; 
Eph=A1*sin(dk/(2*bz/pi)) + A3*sin(3*dk/(2*bz/pi)) + A5*sin(5*dk/(2*bz/pi));
p(2)=plot(dk,Eph,'Color',cvec(2,:));
A1=18; A3=1.5; A5=0.3; 
Eph=A1*sin(dk/(2*bz/pi)) + A3*sin(3*dk/(2*bz/pi)) + A5*sin(5*dk/(2*bz/pi));
p(3)=plot(dk,Eph,'Color',cvec(3,:));

%GX direction
filelist=[012062]; x=[10.1 -2.11 -3.07]; dksc = [1.1601    0.2077    0.1008];
filelist=[012027]; x=[3.93 -0.61 ]; dksc = [0.4044    0.0596  ];
filelist=[012028]; x=[x, 5.54 -0.91 -1.5 ]; dksc = [dksc, 0.5741    0.0874    0.0335  ]; 
filelist=[012029]; x=[x, 6.99 -1.15 -1.79 ]; dksc = [dksc, 0.7316    0.1192    0.0584  ];
filelist=[012030]; x=[x, 8.4 -1.45 -2.17 ]; dksc = [dksc, 0.8858    0.1440    0.0724   ]; 
filelist=[012031]; x=[x, 9.36 -1.71 -2.53 ]; dksc = [dksc, 1.0122    0.1710    0.0858]; %RW stronger than specular!
filelist=[012032]; x=[x, 9.85 -1.99 -2.89 ]; dksc = [dksc, 1.1105    0.1945    0.0965]; %RW stronger than specular!
filelist=[012033]; x=[x, 10.53 -2.23 -3.18 ]; dksc = [dksc, 1.2210    0.2205    0.1127]; %RW stronger than specular!
filelist=[012034]; x=[x,10.3 -2.35 -3.37 ]; dksc = [dksc,1.2413    0.2330    0.1145]; %RW stronger than specular!
filelist=[012035]; x=[x,10.53 -2.41 -3.49 ]; dksc = [dksc,1.2901    0.2518    0.1242]; %RW stronger than specular!
filelist=[012036]; x=[x,10.61 -2.53 -3.58 ]; dksc = [dksc,1.3297    0.2637    0.1372]; %RW stronger than specular!
filelist=[012037]; x=[x,10.78 -2.65 -3.71 ]; dksc = [dksc, 1.3750    0.2751    0.1448]; %RW stronger than specular!
filelist=[012038]; x=[x,10.8 -2.77 -3.85 ]; dksc = [dksc,1.4108    0.2862    0.1505]; %RW stronger than specular!
filelist=[012039]; x=[x,10.78 -2.89 -3.97]; dksc = [dksc,1.4441    0.2969    0.1583]; %RW stronger than specular!
filelist=[012040]; x=[x,10.8 -3.01 -4.21]; dksc = [dksc, 1.4800    0.3071    0.1487]; %RW stronger than specular!
filelist=[012041]; x=[x,10.9 -3.13 -4.35]; dksc = [dksc, 1.5210    0.3169    0.1522]; %RW stronger than specular!
filelist=[012042]; x=[x,11.0 -3.25 ]; dksc = [dksc, 1.5621    0.3263]; %RW stronger than specular!
filelist=[012043]; x=[x,10.8 -3.31 -4.65]; dksc = [dksc, 1.5835    0.3432    0.1539]; %RW stronger than specular!
filelist=[012044]; x=[x,0.67 ]; dksc = [dksc, -0.0611  ]; 
filelist=[012045]; x=[x,2.28 1.33 ]; dksc = [dksc,-0.0668   -0.1292  ]; 
filelist=[012046]; x=[x,3.43 2.01 ]; dksc = [dksc,-0.1167   -0.2017   ]; 
filelist=[012047]; x=[x,4.73 2.72  ]; dksc = [dksc,-0.1687   -0.2781     ]; 
filelist=[012048]; x=[x,6.21 3.53  ]; dksc = [dksc,-0.2228   -0.3548   ]; 
filelist=[012048]; x=[x,6.21 3.53  ]; dksc = [dksc,-0.2228   -0.3548   ]; 
filelist=[012049]; x=[x,7.7  4.25 ]; dksc = [dksc,-0.2879   -0.4418 ]; 
filelist=[012050]; x=[x,8.3  4.71 ]; dksc = [dksc,-0.3301   -0.4827 ]; 
filelist=[012051]; x=[x,9.14  5.16 ]; dksc = [dksc,-0.3652   -0.5256 ]; 
filelist=[012052]; x=[x,9.81  5.59 ]; dksc = [dksc,-0.4091   -0.5709 ]; 
filelist=[012053]; x=[x,7.34  6.05 ]; dksc = [dksc,-0.5672   -0.6164 ]; 
filelist=[012054]; x=[x,8.09  6.49 ]; dksc = [dksc,-0.6063   -0.6642];
filelist=[012055]; x=[x,8.91  6.92 ]; dksc = [dksc,-0.6456   -0.7137];
filelist=[012056]; x=[x,9.45  7.49 ]; dksc = [dksc,-0.6963   -0.7597];
filelist=[012057]; x=[x,9.9  7.89 ]; dksc = [dksc,-0.7512   -0.8129];
filelist=[012058]; x=[x,10.53  8.32 ]; dksc = [dksc,-0.8022   -0.8664];
filelist=[012059]; x=[x,10.8  8.75 ]; dksc = [dksc,-0.8645   -0.9210];
filelist=[012060]; x=[x,11.3  9.13 ]; dksc = [dksc, -0.9217   -0.9781];
filelist=[012061]; x=[x,11.95  9.55 ]; dksc = [dksc,-0.9765   -1.0352];
filelist=[012062]; x=[x, 10.1 -2.11 -3.07]; dksc = [dksc, 1.1601    0.2077    0.1008];
filelist=[012063]; x=[x,10.2 -2.15 -3.07]; dksc = [dksc,1.19    0.2077    0.1008];


%folding back into first bz along GM
bz=2.182/2*sqrt(2);
dks=dksc;
dks = mod(dks,2*bz); 
dks(dks>bz)=abs(dks(dks>bz)-2*bz);
plot(-abs(dks),abs(x),'Marker','o','MarkerSize',7,'LineStyle','None','Color',cvec(1,:),'MarkerFaceColor',cvec2(1,:));
hold on;

dk=linspace(0,bz,101); cm2meV=1/8.065;
A1=87*cm2meV; A3=-2.5*cm2meV; A5=-1*cm2meV; 
Eph=A1*sin(dk/(2*bz/pi)) + A3*sin(3*dk/(2*bz/pi)) + A5*sin(5*dk/(2*bz/pi));
plot(-dk,Eph,'-','Color',cvec(4,:));
A1=14; A3=0; A5=0; 
Eph=A1*sin(dk/(2*bz/pi)) + A3*sin(3*dk/(2*bz/pi)) + A5*sin(5*dk/(2*bz/pi));
plot(-dk,Eph,'Color',cvec(2,:));
A1=20; A3=1.5; A5=0.3; 
Eph=A1*sin(dk/(2*bz/pi)) + A3*sin(3*dk/(2*bz/pi)) + A5*sin(5*dk/(2*bz/pi));
plot(-dk,Eph,'Color',cvec(3,:));

xlim([-2.182/2*sqrt(2),2.182/2]); ylim([0,13])
line([0,0],[0,13],'Color',[0 0 0],'LineWidth',1.5);

xtickformat('%.1f');
text(-2.182/2*sqrt(2),13,'$\overline{\mathrm{M}}$','BackgroundColor',[1 1 1],'VerticalAlignment','bottom','HorizontalAlignment','center','FontSize',18);
text(0,13,'$\overline{\mathrm{\Gamma}}$','BackgroundColor',[1 1 1],'VerticalAlignment','bottom','HorizontalAlignment','center','FontSize',18); 
text(2.182/2,13,'$\overline{\mathrm{X}}$','BackgroundColor',[1 1 1],'VerticalAlignment','bottom','HorizontalAlignment','center','FontSize',18); 

xlabel('Wavevector $(\mathrm{\AA}^{-1})$'); ylabel('Phonon energy (meV)');
l1=legend(p,'Rayleigh mode','longitudinal resonance','bulk band edge');
set(l1,'Interpreter','LaTeX');