#! /usr/bin/env python
#
def gegenbauer_cc1 ( n, lam, f ):

#*****************************************************************************80
#
## GEGENBAUER_CC1 estimates the Gegenbauer integral of a function.
#
#  Discussion:
#
#     value = integral ( -1 <= x <= + 1 ) ( 1 - x^2 )^(lambda-1/2) * f(x) dx
#
#     The approximation uses the practical abscissas, that is, the extreme
#     points of Tn(x).
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    14 January 2016
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    D B Hunter, H V Smith,
#    A quadrature formula of Clenshaw-Curtis type for the Gegenbauer weight function,
#    Journal of Computational and Applied Mathematics,
#    Volume 177, 2005, pages 389-400.
#
#  Parameters:
#
#    Input, integer N, the number of points to use.
#    1 <= N.
#
#    Input, real LAM, used in the exponent of (1-x^2).
#    -0.5 < LAM.
#
#    Input, real F(x), the handle for the function to be integrated with the
#    Gegenbauer weight.
#
#    Output, real WEIGHT, the estimate for the Gegenbauer integral of F.
#
  import numpy as np
  from scipy.special import gamma
  from chebyshev_even1 import chebyshev_even1

  value = 0.0

  s = ( n // 2 )
  sigma = ( n % 2 )

  a2 = chebyshev_even1 ( n, f )

  rh = s
  u = 0.5 * ( sigma + 1.0 ) * a2[rh]
  for rh in range ( s - 1, 0, -1 ):
    u = ( rh - lam ) / ( rh + lam + 1 ) * u + a2[rh]
  u = - lam * u / ( lam + 1.0 ) + 0.5 * a2[0]

  value = gamma ( lam + 0.5 ) * np.sqrt ( np.pi ) * u / gamma ( lam + 1.0 );

  return value

def gegenbauer_cc1_test ( ):

#*****************************************************************************80
#
## GEGENBAUER_CC1_TEST tests GEGENBAUER_CC1.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    14 January 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from scipy.special import gamma
  from scipy.special import jv

  print ( '' )
  print ( 'GEGENBAUER_CC1_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  GEGENBAUER_CC1 estimates the Gegenbauer integral of' )
  print ( '  a function f(x) using a Clenshaw-Curtis type approach' )
  print ( '  based on the extreme points of Tn(x).' )

  lam = 0.75
  a = 2.0
  n = 6

  value = gegenbauer_cc1 ( n, lam, f )

  print ( '' )
  print ( '  Value = %g' % ( value ) )

  exact = gamma ( lam + 0.5 ) * np.sqrt ( np.pi ) * jv ( lam, a ) / ( 0.5 * a ) ** lam

  print ( '  Exact = %g' % ( exact ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'GEGENBAUER_CC1_TEST' )
  print ( '  Normal end of execution.' )
  return

def f ( x ):

#*****************************************************************************80
#
## F defines the function whose Fourier coefficients are desired.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    14 January 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, real X, the argument.
#
#    Output, real VALUE, the value.
#
  import numpy as np

  a = 2.0
  value = np.cos ( a * x )

  return value

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