function [signal, this] = sqw_phonons_template_dho(HKL, t, FREQ, POLAR, is_event, resize_me, Amplitude, Gamma, Bkg, T)
% sqw_phonon_template_dho: code to build signal from DHO's at each mode FREQ, for all HKL locations
%   when POLAR is not empty, the phonon intensity (Q.e)^2 form factor is computed
%   when resize_me is specified [size], the signal is reshaped to it
% Reference: B. Fak, B. Dorner / Physica B 234-236 (1997) 1107-1108
%            H. Schober, J Neut. Res. 17 (2014) 109

% size of 'w' is [ numel(x) numel(t) ]. the energy for which we evaluate the model

this.UserData.maxFreq= max(FREQ); % for each mode
nt = numel(t);
% store when not too large
if numel(HKL) <= 1e5, 
  this.UserData.FREQ = FREQ; 
  this.UserData.HKL  = HKL;
  this.UserData.POLAR = POLAR; 
else
  this.UserData.FREQ = []; 
  this.UserData.HKL  = [];
  this.UserData.POLAR = [];
end

% test for unstable modes
wrong_w = numel(find(FREQ(:) < 0 | ~isreal(FREQ(:)) | ~isfinite(FREQ(:))));
if wrong_w, 
  disp([ 'WARNING: found ' num2str(wrong_w) ...
         ' negative/imaginary phonon frequencies (' ...
         num2str(wrong_w*100/numel(FREQ)) '% of total) in ' this.Name ]); 
end
clear wrong_w

% we compute the Q vector in [Angs-1]. Search for the B=rlu2cartesian matrix
UD = this.UserData; B=[];
if isfield(UD, 'reciprocal_cell')
  B = UD.reciprocal_cell;
elseif isfield(UD, 'properties') && isfield(UD.properties, 'reciprocal_cell')
  B = UD.properties.reciprocal_cell;
else
  B = eye(3); % assume cubic, a=b=c=2*pi, 90 deg, a*=2pi/a=1
end
q_cart = B*HKL';
qx=q_cart(1,:); qy=q_cart(2,:); qz=q_cart(3,:);
Q=[ qx(:) qy(:) qz(:) ];  % in Angs-1
clear q_cart qx qy qz

% check if we have everything needed for the intensity estimate: b_coh, positions
b_coh = []; positions = []; masses=1;
if isfield(UD, 'positions')
  positions = UD.positions;
elseif isfield(UD, 'properties') && isfield(UD.properties, 'positions')
  positions = UD.properties.positions;
end

if isfield(UD, 'masses')
  masses = UD.positions;
elseif isfield(UD, 'masses') && isfield(UD.properties, 'masses')
  masses = UD.properties.masses;
end
if isscalar(masses) && size(positions,1) > 1
  for index=1:size(positions,1)
    if index > numel(masses), masses(index) = 1; end
    if isnan(masses(index)) || masses(index) <= 0, masses(index) = 1; end
  end
end

if isfield(UD, 'b_coh')
  b_coh = UD.b_coh;
elseif isfield(UD, 'properties') && isfield(UD.properties, 'b_coh')
  b_coh = UD.properties.b_coh;
end

if isempty(b_coh) || any(b_coh == 0 | isnan(b_coh))
  if isfield(UD, 'sigma_coh')
    sigma_coh = UD.sigma_coh;
  elseif isfield(UD, 'properties') && isfield(UD.properties, 'sigma_coh')
    sigma_coh = UD.properties.sigma_coh;
  else sigma_coh = [];
  end
  if ~isempty(sigma_coh)
    b_coh = sqrt(abs(sigma_coh)*100/4/pi); % in [fm]
    b_coh = b_coh .* sign(sigma_coh);
  end
end

if ~isempty(POLAR)
  if (isempty(b_coh) || any(isnan(b_coh) | b_coh == 0)) && ~isempty(positions)
    disp([ 'WARNING: Unspecified/invalid coherent neutron scattering length specification for the material ' ...
      UD.properties.chemical_formula '. Using b_coh=1 [fm] for atoms (sigma_coh=0.126 barns). ' ...
      'Specify model.UserData.properties.b_coh as a vector with ' num2str(size(positions,1)) ' value(s) in ' this.Name ]);
    for index=1:size(positions,1)
      if index > numel(b_coh), b_coh(index) = 1; end
      if isnan(b_coh(index)) || b_coh(index) == 0, b_coh(index) = 1; end
    end
  end

  if isscalar(b_coh) && size(positions,1) > 1
    disp([ 'WARNING: The material ' UD.properties.chemical_formula ...
      ' has ' num2str(size(positions,1)) ...
      ' atoms in the cell, but only one coherent neutron scattering length is defined. ' ...
      'Using the same value for all (may be wrong). ' ...
      'Specify model.UserData.properties.b_coh as a vector in ' this.Name ]);
    b_coh = b_coh * ones(1, size(positions,1));
  end

  if ~isempty(b_coh) && ~isempty(positions) && numel(b_coh) ~= size(positions,1)
    disp([ 'WARNING: Inconsistent coherent neutron scattering length specification: has ' ...
      num2str(numel(b_coh)) ' but the material ' UD.properties.chemical_formula ' has ' ...
      num2str(size(positions,1)) ' atoms in the cell. Will not compute phonon intensities. ' ...
      'Specify model.UserData.properties.b_coh as a vector in ' this.Name ]);
      b_coh = [];
  end
end

if ~isempty(b_coh)
  this.UserData.properties.b_coh     = b_coh;
  this.UserData.properties.sigma_coh = 4*pi.*b_coh.*b_coh/100;
end

if is_event, w = t(:);
else         w = ones(size(FREQ,1),1) * t(:)'; end
signal               = zeros(size(w));

% the Bose factor is negative for w<0, positive for w>0
% (n+1) converges to 0 for w -> -Inf, and to 1 for w-> +Inf. It diverges at w=0
if ~isempty(T) && T > 0, 
  n         = 1./(exp(w/T)-1);
  n(w == 0) = 0;  % avoid divergence
else n=0; end

if Gamma<=0, Gamma=1e-4; end

% compute Gamma point modes (IR/Raman) 
index = find(sum(abs(HKL),2) == 0);
if ~isempty(index) && (~isfield(UD,'properties') ...
  || ~isfield(UD.properties,'vibrational_energies') ...
  || isempty(UD.properties.vibrational_energies))
  disp([ strtok(this.Name) ': Gamma point energies (IR/Raman):' ]);
  this.UserData.properties.vibrational_energies = FREQ(index(1),:)';
  f = this.UserData.properties.vibrational_energies(:);
  disp(' [meV]      [THz]     [cm-1]')
  disp(num2str([ f f*.2418 f*8.0657 ],'%10.3f'));
end

% convert negative frequencies into imaginary
index=find(real(FREQ) < 0);
if ~isempty(index), 
  FREQ(index) = imag(FREQ(index)) - i*real(FREQ(index));
end

for index=1:size(FREQ,2)  % loop on modes
  % % size of 'w0' is [ numel(x) numel(t) ]. the energy of the modes (columns) at given x=q (rows)
  % we assume Gamma(w) = Gamma w/w0
  % W0 is the renormalized phonon frequency W0^2 = w0^2+Gamma^2
  Gamma2 = Gamma^2 + imag(FREQ(:,index)).^2;  % imaginary frequency goes in the damping
  W02    = Gamma2  + real(FREQ(:,index)).^2;  % shifts unstable mode energies (soft modes)
  
  % phonon form factor (intensity) |Q.e|^2 for neutron scattering
  % POLAR is a set of [xyz] vectors for each atom in the cell POLAR(HKL,mode,atom,xyz)
  % polarisation vectors already contain the 1/mass factor: ASE/phonons.py:band_structure
  % ZQ=|F(Q)|^2=|sum_atom[ bcoh(atom). exp(-WQ) .* (Q.*POLAR(HKL,index,atom,:)) .* exp(-i*Q.*pos(atom)) ]|^2
  ZQ = 1;
  try
    if numel(POLAR) > 1 && ~isempty(positions)
      ZQ = 0;
      for atom=1:size(positions,1)
        % one-phonon structure/form factor in a Bravais lattice
        if ~isempty(T) && T > 0
          nw0       = 1./(exp(FREQ(:,index)/T)-1);
          nw0(FREQ(:,index) == 0) = 0;
        else nw0=0; end
        DW = abs(sum(Q.*squeeze(POLAR(:,index,atom,:)),2)).^2./FREQ(:,index).*(2*nw0+1)/2; % DW function, Schober (9.103)
        clear nw0
        DW = exp(-DW);  % Debye-Waller factor
        % one phonon form factor: H. Schober, (9.203)
        if ~isempty(b_coh) && all(b_coh)
          ZQ = ZQ + b_coh(atom) / sqrt(masses(atom)) .* DW .* sum(Q.*squeeze(POLAR(:,index,atom,:)),2) .* exp(-i*Q*positions(atom,:)');
        else  % when no b_coh, we only show DW*exp(-iQr)
          ZQ = ZQ + DW .* exp(-i*Q*positions(atom,:)');
        end
      end
      clear DW
      ZQ = abs(ZQ).^2;
    end
  catch
    disp('Failed intensity estimate (mode polarisation)');
    ZQ = 1;
  end
  if ~is_event
    if ~isscalar(ZQ), ZQ = ZQ*ones(1,nt); end
    W02    = W02   *ones(1,nt);
    Gamma2 = Gamma2*ones(1,nt);
  end
  % sum-up all contributions to signal: Fak Dorner 1997 Eq (2)
  signal = signal+ (n+1).*ZQ*4.*w.*sqrt(Gamma2)/pi ./ ((w.^2-W02).^2 + 4*w.^2.*Gamma2);
end % for mode index
clear POLAR W02 Gamma2 ZQ n w

% compute the phonon DOS: total and partials. 
% This requires that the Q range is large enough.
Qrange1 = (max(HKL(:,1))-min(HKL(:,1)) >= 1 & ...
           max(HKL(:,2))-min(HKL(:,2)) >= 1 & ...
           max(HKL(:,3))-min(HKL(:,3)) >= 1);

if Qrange1 && (~isfield(this.UserData,'DOS') || isempty(this.UserData.DOS)) && ...
    size(HKL,1) >=500
  disp([ strtok(this.Name) ': Computing the phonon DOS on ' num2str(size(HKL,1)) ' HKL points.' ]);
  f = FREQ;
  index= find(imag(f) == 0);
  [dos_e,omega_e]=hist(f(index),100);
  dos_factor = size(FREQ,2) / trapz(omega_e(:), dos_e(:));
  dos_e = dos_e * dos_factor ; % 3n modes per unit cell
  DOS=iData(omega_e,dos_e);
  DOS.Title = [ 'Total DOS ' this.Name ];
  DOS.Label = '';
  xlabel(DOS,'Energy [meV]'); 
  ylabel(DOS,[ 'Total DOS/unit cell ' strtok(this.Name) ]);
  DOS.Error=0; this.UserData.DOS=DOS;
  % partial phonon DOS (per mode)
  pDOS = [];
  for mode=1:size(FREQ,2)
    f = FREQ(:,mode);
    index = find(imag(f) == 0);
    dos_e = hist(f(index),omega_e);
    dos_e = dos_e * dos_factor ; % normalize to the total DOS
    DOS=iData(omega_e,dos_e);
    DOS.Title = [ 'Mode [' num2str(mode) '] DOS ' this.Name ]; 
    DOS.Label = '';
    xlabel(DOS,'Energy [meV]'); 
    ylabel(DOS,[ 'Partial DOS[' num2str(mode) ']/unit cell ' strtok(this.Name) ]);
    DOS.Error = 0;
    pDOS = [ pDOS DOS ];
  end
  clear f index dos_e omega_e dos_factor DOS
  this.UserData.DOS_partials=pDOS;
  clear pDOS
end

clear HKL FREQ Qrange1

% Amplitude
if Amplitude
  signal = signal*Amplitude + Bkg;
end
if ~isempty(resize_me) && prod(resize_me) == numel(signal)
  signal = reshape(signal, resize_me); % initial 4D cube dimension = [ size(x) numel(t) ]
end