Int=integral_simpson(x, f)

% Calculates integral $ Int = \int f(x)*dx $  from min(x) to max(x) on a non-uniform grid
% USAGE: Int=integral_simpson(x, f)
% 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
% OUTPUTs: Int is the result of the integration
% NOTES: f(x) is approximated by 2nd-order polynomial on each double-interval of the
%        numerical integration (this parabolic approximation is the Simpson integration
%        method). Next, this polynomial is analytically integrated on the double-interval. 
%        The full integral is the sum of the integrals over all double-intervals (and over
%        the last single interval if the later exists).
% Last modified: Apr 2007.

function Int=integral_simpson(x, f)

% 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 constants
Int=0;
inv_3=1/3;

% 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)
  Int=Int+(x(n2)-x(n0))*(c+0.5*b*(x(n2)+x(n0))+inv_3*a*(x(n2)^2+x(n2)*x(n0)+x(n0)^2));
  % 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)
  Int=Int+(x(n2)-x(n1))*(c+0.5*b*(x(n2)+x(n1))+inv_3*a*(x(n2)^2+x(n2)*x(n1)+x(n1)^2));
  % advance n
  n=n+1;
end