function [w,u,v]=empirical_dispersion(h,k,l,mode,pars)
% function [w,u,v]=empirical_dispersion(h,k,l,mode,pars)
%
% The energy of a magnon mode in a spin-1/2 triangular antiferromagnet at a
% given wave vector, calculated using the analytic LSWT dispersions with 
% empirical renormalizations that parametrize the experimental neutron 
% scattering data collected for Ba3CoSb2O9 (see Appendix).
%
% INPUT
% h, k, l - arrays of wave vector components (r.l.u.)
% mode    - mode label [1 -> w(k), 2 -> w(k-Q), 3 -> w(k+Q)]
% pars    - array of empirical model parameters:
%
%           [J1 Jz Delta B +ch/-ch alphaM betaM gammaM deltaM kappaM ...
%            1  2    3   4    5      6      7      8      9     10 ...
%           ... alphaK/2 betaK/2 gammaK/2 deltaK/2 epsilon]
%           ...    11      12       13       14       15
%            
%           +ch/-ch - >0 for +ve chirality domain, <0 for -ve chirality
%           Other parameters have the same names as in Appendix in Table I.
%
% N.B. If only parameters 1--5 specified, no renormalization is applied.
%
% OUTPUT
% w    - energy (meV)
% u, v - the u(k) and v(k) parameters in the correlation functions
%
% See empirical_model_example.m for an example script using this function.

% Define reciprocal lattice vectors (1/Angs.)
ar=[1.0768 0.62169 0];
br=[0 1.2434 0];

% Define ordering wave vector (r.l.u.)
Q=sign(pars(5))*[1/3 1/3 1];

switch mode
    % w(k) mode by default
    case {2}, h = h - Q(1); k = k - Q(2); l = l - Q(3);  % w(k-Q) mode
    case {3}, h = h + Q(1); k = k + Q(2); l = l + Q(3);  % w(k+Q) mode
end

% LSWT calculation of dispersion
[w,u,v]=lswt_dispersion(pars(1:4),h,k,l,Q);

% Apply renormalization if relevant parameters specified
if numel(pars)>5    
    % Fold (h,k) onto 1st Brillouin zone
    [hh,kk]=fold_trbz(h,k);
    
    % Add softening at M points (see Appendix B)
    pt=[0 1/2 0; 1/2 0 -pi/3; 1/2 1/2 pi/3; -1/2 1/2 pi/3];
    kappa=pars(10);
    for i=1:size(pt,1)
        alpha=pars(6);
        beta=pars(7);
        gamma=pars(8);
        delta=pars(9)*w;
        h1=hh-pt(i,1);
        k1=kk-pt(i,2);
        Qx=h1*norm(ar)*sqrt(3)/2;
        Qy=h1*norm(ar)*1/2+k1*norm(br);
        t=Qx*cos(pt(i,3))+Qy*sin(pt(i,3));
        Qy=-Qx*sin(pt(i,3))+Qy*cos(pt(i,3));
        Qx=t;
        w2=alpha+beta*cos(pi*l)+gamma*(Qx.^2+kappa*Qy.^2); 
        w=(w+w2)/2-sqrt(((w-w2)/2).^2+delta.^2); 
    end
    
    % Add softening at K/2 points (see Appendix B)
    epsilon=pars(15);
    pt=[epsilon*[1 1]; epsilon*[-1 2]; [0 1]+epsilon*[1 -2]; ...
        [1 0]+epsilon*[-2 1]; epsilon*[2 -1]; epsilon*[-2 1]];
    for i=1:size(pt,1)
        alpha=pars(11);
        beta=pars(12);
        gamma=pars(13);
        delta=pars(14)*w;
        h1=hh-pt(i,1);
        k1=kk-pt(i,2);
        w2=alpha+beta*cos(pi*l)+gamma*(h1.^2*dot(ar,ar)...
            +k1.^2*dot(br,br)+2*h1.*k1*dot(ar,br));
        w=(w+w2)/2-sqrt(((w-w2)/2).^2+delta.^2);
    end
end
end

function [w,u,v]=lswt_dispersion(pars,h,k,l,Q)
% Analytic form of LSWT dispersion, as presented in Appendix A.
Jq=@(qh,qk,ql) pars(1)*(cos(2*pi*qh)+cos(2*pi*qk)+cos(2*pi*(qh+qk)))...
    +pars(2)*cos(pi*ql); 

theta=asin(4*pars(4)/17.276/(pars(3)*Jq(0,0,0)-Jq(Q(1),Q(2),Q(3))));

b=(Jq(h-Q(1),k-Q(2),l-Q(3))+Jq(h+Q(1),k+Q(2),l+Q(3)))/2-Jq(Q(1),Q(2),Q(3));
a=(pars(3)*Jq(h,k,l)-Jq(Q(1),Q(2),Q(3)))*cos(theta)^2+b*sin(theta)^2;

A=(a+b)/2;
B=(a-b)/2;
C=(Jq(h+Q(1),k+Q(2),l+Q(3))-Jq(h-Q(1),k-Q(2),l-Q(3)))*sin(theta)/2;  % Eq.(A1)

w=sqrt(A.^2-B.^2)+C;  % Eq.(A2)
w=real(w);

thq=atanh(B./A)/2;
u=cosh(thq);  
v=sinh(thq);
end

function [h,k]=fold_trbz(hh,kk)
% Reduce (h,k) to large BZ [0,1), [0,1) in r.l.u.
h=hh-floor(hh);
k=kk-floor(kk);
M=[1 1; -1 1];
N=transpose(inv(M));
hp=N(1,1)*h+N(1,2)*k;
kp=N(2,1)*h+N(2,2)*k;
kp=abs(kp);
hp=1/2-abs(hp-1/2);
N=transpose(M);
h=N(1,1)*hp+N(1,2)*kp;
k=N(2,1)*hp+N(2,2)*kp;

M=[1 -1/2; 0 1];
N=transpose(inv(M));
hp=N(1,1)*h+N(1,2)*k;
kp=N(2,1)*h+N(2,2)*k;
kp=1/2-abs(kp-1/2);
N=transpose(M);
h=N(1,1)*hp+N(1,2)*kp;
k=N(2,1)*hp+N(2,2)*kp;

M=[1 1; -1 1];
N=transpose(inv(M));
hp=N(1,1)*h+N(1,2)*k;
kp=N(2,1)*h+N(2,2)*k;
kp=abs(kp);
N=transpose(M);
h=N(1,1)*hp+N(1,2)*kp;
k=N(2,1)*hp+N(2,2)*kp;
end
