function [ theta, weight, x ] = glpairs ( n, k )

%*****************************************************************************80
%
%% GLPAIRS computes the K-th pair of an N-point Gauss-Legendre rule.
%
%  Discussion:
%
%    This routine is intended for cases were 100 < N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2016
%
%  Author:
%
%    Original C++ version by Ignace Bogaert.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Ignace Bogaert,
%    Iteration-free computation of Gauss-Legendre quadrature nodes and weights,
%    SIAM Journal on Scientific Computing,
%    Volume 36, Number 3, 2014, pages A1008-1026.
%
%  Parameters:
%
%    Input, integer N, the number of points in the given rule.
%    1 <= N.
%
%    Input, integer K, the index of the point to be returned.
%    1 <= K <= N.
%
%    Output, real THETA, WEIGHT, X, the theta coordinate, weight, 
%    and x coordinate of the point.
%
  if ( n < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GLPAIRS - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of N.\n' );
    error ( 'GLPAIRS - Fatal error!' );
  end

  if ( k < 1 || n < k )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GLPAIRS - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of K.\n' );
    error ( 'GLPAIRS - Fatal error!' );
  end

  if ( n < ( 2 * k - 1 ) )
    kcopy = n - k + 1;
  else
    kcopy = k;
  end
%
%  Get the Bessel zero.
%
  w = 1.0E+00 / ( n + 0.5E+00 );
  nu = besseljzero ( kcopy );
  theta = w * nu;
  y = theta ^ 2;
%
%  Get the asymptotic BesselJ(1,nu) squared.
%
  b = besselj1squared ( kcopy );
%
%  Get the Chebyshev interpolants for the nodes.
%
  sf1t = ((((( ...
    - 1.29052996274280508473467968379E-12 * y ...
    + 2.40724685864330121825976175184E-10 ) * y ... 
    - 3.13148654635992041468855740012E-08 ) * y ...
    + 0.275573168962061235623801563453E-05 ) * y ... 
    - 0.148809523713909147898955880165E-03 ) * y ...
    + 0.416666666665193394525296923981E-02 ) * y ...
    - 0.416666666666662959639712457549E-01;

  sf2t = ((((( ...
    + 2.20639421781871003734786884322E-09 * y ...
    - 7.53036771373769326811030753538E-08 ) * y ...
    + 0.161969259453836261731700382098E-05 ) * y ...
    - 0.253300326008232025914059965302E-04 ) * y ...
    + 0.282116886057560434805998583817E-03 ) * y ...
    - 0.209022248387852902722635654229E-02 ) * y ...
    + 0.815972221772932265640401128517E-02;

  sf3t = ((((( ...
    - 2.97058225375526229899781956673E-08 * y ...
    + 5.55845330223796209655886325712E-07 ) * y ...  
    - 0.567797841356833081642185432056E-05 ) * y ... 
    + 0.418498100329504574443885193835E-04 ) * y ...
    - 0.251395293283965914823026348764E-03 ) * y ...
    + 0.128654198542845137196151147483E-02 ) * y ...
    - 0.416012165620204364833694266818E-02;
%
%  Get the Chebyshev interpolants for the weights.
%
  wsf1t = (((((((( ...
    - 2.20902861044616638398573427475E-14 * y ...
    + 2.30365726860377376873232578871E-12 ) * y ...
    - 1.75257700735423807659851042318E-10 ) * y ...
    + 1.03756066927916795821098009353E-08 ) * y ...
    - 4.63968647553221331251529631098E-07 ) * y ...
    + 0.149644593625028648361395938176E-04 ) * y ...
    - 0.326278659594412170300449074873E-03 ) * y ...
    + 0.436507936507598105249726413120E-02 ) * y ...
    - 0.305555555555553028279487898503E-01 ) * y ...
    + 0.833333333333333302184063103900E-01;

  wsf2t = ((((((( ...
    + 3.63117412152654783455929483029E-12 * y ...
    + 7.67643545069893130779501844323E-11 ) * y ...
    - 7.12912857233642220650643150625E-09 ) * y ...
    + 2.11483880685947151466370130277E-07 ) * y ...
    - 0.381817918680045468483009307090E-05 ) * y ...
    + 0.465969530694968391417927388162E-04 ) * y ...
    - 0.407297185611335764191683161117E-03 ) * y ...
    + 0.268959435694729660779984493795E-02 ) * y ...
    - 0.111111111111214923138249347172E-01;

  wsf3t = ((((((( ...
    + 2.01826791256703301806643264922E-09 * y ...
    - 4.38647122520206649251063212545E-08 ) * y ...
    + 5.08898347288671653137451093208E-07 ) * y ...
    - 0.397933316519135275712977531366E-05 ) * y ...
    + 0.200559326396458326778521795392E-04 ) * y ...
    - 0.422888059282921161626339411388E-04 ) * y ...
    - 0.105646050254076140548678457002E-03 ) * y ...
    - 0.947969308958577323145923317955E-04 ) * y ...
    + 0.656966489926484797412985260842E-02;
%
%  Refine with the paper expansions.
%
  nuosin = nu / sin ( theta );
  bnuosin = b * nuosin;
  winvsinc = w * w * nuosin;
  wis2 = winvsinc * winvsinc;
% 
%  Finally compute the node and the weight.
%
  theta = w * ( nu + theta * winvsinc ...
    * ( sf1t + wis2 * ( sf2t + wis2 * sf3t ) ) );
  deno = bnuosin + bnuosin * wis2 * ( wsf1t + wis2 * ( wsf2t + wis2 * wsf3t ) );
  weight = ( 2.0E+00 * w ) / deno;

  if ( n < ( 2 * k - 1 ) )
    theta = pi - theta;
  end

  x = cos ( theta );

  return
end
