function retval = GaussInt(fn_,nodes,n)
% usage retval = GaussInt(fcn,nodes,n)
% description:  inputs a function of one variable fcn, ordered
% vector nodes and degree n, 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=2.
% Function name fcn is a string, so must be quoted, or use a handle.
% Degree n is optional and default is n=2.

% 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

if (nargin ~= 3)
  n = 2;
end
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
    error('This node number is not supported');
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;