function [x,wt] = GaussPts(nodes,n)
% usage: [x,wt] = GaussPts(nodes,n)
% description:  Given input of ordered vector nodes and
% degree n, returns returns the abscissas x and weights
% w for an n-point Gaussian quadrature method on on each
% subinterval of adjacent nodes. Only 2<=n<=8 are
% implemented. Degree n is optional and default is n=5.

% 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 ~= 2)
  n = 5;
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
nn = length(nodes);
x = zeros(1,n*(nn-1));
wt = repmat(wts,1,nn-1);
for k = 1:nn-1
  x(1,1+(k-1)*n:k*n) = nodes(k) + (nodes(k+1)-nodes(k))/2*(xabsc+1);
end;