#! /usr/bin/env python3
#
def gauss_quad ( func, xn, nq ):

#*****************************************************************************80
#
## gauss_quad() applies nq-point Gauss rule to integral of func() over a mesh xn.
#
#  Input:
#    func(): a function to be integrated over the mesh.
#    xn(): points which define a mesh.
#    nq: the number of quadrature points to use.
#  Output:
#    q: the integral as estimated by Gauss quadrature.
#
  from numpy.polynomial.legendre import leggauss
  import numpy as np

  xq, wq = leggauss ( nq )

  q = 0.0
  nn = len ( xn )

  for i in range ( 0, nn - 1 ):
    xab = ( xn[i] * ( 1.0 - xq ) + xn[i+1] * ( xq + 1.0 ) ) / 2.0
    wab = ( xn[i+1] - xn[i] ) * wq / 2.0
    fab = func ( xab )
    qab = np.sum ( wab * fab )
    q = q + qab

  return q

def gauss_quad_test ( ):
  import numpy as np

  print ( 'gauss_quad_test:' )

  func = humps
  a = 0.0
  b = 2.0
  exact = humps_antideriv ( b ) - humps_antideriv ( a )
  for nn in [ 2, 4, 8, 16, 32 ]:
    xn = np.linspace ( 0.0, 2.0, nn )
    print ( "" )
    for nq in [ 1, 2, 4, 8, 16 ]:
      q = gauss_quad ( func, xn, nq )
      err = abs ( q - exact )
      print ( "(nq,nn) = ", nq, nn, "  err = ", err )

  return

def humps_antideriv ( x ):

#*****************************************************************************80
#
## humps_antideriv evaluates the antiderivative of the humps function.
#
#  Discussion:
#
#    y = 1.0 / ( ( x - 0.3 )**2 + 0.01 ) \
#      + 1.0 / ( ( x - 0.9 )**2 + 0.04 ) \
#      - 6.0
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    24 September 2021
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real x(): the argument.
#
#  Output:
#
#    real y(): the value of the antiderivative at x.
#
  import numpy as np

  y = 10.0 * np.arctan ( 10.0 * ( x - 0.3 ) ) \
    +  5.0 * np.arctan (  5.0 * ( x - 0.9 ) ) \
    - 6.0 * x

  return y

def humps ( x ):

#*****************************************************************************80
#
## humps evaluates the humps function.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 August 2019
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real x(): the evaluation points.
#
#  Output:
#
#    real y(): the function values.
#
  y = 1.0 / ( ( x - 0.3 )**2 + 0.01 ) \
    + 1.0 / ( ( x - 0.9 )**2 + 0.04 ) \
    - 6.0

  return y

if ( __name__ == '__main__' ):
  gauss_quad_test ( )