#! /usr/bin/env python
#
def gegenbauer_ek_compute ( n, alpha ):

#*****************************************************************************80
#
## GEGENBAUER_EK_COMPUTE: Elhay-Kautsky method for Gauss-Gegenbauer rule.
#
#  Discussion:
#
#    The integral:
#
#      integral ( -1 <= x <= 1 ) (1-x^2)^(alpha-1/2) * f(x) dx
#
#    The quadrature rule:
#
#      sum ( 1 <= i <= n ) w(i) * f ( x(i) )
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    24 November 2015
#
#  Author:
#
#    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Sylvan Elhay, Jaroslav Kautsky,
#    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
#    Interpolatory Quadrature,
#    ACM Transactions on Mathematical Software,
#    Volume 13, Number 4, December 1987, pages 399-415.
#
#  Parameters:
#
#    Input, integer N, the order.
#
#    Input, real ALPHA, determines the exponent of (1-X^2)^(ALPHA-1/2).
#
#    Output, real X(N), the abscissas.
#
#    Output, real W(N), the weights.
#
  import numpy as np
  from imtqlx import imtqlx
  from r8_gamma import r8_gamma
#
#  Define the zero-th moment.
#
  zemu = 2.0 ** ( 2.0 * alpha + 1.0 ) \
    * r8_gamma ( alpha + 1.0 ) \
    * r8_gamma ( alpha + 1.0 ) \
    / r8_gamma ( 2.0 * alpha + 2.0 )
#
#  Define the Jacobi matrix.
#
  x = np.zeros ( n )

  bj = np.zeros ( n )

  bj[0] = 4.0 * ( alpha + 1.0 ) ** 2 \
    / ( ( 2.0 * alpha + 3.0 ) * ( 2.0 * alpha + 2.0 ) ** 2 )

  for i in range ( 2, n ):
    abi = 2.0 * ( alpha + i )
    bj[i-1] = 4.0 * i * ( alpha + i ) ** 2 * ( 2.0 * alpha + i ) \
      / ( ( abi - 1.0 ) * ( abi + 1.0 ) * abi * abi )

  for i in range ( 0, n ):
    bj[i] = np.sqrt ( bj[i] )

  w = np.zeros ( n )

  w[0] = np.sqrt ( zemu )
#
#  Diagonalize the Jacobi matrix.
#
  x, w = imtqlx ( n, x, bj, w )

  for i in range ( 0, n ):
    w[i] = w[i] ** 2

  return x, w

def gegenbauer_ek_compute_test ( ):

#*****************************************************************************80
#
## GEGENBAUER_EK_COMPUTE_TEST compares GEGENBAUER_EK_COMPUTE.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    24 November 2015
#
#  Author:
#
#    John Burkardt
#
  import platform

  alpha = 0.5

  print ( '' )
  print ( 'GEGENBAUER_EK_COMPUTE_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  GEGENBAUER_EK_COMPUTE computes Gauss-Gegenbauer rules' )
  print ( '' )
  print ( '  Abscissas and weights for a generalized Gauss Gegenbauer rule' )
  print ( '  with ALPHA = %g' % ( alpha ) )
  print ( '  Integration interval is [-1,+1]' )

  for n in range ( 1, 11 ):

    x, w = gegenbauer_ek_compute ( n, alpha )

    print ( '' )

    for i in range ( 0, n ):
      print ( '  %2d  %24.16g  %24.16g' % ( i, w[i], x[i] ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'GEGENBAUER_EK_COMPUTE_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  gegenbauer_ek_compute_test ( )
  timestamp ( )