%% Simulated intensity and dispersion for best-fit empirical model
% Similar wave vector path to Fig. 6(b)
G = [0 1 0];      % Gamma point
M = [0 1/2 0];    % M point
K = [1/3 1/3 0];  % K point
L = [0 0 1];
QQ = [K; M; G; K; K+L; M+L; G+L];

% Define reciprocal lattice vectors (r.l.u.)
ar=[1.0768 0.62169 0];
br=[0 1.2434 0];
cr=[0 0 0.43487];

% Load best-fit dispersion parameters (zero-field, +ve chirality)
load('best_fit_pars.txt', '-ascii')  % see Table I

% Define intensity parameters
amp = 100;       % overall intensity scale
Zx = 1; Zy = 1;  % unchanged relative Sxx/Syy intensity
bkg = 0;         % no background
fwhm = 0.1;      % Gaussian broadening of FWHM 0.1 meV
T = 1;           % Bose factor for T = 1 K
mff = 1;         % include magnetic form factor for Co^{2+} ions
pol = 1;         % include neutron polarization factor
int_pars = [amp Zx Zy bkg fwhm T mff pol];

% Define x points (modQ) of plot
n_plot = 50;            % no. of plotted points in each panel
n_ticks = size(QQ, 1);  % no. of symmetry points (x ticks)
prev_Q = QQ(1,:);
[x, qq(:,1), qq(:,2), qq(:,3)] = deal(NaN((n_ticks-1)*n_plot+1, 1));
ticks = zeros(n_ticks, 1);
for i = 2:n_ticks
    qq_i = [linspace(prev_Q(1), QQ(i,1), n_plot+1)' ...
            linspace(prev_Q(2), QQ(i,2), n_plot+1)' ...
            linspace(prev_Q(3), QQ(i,3), n_plot+1)'];
    qq(1+(i-2)*n_plot:(i-1)*n_plot,:) = qq_i(1:end-1,:);
    qq_diff = qq_i - repmat(prev_Q, n_plot+1, 1);
    modq = sqrt(sum((qq_diff(:,1)*ar + qq_diff(:,2)*br + qq_diff(:,3)*cr).^2, 2));
    x(1+(i-2)*n_plot:(i-1)*n_plot) = ticks(i-1) + modq(1:end-1);
    ticks(i) = ticks(i-1) + modq(end);
    prev_Q = QQ(i,:);
end
qq(end,:) = qq_i(end,:); x(end) = ticks(end);

% Define energy and wave vector grids for calculation
en = 0:0.01:2;                  
[H,E]=meshgrid(qq(:,1),en);     
K=repmat(qq(:,2)',numel(en),1);
L=repmat(qq(:,3)',numel(en),1);
X=repmat(x',numel(en),1);

% Calculate intensity grid
int=empirical_intensity(H,K,L,E,[1 2 3],best_fit_pars,int_pars);

% Plot calculated intensities
figure
pcolor(X,E,int)

box on
shading flat

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

xticks(ticks)
xticklabels({'K_1' 'M_2' '\Gamma' 'K_1' 'K''_1' 'M''_2' '\Gamma'''})
set(gca,'TickDir','out')

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

% Calculate dispersions
w1=empirical_dispersion(qq(:,1),qq(:,2),qq(:,3),1,best_fit_pars);  % w(k)
w2=empirical_dispersion(qq(:,1),qq(:,2),qq(:,3),2,best_fit_pars);  % w(k-Q)
w3=empirical_dispersion(qq(:,1),qq(:,2),qq(:,3),3,best_fit_pars);  % w(k+Q)

% Overplot dispersions
hold on
h1=plot(x,w1,'m-','LineWidth',2);
h2=plot(x,w2,'m--','LineWidth',2);
h3=plot(x,w3,'m:','LineWidth',2);
hold off

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

savefig('empirical_model_example')
