function retval = GaussInt(fn_,nodes,n)
% usage retval = GaussInt(fcn,nodes,n)
% description:  inputs function handle of one variable fcn and
% ordered vector nodes and returns an approximate value for the
% definite integral of fcn(x) on the interval enclosing the
% coordinates of nodes, calculated by using a n-point  Gaussian
% quadrature method between adjacent nodes. The function fcn
% must handle vector arguments correctly: evaluate each coordinate.
% Only 2<=n<=8 are implemented.  This routine is useful for integrating 
% pp functions exactly, e.g., for cubic splines use spline nodes and n=3.

% local variables:
% xabsc: Gaussian abscissas 
% t: translates of Gaussian abscissas and vals
% wts: Gaussian weights
% h: width of subinterval of Gauss integration
% k: index variables

switch n %decide on size of interpolation
  case 2
    xabsc = [ -0.577350269189626, 0.577350269189626];
    wts = [ 1.0, 1.0 ];
  case 3
    xabsc = [ -0.774596669241483 0 0.774596669241483 ];
    wts = [ 0.555555555555556 0.888888888888889 0.555555555555556 ];
  case 4
    xabsc = [  -0.861136311594052 -0.339981043584856 0.339981043584856 ...
0.861136311594052]; 
wts = [  0.347854845137454 0.652145154862546 ... 
0.652145154862546 0.347854845137454];
  case 5
    xabsc = [ -0.9061798459386639928, -0.5384693101056830910, ... 
0.,0.5384693101056830910, 0.9061798459386639928];      
    wts = [0.236926885056189, 0.4786286704993667, ...
          0.5688888888888887, 0.4786286704993664, ...
          0.2369268850561892];
  case 6
    xabsc = [-0.932469514203151 -0.661209386466264 -0.238619186083197 ...
0.238619186083197 0.661209386466264 0.932469514203151];
    wts = [  0.17132449237917 0.360761573048138 0.467913934572691 ...	
0.467913934572691 0.360761573048138 0.17132449237917];
  case 7
    xabsc = [-0.949107912342758 -0.741531185599394 -0.405845151377397 0 ...
0.405845151377397 0.741531185599394 0.949107912342758];
    wts = [ 0.129484966168879 0.279705391489274 0.381830050505119 ...
0.417959183673469 0.381830050505119 0.279705391489274 0.129484966168879];
  case 8
    xabsc = [  -0.960289856497535 -0.796666477413627 -0.52553240991633 ...
-0.183434642495649 0.183434642495649 ...
0.52553240991633 0.796666477413627 0.960289856497535];
    wts = [ 0.101228536290382 0.22238103445337 0.313706645877889 ...
0.362683783378362 0.362683783378362 ...
0.313706645877889 0.22238103445337 0.101228536290382];
  otherwise
    disp('This node number is not supported');
    return;
end
retval = 0;
nn = length(nodes);
for k = 1:nn-1
  t = nodes(k) + (nodes(k+1)-nodes(k))/2*(xabsc+1);
  t = feval(fn_,t);
  retval = retval + (nodes(k+1)-nodes(k))/2*wts*t';
end;