#! /usr/bin/env python3
#
def exercise6 ( ):

#*****************************************************************************80
#
## exercise6() uses gauss_quad() to estimate the integral of a function.
#
  from gauss_quad import gauss_quad
  import numpy as np

  print ( 'exercise 6:' )

  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__' ):
  exercise6 ( )