## Copyright (C) 2009 Carlo de Falco
## 
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA

## usage:
## I = integ1samp (x, y, order)

## Author: Carlo de Falco <carlo@guglielmo.local>
## Created: 2009-06-07

function I = integ1samp (x, y, order)

  sample_length = length(x);
  if (length (y) != sample_length)
    error ("x and y must be of the same size")
  endif
  
  if (mod (sample_length-1, order) != 0)
    error ("the number of sub intervals must be a multiple of the interpolation polynomial order")
  endif
  
  X = __rearrange__ (x, order, sample_length);
  Y = __rearrange__ (y, order, sample_length);

  I = 0;
  for ii=1:rows (X)
    c = polyfit (X(ii, :), Y(ii, :), order);
    c = polyint (c);
    I += diff(polyval (c, X(ii, [1 end])));
  endfor

endfunction

function X = __rearrange__ (x, order, sample_length)
  X = zeros ((sample_length-1)/order, order+1);
  X([1 end])  = x([1 end]);
  X(2:end, 1) = X(1:end-1, end) = x(order+1:order:end-order);
  x(1:order:end) = [];
  X(:,2:end-1) = reshape (x, order-1, [])';
endfunction

%!demo
%! f = @(x) 1./(1+25*x.^2);
%! Iex = quad (f, -1, 1);
%! printf ("  N \t\t error:\n \t\t order 2 \t\t order 3 \t\t order N-1\n")
%! for ii = 1:7
%!  N(ii) = 3 * 2 * ii + 1;
%!  x = linspace (-1, 1, N(ii));
%!  I2 = integ1samp (x, f(x), 2);
%!  I3 = integ1samp (x, f(x), 3);
%!  IN = integ1samp (x, f(x), N(ii)-1);
%!  printf ("% 3d \t\t %03.2e \t\t %03.2e \t\t %03.2e\n", N(ii), abs(I2-Iex), abs(I3-Iex), abs(IN-Iex))
%! endfor
%!demo
%! f = @(x) exp (-5 * x.^2);
%! Iex = quad (f, -pi, pi);
%! printf ("  N \t\t error:\n \t\t order 2 \t\t order 3 \t\t order N-1\n")
%! for ii = 1:7
%!  N(ii) = 3 * 2 * ii + 1;
%!  x = linspace (-pi, pi, N(ii));
%!  I2 = integ1samp (x, f(x), 2);
%!  I3 = integ1samp (x, f(x), 3);
%!  IN = integ1samp (x, f(x), N(ii)-1);
%!  printf ("% 3d \t\t %03.2e \t\t %03.2e \t\t %03.2e\n", N(ii), abs(I2-Iex), abs(I3-Iex), abs(IN-Iex))
%! endfor
%!demo
%! f = @(x) x.^5;
%! Iex = quad (f, 0, 1);
%! printf ("  N \t\t error:\n \t\t order 2 \t\t order 3 \t\t order N-1\n")
%! for ii = 1:7
%!  N(ii) = 3 * 2 * ii + 1;
%!  x = linspace (0, 1, N(ii));
%!  I2 = integ1samp (x, f(x), 2);
%!  I3 = integ1samp (x, f(x), 3);
%!  IN = integ1samp (x, f(x), N(ii)-1);
%!  printf ("% 3d \t\t %03.2e \t\t %03.2e \t\t %03.2e\n", N(ii), abs(I2-Iex), abs(I3-Iex), abs(IN-Iex))
%! endfor
%!demo
%! f = @(x) abs(mod(x, 1) - 1/2);
%! Iex = quad (f, 0, 3);
%! printf ("  N \t\t error:\n \t\t order 1 \t\t order 3 \t\t order N-1\n")
%! for ii = 1:7
%!  N(ii) = 2 * 3 * ii + 1;
%!  x = linspace (0, 3, N(ii));
%!  I2 = integ1samp (x, f(x), 1);
%!  I3 = integ1samp (x, f(x), 3);
%!  IN = integ1samp (x, f(x), N(ii)-1);
%!  printf ("% 3d \t\t %03.2e \t\t %03.2e \t\t %03.2e\n", N(ii), abs(I2-Iex), abs(I3-Iex), abs(IN-Iex))
%! endfor
