function signal=sqw_egelstaff(varargin)
% model = sqw_egelstaff(p, q ,w, {signal}) : Recoil dispersion(Q) for a single harmonic oscillator
%
%   iFunc/sqw_egelstaff: a 2D S(q,w) with a recoil dispersion, isotropic harmonic oscillator 
%
%   This is a 2D recoil model for a single harmonic oscillator particle with given mass.
%   The mass defines the dispersion energy, while the harmonic oscillator energy
%   defines the Debye-Waller function, i.e. intensity vs Q.
%   This model is also known as 'impulse approximation' or 'short time' Gaussian. 
%
%   The dispersion has the form:
%      S(q,w) = 1/sqrt(2*pi*delta2).*exp( -(w - Er).^2/2./delta2 )
%   where
%      Er     = h2 q2/2/m
%      delta2 = Er.*coth(w0/2/kT)
%
%   You can build a recoil model for a given mass and oscillator energy with syntax:
%      sqw = sqw_egelstaff(mass)
%   or
%      sqw = sqw_egelstaff([ mass energy ])
%
%   You can of course tune other parameters once the model object has been created.
%   To model more than one oscillator/particle recoil, just add these models together.
%
%   Evaluate the model with syntax:
%     sqw(p, q, w)
%
% input:  p: sqw_egelstaff model parameters (double)
%             p(1)= Mass       Mass of the scattering unit [g/mol]
%             p(2)= Ex_Energy  Harmonic Excitation Energy [meV]
%             p(3)= Amplitude 
%             p(4)= Temperature  Temperature [K]  
%         q:  axis along wavevector/momentum in Angs-1 (row,double)
%         w:  axis along energy in meV (column,double)
% output: signal: model value
%
%
% Example:
%   s=sqw_egelstaff(2); % Deuterium mass
%   plot(iData(s, [], 0:.1:20, (0:100)))  % q=0:20, energy=0:100
%
% Reference: Schober, JNR 17 (2014) 109–357 - DOI 10.3233/JNR-140016
%
% Version: $Date$
% See also iData, iFunc
%   <a href="matlab:doc(iFunc,'Models')">iFunc:Models</a>
% (c) E.Farhi, ILL. License: EUPL.

signal.Name           = [ 'sqw_egelstaff dispersion(Q) for given mass and excitation energy [' mfilename ']' ];
signal.Description    = 'A 2D S(q,w) with a recoil dispersion.';

signal.Parameters     = {  ...
'Mass         Mass of the scattering unit [g/mol]' ...
'Diffusion    Diffusion coefficient D [mm2/s]' ...
'Amplitude' ...
'Temperature  Temperature [K]' ...
};
  
signal.Dimension      = 2;         % dimensionality of input space (axes) and result
signal.Guess = [ 1 10 1 10 ];

signal.Expression     = { ...
 'm = p(1); D= p(2); T = p(4);' ...
 'q = x; w = y;' ...
 'if isvector(q) && isvector(w) && numel(q) ~= numel(w), [q,w] = meshgrid(q,w); end' ...
 'kb      = 1.38064852E-023;  % Boltzmann [J/K]' ...
 'mn      = 1.674927471E-027; % neutron mass [kg]' ...
 'e       = 1.60217662E-019;  % [C]' ...
 'HBAR    = 1.05457168e-34;   % Plank/2PI [m2 kg / s]' ...
 'meVtoK  = e/1000/kb;        % = 11.6045 = [meV] to [K]' ...
 'meV2Hz  = e/HBAR/1000/2/pi; % 1 meV = 241.80e9 GHz ' ...
 'kT      = T/meVtoK;                   % T [K] -> [meV] = 11.6045' ...
 'DQ2     = D*q.^2*1e14 / meV2Hz/2/pi;  % converted to meV' ...
 'tau     = m*D/kT; % inverse of an energy [g/mol*cm2/s]/[meV-1]' ...
 't0      = 1/2/kT;' ...
 'sq      = sqrt( w.^2+DQ2.^2 );
 'lambda  = sqrt( 1 + (1/2/m/D)^2 );' ...
 'signal  = DQ2*tau*lambda/pi .* exp(w/2/kT + tau*DQ2);' ...
 'signal  = signal .* besselk(1, tau*lambda*sq) ./ sq;' ...
 };
 
 [mm^2/s] *[1e20/m^2] = [1e-6*m^2/s] *1e20/m^2 = 1e14*Hz=1e14/meV2Hz/2/pi
 
 Da         = 1.660538921e-27;
 m*D/kb/T = [kg*Da]*[mm^2/s]/kb/[J]= [kg] * 1e-6*[m^2/s/kb/J]*Da
          = [kg m2/s]*1e-6*Da/kb / J
 
 
 meV2Hz = HBAR/2pi in meV
 
 kT = kb * J/K * K = kb * J
 1 meV = e/1000 J
 
 
 dq2 = [cm2/s*Angs-2] 
 2*pi*(1/s)/(1e10/m)/(1e10/m)/1e-14
 2*pi*(1/s)/(1e20/m^2)/1e-14
 2*pi/s*m^2/1e6
 
 % D   = 2*PI*w_c*meV2Hz/q_min/q_min/1e14*RMS2FWHM/2; /* meV*Angs^2 -> mm^2/s */

signal= iFunc(signal);
signal= iFunc_Sqw2D(signal); % overload Sqw2D flavour

if nargin == 1 && isnumeric(varargin)
  p = varargin{1};
  if numel(p) >= 1, signal.ParameterValues(1) = p(1); end
  if numel(p) >= 2, signal.ParameterValues(2) = p(2); end
end