function int=empirical_intensity(h,k,l,w,mode,disp_pars,int_pars)
% function int=empirical_intensity(h,k,l,w,mode,disp_pars,int_pars)
%
% The total dynamical structure factor for specified one-magnon modes in a
% spin-1/2 triangular antiferromagnet in the lab frame, calculated using 
% the analytic LSWT result with empirical dispersion renormalizations that
% parametrize the experimental neutron scattering data collected for 
% Ba3CoSb2O9 (see Appendix).
%
% In this code, the following notation convention is used: 
% (x,y,z) - rotating frame coordinates; (x0,y0,z0) - lab frame coordinates
%
% INPUT
% h, k, l   - arrays of wave vector components (r.l.u.)
% w         - array of energies (meV) 
%             [if size of w doesn't match h,k,l, S(k) calculated instead]             
% mode      - array of mode labels [1 -> w(k), 2 -> w(k-Q), 3 -> w(k+Q)]
% disp_pars - array of dispersion parameters, as in empirical_dispersion.m
% int_pars  - array of intensity parameters: 
%             [amp Zx Zy bkg fwhm T mff pol]
%               1  2  3   4    5  6  7   8
%             amp    - overall intensity scale factor
%             Zx, Zy - Sxx anf Syy relative scale factors
%             bkg    - flat background
%             fwhm   - full-width half-maximum (meV) of Gaussian broadening
%             T      - temperature (K) used in Bose factor
%             mff    - magnetic form factor flag (included if true)
%             pol    - neutron polarization flag (include if true)
%
% N.B. If int_pars not specified, assumed int_pars = [1 1 1 0 0.01 0 0 0].
%
% OUTPUT
% int    - total dynamical structure factor in lab frame for chosen modes
%
% See empirical_model_example.m for an example script using this function.

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

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

% Specify default intensity parameters
if ~exist('int_pars','var') || isempty(int_pars)
    int_pars=[1 1 1 0 0.01 0 0 0];
end

% Extract intensity parameters
amp=int_pars(1);
Zx=int_pars(2);
Zy=int_pars(3);
bkg=int_pars(4);
fwhm=int_pars(5);
T=int_pars(6);
mff_flag=int_pars(7);
pol_flag=int_pars(8);

% Define Fourier transform of exchange interactions [see Eq.(A13)]
Jq=@(qh,qk,ql) disp_pars(1)*(cos(2*pi*qh)+cos(2*pi*qk)+cos(2*pi*(qh+qk)))...
    +disp_pars(2)*cos(pi*ql);

% Define canting angle in applied field [see Eq.(A15)]
theta=asin(4*disp_pars(4)/17.276/(disp_pars(3)*Jq(0,0,0)-Jq(Q(1),Q(2),Q(3))));

% Calculate |Q|=|(h,k,l)| in 1/Angs.
modQ=sqrt((h*ar(1)+k*br(1)+l*cr(1)).^2+(h*ar(2)+k*br(2)+l*cr(2)).^2 ...
    +(h*ar(3)+k*br(3)+l*cr(3)).^2);  

% Apply neutron polarization factor
if pol_flag 
    pol=(h*ar(3)+k*br(3)+l*cr(3))./modQ;  % assume spins ordered in ab plane
    pz0=1-pol.^2;
    px0=1+pol.^2;  % including y0 component
else  % no polarization factor
    pz0=1;
    px0=2;  % including y0 component
end

% Sum intensity contributions
int=zeros(size(h));
if ismember(1,mode)  % add w(k)
    int=int+Zy*cos(theta)^2*Syy(h,k,l,disp_pars)...
        .*lineshape(h,k,l,w,disp_pars,fwhm).*pz0;  % Eq.(A24)
end
if ismember(2,mode)  % add w(k-Q)
    int=int+(1/4)*(Zx*Sxx(h-Q(1),k-Q(2),l-Q(3),disp_pars) + ...
        Zy*sin(theta)^2*Syy(h-Q(1),k-Q(2),l-Q(3),disp_pars) - 2i*sin(theta)*Sxy) ...
        .*lineshape(h-Q(1),k-Q(2),l-Q(3),w,disp_pars,fwhm).*px0;  % Eq.(A26)
end
if ismember(3,mode)  % add w(k+Q)
    int=int+(1/4)*(Zx*Sxx(h+Q(1),k+Q(2),l+Q(3),disp_pars) + ...
        Zy*sin(theta)^2*Syy(h+Q(1),k+Q(2),l+Q(3),disp_pars) + 2i*sin(theta)*Sxy) ...
        .*lineshape(h+Q(1),k+Q(2),l+Q(3),w,disp_pars,fwhm).*px0;  % Eq.(A27)
end

% Apply finite-temperature Bose factor
if T>0
    kB=8.617e-2;   % Boltzmann constant (meV/K)
    w(w==0)=0.01;  % remove divergent terms
    int=int./(1-exp(-w/(kB*T)));
end

% Apply spherical magnetic form factor for Co^{2+} ions
% (using values in https://www.ill.eu/sites/ccsl/ffacts/ffactnode5.html)
if mff_flag
    z=(modQ/4/pi).^2;  % define exponent
    fQ=0.4332*exp(-14.3553*z)+0.5857*exp(-4.6077*z)-0.0382*exp(-0.1338*z)+0.0179;
    int=int.*(fQ.^2);
end

% Apply scale factor and add flat background
int=amp*int+bkg;
int=real(int);

end

function int=Sxx(h,k,l,disp_pars)
% Sxx(k) structure factor, see Eq.(A17)
[~,u,v]=empirical_dispersion(h,k,l,1,disp_pars);
int=(1/4)*(u+v).^2;
end

function int=Syy(h,k,l,disp_pars)
% Syy(k) structure factor, see Eq.(A18)
[~,u,v]=empirical_dispersion(h,k,l,1,disp_pars);
int=(1/4)*(u-v).^2;
end

function int=Sxy
% Sxy(k) structure factor, see Eq.(A19)
int=1i/4;
end

function y=lineshape(h,k,l,w,disp_pars,fwhm)
% Intensity lineshape in energy for modes
wq=empirical_dispersion(h,k,l,1,disp_pars);
if numel(h)==numel(w)  % use Gaussian broadening
    y=gssn(w,wq,fwhm);
else  % integrate over energy to give S(k)
    y=1;
end
end

function y=gssn(x,x0,fwhm)
% Normalized Gaussian function
sig=fwhm/(2*sqrt(2*log(2)));
y=exp(-(x-x0).^2/2/sig^2)/(sqrt(2*pi)*sig);
end
