function [xgauss,wgauss] = GaussGrid(xknots,ngauss)
%[xgauss,wgauss] = GaussGrid(xknots,ngauss)

nx = length(xknots);

% check lengths

[xg,wg] = GaussianPts(ngauss); % gaussian quadrature points on each interval for [-1 1]
% (GaussianPts function taken from one of S. Brunner's files)

wgi = zeros(ngauss,nx-1); % allocate
xgi = zeros(ngauss,nx-1);
for ix = 1:nx-1
    % modify gaussian parameters for this integration interval
    dig = diff(xknots(ix:ix+1));
    wgi(:,ix) = dig/2 *wg; % modified weights
    xgi(:,ix) = dig/2*xg + mean(xknots(ix:ix+1));
end

xgauss = reshape(xgi,numel(xgi),1);
wgauss = reshape(wgi,numel(wgi),1);