function gausleg = lgndre ( NP )

%*****************************************************************************80
%
%MATLAB version of SUBROUTNE LGNDRE on page 316 of conte/deboor
% returns, in the structure  gausleg , the points and weights of the
% Gauss-LEgendre quadrature involving  NP  points.
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
  if NP>10
    disp('LGNDRE only knows the rules for  NP<= 10')
    return
  end

  point = zeros(1,NP); 
  weight = point;

  switch NP
  case 1
    point = 0; 
    weight = 1;
  case 2
    point(2) = 0.5773502691896257; 
    weight(2) = 1;
  case 3
    point(3) = 0.7745966692414834; 
    weight(2:3) = [5,8]/9;
  case 4
    point(3:4) = [0.3399810435848563, 0.8611363115940526];
    weight(3:4) = [0.6521451548625461, 0.3478548451374538];
  case 5
    point(4:5) = [0.5384693101056831, 0.9061798459386640];
    weight(3:5) = [0.5688888888888889, 0.4786286704993665, 0.2369268850561891];
  case 6
    point(4:6) = [0.2386191860831969, 0.6612093864662645, 0.9324695142031521];
    weight(4:6) = [0.4679139345726910, 0.3607615730481386, 0.1713244923791704];
  case 7
    point(5:7) = [ ...
      0.405845151377397166906606412077, ...
      0.741531185599394439863864773281, ...
      0.949107912342758524526189684048 ];
    weight(4:7) = [ ...
      0.417959183673469387755102040816, ...
      0.381830050505118944950369775489, ...
      0.279705391489276667901467771424, ...
      0.129484966168869693270611432679 ];
  case 8
    point(5:8) = [ ...
      0.183434642495649804939476142360, ...
      0.525532409916328985817739049189, ...
      0.796666477413626739591553936476, ...
      0.960289856497536231683560868569 ];
    weight(5:8) = [ ...
      0.362683783378361982965150449277, ...
      0.313706645877887287337962201987, ...
      0.222381034453374470544355994426, ...
      0.101228536290376259152531354310 ];
  case 9
    point(6:9) = [ ...
      0.324253423403808929038538015, ...
      0.613371432700590397308702039, ...
      0.836031107326635794299429788, ...
      0.968160239507626089835576203 ];
    weight(5:9) = [ ...
      0.33023935500125976316452506929, ...
      0.31234707704000284006863040658, ...
      0.26061069640293546231874286942, ...
      0.18064816069485740405847203124, ...
      0.081274388361574411971892158111 ];
  case 10
    point(6:10) = [ ...
      0.148874338981631210884826001, ...
      0.433395394129247190799265943, ...
      0.679409568299024406234327365, ...
      0.865063366688984510732096688, ...
      0.973906528517171720077964012 ];
    weight(6:10) = [ ...
      0.29552422471475287017389299465, ...
      0.26926671930999635509122692157, ...
      0.21908636251598204399553493423, ...
      0.14945134915058059314577633966, ...
      0.066671344308688137593568809893 ] ;
  end

  ii = 1:NP/2; 
  point(ii) = -point(NP+1-ii);
  weight(ii) = weight(NP+1-ii);
  gausleg = struct ( 'points', point, 'weights', weight ); 

  return
end
