function [xg, wg] = GaussianPts(kg)
% function [xg, wg] = GaussianPts(kg)
% Provides Gaussian integration points and weights.
%
% Integration points: xg = Zeros of Legendre polynomial L_n of
%                         degree n = kg
% Weights: wg = 2/(1-xg^2)(d/dx L_n(xg))^2

Lc = Legendre(kg);
xg = sort(roots(Lc));
wg = 2./((1-xg.^2).*polyval(polyder(Lc), xg).^2);
 
% % use lookuptable
% switch kg
%     case 1
%         xg = 1; 
%         wg=2;
%     case 2
%         xg = 1/sqrt(3)* [-1 1]; 
%         wg = [1 1];
%     case 3
%         xg = [-sqrt(15)/5,0,sqrt(15)/5]; wg = [5/9 8/9 5/9];
%     case 4
%         xg = [sqrt((3+2*sqrt(6/5))/7),sqrt((3-2*sqrt(6/5))/7)]*[-1 0 0 1;0 -1 1 0];
%         wg = [(18-sqrt(30))/36,(18+sqrt(30))/36]*[1 0 0 1;0 1 1 0];
% end
%------------------------------

function Lc = Legendre(n)

% Provides coefficients of Legendre polynomial of degree n
%
% Ln(x) = (1/norm) d^n/dx^n (x^2-1)^n
%
% with norm = 2^n * n!

pc = [1 0 -1];
if n == 0,
  Lc = 1;
else
  Lc = pc;
end
for il = 1:n-1,
  Lc = conv(Lc, pc);
end

for il = 1:n,
  Lc = polyder(Lc);
end

Lc = Lc / (2^n * factorial(n));