function [xGL,wGL] = g_data(n)
%
% [xGL,wGL] = g_data(n)
% 
% Returns Gauss-Legendre nodes and weights for n points using 
% m-files of Walter Gautschi.
%
% Acknowlegdement
%
% File supplied by O. Ernst, 200?
%
% SFemLib function 02/10/07

if n == 1,
    xGL = 0;
    wGL = 1;
    return;
end

p = n-1;
ab = r_jacobi(p+1,0,0);
xw = gauss(p+1,ab);
xGL = xw(:,1);
wGL = xw(:,2);

return

%--------------------------------------------------------------------------
% 
% The following utility functions, obtained from Walter Gautschi's web site
% 
% http://www.cs.purdue.edu/archives/2002/wxg/codes/ ,
% 
% compute recurrence coefficients and Gaussian quadrature rules.
%
%--------------------------------------------------------------------------

% R_JACOBI Recurrence coefficients for monic Jacobi polynomials.
%
%    ab=R_JACOBI(n,a,b) generates the first n recurrence 
%    coefficients for monic Jacobi polynomials with parameters 
%    a and b. These are orthogonal on [-1,1] relative to the
%    weight function w(t)=(1-t)^a(1+t)^b. The n alpha-coefficients
%    are stored in the first column, the n beta-coefficients in
%    the second column, of the nx2 array ab. The call ab=
%    R_JACOBI(n,a) is the same as ab=R_JACOBI(n,a,a) and
%    ab=R_JACOBI(n) the same as ab=R_JACOBI(n,0,0).
%
%    Supplied by Dirk Laurie, 6-22-1998; edited by Walter
%    Gautschi, 4-4-2002.
%
function ab=r_jacobi(N,a,b)
if nargin<2, a=0; end;  if nargin<3, b=a; end
if((N<=0)|(a<=-1)|(b<=-1)) error('parameter(s) out of range'), end
nu=(b-a)/(a+b+2);
mu=2^(a+b+1)*gamma(a+1)*gamma(b+1)/gamma(a+b+2);
if N==1, ab=[nu mu]; return, end 
N=N-1; n=1:N; nab=2*n+a+b;
A=[nu (b^2-a^2)*ones(1,N)./(nab.*(nab+2))];
n=2:N; nab=nab(n);
B1=4*(a+1)*(b+1)/((a+b+2)^2*(a+b+3));
B=4*(n+a).*(n+b).*n.*(n+a+b)./((nab.^2).*(nab+1).*(nab-1));
ab=[A' [mu; B1; B']];

%----------------------------------------------------------------------

% LOBATTO Gauss-Lobatto quadrature rule.
%
%    Given a weight function w encoded by the (n+2)x2 array ab of
%    the first n+2 recurrence coefficients for the associated
%    orthogonal polynomials, the first column of ab containing the
%    n+2 alpha-coefficients and the second column the n+2 beta-
%    coefficients, the call xw=LOBATTO(n,ab,endl,endr) generates 
%    the nodes and weights xw of the (n+2)-point Gauss-Lobatto 
%    quadrature rule for the weight function w having two
%    prescribed nodes endl, endr (typically the left and right end
%    points of the support interval of w, or points to the left
%    resp. to the right therof). The n+2 nodes, in increasing 
%    order, are stored in the first column, the n+2 corresponding
%    weights in the second column, of the (n+2)x2 array xw.
%        
%    For Jacobi weight functions, see also LOBATTO_JACOBI.
%
function xw=lobatto(N,ab,endl,endr)
N0=size(ab,1); if N0<N+2, error('input array ab too short'), end
ab0=ab;
p0l=0; p0r=0; p1l=1; p1r=1;
for n=1:N+1
  pm1l=p0l; p0l=p1l; pm1r=p0r; p0r=p1r;
  p1l=(endl-ab0(n,1))*p0l-ab0(n,2)*pm1l;
  p1r=(endr-ab0(n,1))*p0r-ab0(n,2)*pm1r;
end
det=p1l*p0r-p1r*p0l;
ab0(N+2,1)=(endl*p1l*p0r-endr*p1r*p0l)/det;
ab0(N+2,2)=(endr-endl)*p1l*p1r/det;
xw=gauss(N+2,ab0);

%----------------------------------------------------------------------

% GAUSS Gauss quadrature rule.
%
%    Given a weight function w encoded by the nx2 array ab of the 
%    first n recurrence coefficients for the associated orthogonal
%    polynomials, the first column of ab containing the n alpha-
%    coefficients and the second column the n beta-coefficients, 
%    the call xw=GAUSS(n,ab) generates the nodes and weights xw of
%    the n-point Gauss quadrature rule for the weight function w.
%    The nodes, in increasing order, are stored in the first 
%    column, the n corresponding weights in the second column, of
%    the nx2 array xw.
%
function xw=gauss(N,ab)
N0=size(ab,1); if N0<N, error('input array ab too short'), end
J=zeros(N);
for n=1:N, J(n,n)=ab(n,1); end
for n=2:N
  J(n,n-1)=sqrt(ab(n,2));
  J(n-1,n)=J(n,n-1);
end
[V,D]=eig(J);
[D,I]=sort(diag(D));
V=V(:,I);
xw=[D ab(1,2)*V(1,:)'.^2];