Int=integral_sin(x, f, k)

% Calculates integral $ Int(k)= \int f(x)*sin(k*x)*dx $  from min(x) to max(x)
% USAGE: Int=integral_sin(x, f, k)
% INPUTs:  x is the vector of x-coordinate values (not necessary on a uniform grid)
%          f is the vector of f(x) values, must have the same length as vectors x
%          k is the array of k values (this array can be of arbitraly dimentions and size)
% OUTPUTs: Int is the array of integration results, has the same size as array k
% NOTES: f(x) is approximated by 2nd-order polynomial on each double-interval of the
%        numerical integration (i.e. we use the same parabolic approximation, as in 
%        the Simpson integration method). 
%        Next, for |k*x|>=3 this polynomial is analytically (not numerically!) 
%        integrated with sin(k*x) on the double-interval. This approach greately improves
%        accuracy for large values of k, when sin(k*x) oscilates a lot. The full integral
%        is the sum of the integrals over all double-intervals (and over the last single 
%        interval if the number of intervals is odd). 
%        For |k*x|<=3 we MUST take a different approach to avoid numerical errors in the 
%        calculation of numerical differences of near function values. As a result, 
%        for |k*x|<=3 we use Taylor expansion formula for sin(k*x), which is accurate
%        and fast to evaluate (and it is what computer hardware uses anyway!).
%        For any value of k (large and small) our method has the same order of accuracy
%        of numerical integration as the (non-symmetric) Simpson method has.
% Last modified: Apr 2007.

function Int=integral_sin(x, f, k)

% for |k*x|<kx_expansion_radius we use Taylor expansion of sin(k*x),
%   $ sin(k*x) = \sum_{m=0}^M (-1)^m*(kx)^{2m+1}/(2m+1)! + O[(kx)^{2M+3}/(2M+3)!] $, where
%   M=taylor_expansion_order, the optimal choice is 1<=kx_expansion_radius<=10
taylor_expansion_order=14;
kx_expansion_radius=3.0;

% check that x and f are of the same length
N=numel(x);
if numel(f)~=N
  disp(error('x and f are NOT of the same length\n'));
end

% set arrays and constants
Int=zeros(size(k));

% return if all elements of k are zero
if ~any(k)
  return; 
end

% run Simpson integration over double-intervals
n=1;
while n<=N-2
  % set constants
  n0=n;
  n1=n+1;
  n2=n+2;
  % find a*x^2+b*x+c 2nd order polynomial approximation of f(x)
  f_prime_0=(f(n1)-f(n0))/(x(n1)-x(n0));
  f_prime_1=(f(n2)-f(n1))/(x(n2)-x(n1));
  a=(f_prime_1-f_prime_0)/(x(n2)-x(n0));
  b=f_prime_1-a*(x(n2)+x(n1));
  c=f(n1)-x(n1)*(a*x(n1)+b);
  % calculate the integral from x(n0) to x(n2)
  tmp=abs(x(n1)*k);
  index_small=find(tmp<=kx_expansion_radius);
  index_large=find(tmp>kx_expansion_radius);
  clear tmp;
  if ~isempty(index_small)
    % for speed and numerical precision, use Taylor expansion of sin(k*x) for small |k*x|
    k_small=k(index_small);
    k2=k_small.^2;
    Int_small=zeros(size(k_small));
    taylor_factor=ones(size(k_small));
    for m=0:taylor_expansion_order
      mc=2*m+2; mb=2*m+3; ma=2*m+4;
      Int_small=Int_small+taylor_factor.*...
                (c*(x(n2)^mc-x(n0)^(mc))/mc+b*(x(n2)^mb-x(n0)^mb)/mb+a*(x(n2)^ma-x(n0)^ma)/ma);
      taylor_factor=(-1/(mc*mb))*k2.*taylor_factor;
    end
    Int(index_small)=Int(index_small)+k_small.*Int_small;
  end
  if ~isempty(index_large)
    % use sin(k*x) function for large |k*x|; this is the trick -- use parabolic approximation for
    % smooth function f(x)=a*x^2+b*x+c, and calculate integral $ \int (a*x^2+b*x+c)*cos(k*x)*dx $
    % analytically, which gives high accuracy for large values of k (otherwise, we would need
    % a very small step dx to approximate function sin(k*x) that oscilates like crazy).
    k_large=k(index_large);
    inv_k=1./k_large;
    two_inv_k2=2*inv_k.^2;
    xk_2=x(n2)*k_large;
    xk_0=x(n0)*k_large;
    Int(index_large)=Int(index_large)+inv_k.*...
        ( (a*(two_inv_k2-x(n2)^2)-b*x(n2)-c).*cos(xk_2)-...
          (a*(two_inv_k2-x(n0)^2)-b*x(n0)-c).*cos(xk_0)+...
          inv_k.*(2*a*x(n2)+b).*sin(xk_2)-...
          inv_k.*(2*a*x(n0)+b).*sin(xk_0) );
  end
  % advance n
  n=n+2;
end
% make integration over the last interval, if necessary
if n==N-1
    % set constants
    n0=n-1;
    n1=n;
    n2=n+1;
    % find a*x^2+b*x+c 2nd order polynomial approximation of f(x)
    f_prime_0=(f(n1)-f(n0))/(x(n1)-x(n0));
    f_prime_1=(f(n2)-f(n1))/(x(n2)-x(n1));
    a=(f_prime_1-f_prime_0)/(x(n2)-x(n0));
    b=f_prime_1-a*(x(n2)+x(n1));
    c=f(n1)-x(n1)*(a*x(n1)+b);
    % calculate the integral from x(n1) to x(n2)
    tmp=abs(x(n1)*k);
    index_small=find(tmp<=kx_expansion_radius);
    index_large=find(tmp>kx_expansion_radius);
    clear tmp;
    if ~isempty(index_small)
      % for speed and numerical precision, use Taylor expansion of sin(k*x) for small |k*x|
      k_small=k(index_small);
      k2=k_small.^2;
      Int_small=zeros(size(k_small));
      taylor_factor=ones(size(k_small));
      for m=0:taylor_expansion_order
        mc=2*m+2; mb=2*m+3; ma=2*m+4;
        Int_small=Int_small+taylor_factor.*...
                  (c*(x(n2)^mc-x(n1)^(mc))/mc+b*(x(n2)^mb-x(n1)^mb)/mb+a*(x(n2)^ma-x(n1)^ma)/ma);
        taylor_factor=(-1/(mc*mb))*k2.*taylor_factor;
      end
      Int(index_small)=Int(index_small)+k_small.*Int_small;
    end
    if ~isempty(index_large)
      % use sin(k*x) function for large |k*x|
      k_large=k(index_large);
      inv_k=1./k_large;
      two_inv_k2=2*inv_k.^2;
      xk_2=x(n2)*k_large;
      xk_1=x(n1)*k_large;
      Int(index_large)=Int(index_large)+inv_k.*...
          ( (a*(two_inv_k2-x(n2)^2)-b*x(n2)-c).*cos(xk_2)-...
            (a*(two_inv_k2-x(n1)^2)-b*x(n1)-c).*cos(xk_1)+...
            inv_k.*(2*a*x(n2)+b).*sin(xk_2)-...
            inv_k.*(2*a*x(n1)+b).*sin(xk_1) );
    end
    % advance n
    n=n+1;
end