function [ xtab, weight ] = gegenbauer_compute ( order, alpha )

%% GEGENBAUER_COMPUTE computes the abscissa and weights for Gauss-Gegenbauer quadrature.
%
%  Discussion:
%
%    The weight function is w(x) = (1-X^2)^ALPHA.
%
%    The integral to approximate is:
%
%      Integral ( -1 <= X <= 1 ) (1-X^2)^ALPHA * F(X) dX
%
%    The quadrature formula is:
%
%      Sum ( 1 <= I <= NORDER ) WEIGHT(I) * F ( XTAB(I) )
%
%    Thanks to Janiki Raman for pointing out a problem in an earlier
%    version of the code that occurred when ALPHA was -0.5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2008
%
%  Author:
%
%    FORTRAN77 original version by Arthur Stroud, Don Secrest.
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, integer ORDER, the order of the quadrature rule.
%
%    Input, real ALPHA, the exponent of (1-X^2) in the weight.  
%    -1.0 < ALPHA is required.
%
%    Output, real XTAB(NORDER), the abscissas.
%
%    Output, real WEIGHT(NORDER), the weights.
%

%
%  Check ORDER.
%
  if ( order < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GEGENBAUER_COMPUTE - Fatal error!\n' );
    fprintf ( 1, '  1 <= ORDER is required.\n' );
    error ( 'GEGENBAUER_COMPUTE - Fatal error!' );
  end
%
%  Check ALPHA.
%
  if ( alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GEGENBAUER_COMPUTE - Fatal error!\n' );
    fprintf ( 1, '  -1.0 < ALPHA is required.\n' );
    error ( 'GEGENBAUER_COMPUTE - Fatal error!' );
  end
%
%  Set the recursion coefficients.
%
  c(1) = 0.0;
  if ( 2 <= order )
    c(2) = 1.0 / ( 2.0 * alpha + 3.0 );
  end

  for i = 3 : order

    c(i) = ( i - 1 ) ...
      * ( alpha + alpha + i - 1 ) / ...
      ( ( alpha + alpha + 2 * i - 1 ) ...
      * ( alpha + alpha + 2 * i - 3 ) );

  end

  delta = r8_gamma ( alpha         + 1.0 ) ...
        * r8_gamma (         alpha + 1.0 ) ...
        / r8_gamma ( alpha + alpha + 2.0 );

  cc = delta * 2.0^( 2.0 * alpha + 1.0 ) * prod ( c(2:order) );

  for i = 1 : order

    if ( i == 1 )

      an = alpha / order;

      r1 = ( 1.0 + alpha ) * ( 2.78 / ( 4.0 + order * order ) ...
        + 0.768 * an / order );

      r2 = 1.0 + 2.44 * an + 1.281 * an^2;

      x = ( r2 - r1 ) / r2;

    elseif ( i == 2 )

      r1 = ( 4.1 + alpha ) / ...
        ( ( 1.0 + alpha ) * ( 1.0 + 0.156 * alpha ) );

      r2 = 1.0 + 0.06 * ( order - 8.0 ) * ( 1.0 + 0.12 * alpha ) / order;

      r3 = 1.0 + 0.012 * alpha * ...
        ( 1.0 + 0.25 * abs ( alpha ) ) / order;

      x = x - r1 * r2 * r3 * ( 1.0 - x );

    elseif ( i == 3 )

      r1 = ( 1.67 + 0.28 * alpha ) / ( 1.0 + 0.37 * alpha );

      r2 = 1.0 + 0.22 * ( order - 8.0 ) / order;

      r3 = 1.0 + 8.0 * alpha / ( ( 6.28 + alpha ) * order * order );

      x = x - r1 * r2 * r3 * ( xtab(1) - x );

    elseif ( i < order - 1 )

      x = 3.0 * xtab(i-1) - 3.0 * xtab(i-2) + xtab(i-3);

    elseif ( i == order - 1 )

      r1 = ( 1.0 + 0.235 * alpha ) / ( 0.766 + 0.119 * alpha );

      r2 = 1.0 / ( 1.0 + 0.639 * ( order - 4.0 ) ...
        / ( 1.0 + 0.71 * ( order - 4.0 ) ) );

      r3 = 1.0 / ( 1.0 + 20.0 * alpha / ( ( 7.5 + alpha ) * order * order ) );

      x = x + r1 * r2 * r3 * ( x - xtab(i-2) );

    elseif ( i == order )

      r1 = ( 1.0 + 0.37 * alpha ) / ( 1.67 + 0.28 * alpha );

      r2 = 1.0 / ( 1.0 + 0.22 * ( order - 8.0 ) / order );

      r3 = 1.0 / ( 1.0 + 8.0 * alpha / ( ( 6.28 + alpha ) * order * order ) );

      x = x + r1 * r2 * r3 * ( x - xtab(i-2) );

    end

    [ x, dp2, p1 ] = gegenbauer_root ( x, order, alpha, c );

    xtab(i) = x;
    weight(i) = cc / ( dp2 * p1 );

  end
%
%  Reverse the order of the values.
%
  xtab(1:order) = xtab(order:-1:1);
  weight(1:order) = weight(order:-1:1);