function [points,poids]=points_poids_gauss_diago(c,n)

% points_poids_gauss_diago : calcul des poids et des points de Gauss par diagonalisation
% 
% *********************************************************
%
% [points,poids]=points_poids_gauss_diago(c,n) renvoie les points D_i et les points x_i
%   				  des formules d'intégration à l+1 points pour 0<=l<=n, pour 
%                l'une des quatre méthodes : Legendre, Tchebytchev, Hermite ou Laguerre,
%                calculés par diagonalisation d'une matrice tridiagonale.   
%
%	variables d'entrées : 
%    *  n : entier naturel.
%     * c définit la méthode : 
%            c=1 : méthode de Gauss-Legendre,
%            c=2 : méthode de Gauss-Tchebychev,
%            c=3 : méthode de Gauss-Hermite,
%            c=4 : méthode de Gauss-Laguerre;
%
%	variables de sortie
%     * points : tableau de taille (n+1)*(n+1) tel que 
%                pour tout l dans {0,...,n}, points(l+1,1:l+1) contient
%                les l+1 points x_i de la formule (à l+1 points)
%     * poids  : tableau de taille (n+1)*(n+1) tel que 
%                pour tout l dans {0,...,n}, poids(l+1,1:l+1) contient
%                les l+1 poids D_i de la formule (à l+1 points)
%
% 
% 
% ************ Fonctions auxiliaires utilisées ************
%
%       aucune
%
% *********************************************************


% Contrôles d'entrée
% nombre d'arguments
if nargin~=2
   error('nombre d''arguments de la fonction incorrect');
end
% autres tests éventuels 
if (fix(n)~=n) | (n<0)
   error('l''entier n doit être un entier naturel');
end

% Corps d'algorithme.
points=zeros(n+1,n+1);
poids=zeros(n+1,n+1);
switch c
   case 1
      points(1,1)=0;
      poids(1,1)=2;
   case 2
      points(1,1)=0;
      poids(1,1)=pi;
   case 3
      points(1,1)=0;
      poids(1,1)=sqrt(pi);
   case 4
      points(1,1)=1;
      poids(1,1)=1;
end
if (n>=1)
   for k=1:n
      if (c<=3)
         B=zeros(1,k+1);
      else
         B=2*(0:k)+1;
      end
      switch c
         case 1
            gamma=1./(2*sqrt(1-1./(4*(1:k).^2)));
         case 2
            gamma=[1/sqrt(2) ones(1,k-1)/2];
         case 3
            gamma=sqrt((1:k)/2);
         case 4
            gamma=(1:k);
      end
      auxi=gallery('tridiag',gamma,B,gamma);
      if (k<=4) 
         [Q,R]=eigs(auxi,k+1);
      else
         [Q,R]=eig(full(auxi));
      end
      R=diag(R);
      U=(sqrt(sum(Q.^2)));
      switch c
         case 1
            auxik=2;
         case 2
            auxik=pi;
         case 3
            auxik=sqrt(pi);
         case 4
            auxik=1;
      end
      Q=auxik*(Q(1,:)./U).^2;
      [auxipoints,ind]=sort((flipud(R))');
      auxi=fliplr(Q);
      auxipoids=auxi(ind);
      if (c<=3) 
         p=fix(k/2);
         if ~(rem(k,2))
            auxi=(auxipoints(p+2:k+1)-auxipoints(p:-1:1))/2;      
            auxipoints(p+2:k+1)=auxi;
            auxipoints(1:p)=-fliplr(auxi);
            auxipoints(p+1)=0;
            auxi=(auxipoids(p+2:k+1)+auxipoids(p:-1:1))/2;      
            auxipoids(p+2:k+1)=auxi;
            auxipoids(1:p)=fliplr(auxi);
         else
            auxi=(auxipoints(p+2:k+1)-auxipoints(p+1:-1:1))/2;      
            auxipoints(p+2:k+1)=auxi;
            auxipoints(1:p+1)=-fliplr(auxi);
            auxi=(auxipoids(p+2:k+1)+auxipoids(p+1:-1:1))/2;      
            auxipoids(p+2:k+1)=auxi;
            auxipoids(1:p+1)=fliplr(auxi);
         end
      end
      points(k+1,1:k+1)=auxipoints;
      poids(k+1,1:k+1)=auxipoids;
   end
end

