#! /usr/bin/env python
#
def ball01_monomial_integral ( e ):

#*****************************************************************************80
#
## BALL01_MONOMIAL_INTEGRAL returns monomial integrals in the unit ball.
#
#  Discussion:
#
#    The integration region is 
#
#      X^2 + Y^2 + Z^2 <= 1.
#
#    The monomial is F(X,Y,Z) = X^E(1) * Y^E(2) * Z^E(3).
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    21 June 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Gerald Folland,
#    How to Integrate a Polynomial Over a Sphere,
#    American Mathematical Monthly,
#    Volume 108, Number 5, May 2001, pages 446-448.
#
#  Parameters:
#
#    Input, integer E(3), the exponents of X, Y and Z in the 
#    monomial.  Each exponent must be nonnegative.
#
#    Output, real INTEGRAL, the integral.
#
  import numpy as np
  from r8_gamma import r8_gamma

  if ( e[0] < 0 or e[1] < 0 or e[2] < 0 ):
    print ( '' )
    print ( 'BALL01_MONOMIAL_INTEGRAL - Fatal error!' )
    print ( '  All exponents must be nonnegative.' )
    exit ( 'BALL01_MONOMIAL_INTEGRAL - Fatal error!' )
#
#  Integrate over the surface.
#
  if ( e[0] == 0 and e[1] == 0 and e[2] == 0 ):

    integral = 2.0 * np.sqrt ( np.pi ** 3 ) / r8_gamma ( 1.5 )

  elif ( ( e[0] % 2 ) == 1 or ( e[1] % 2 ) == 1 or ( e[2] % 2 ) == 1 ):

    integral = 0.0

  else:

    integral = 2.0

    for i in range ( 0, 3 ):
      integral = integral * r8_gamma ( 0.5 * float ( e[i] + 1 ) )

    integral = integral / r8_gamma ( 0.5 * float ( e[0] + e[1] + e[2] + 3 ) )
#
#  The surface integral is now adjusted to give the volume integral.
#
  r = 1.0
  s = e[0] + e[1] + e[2] + 3

  integral = integral * r ** s / float ( s )

  return integral

def ball01_monomial_integral_test ( ):

#*****************************************************************************80
#
#% BALL01_MONOMIAL_INTEGRAL_TEST tests BALL01_MONOMIAL_INTEGRAL.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    21 June 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from ball01_sample import ball01_sample
  from ball01_volume import ball01_volume
  from i4vec_uniform_ab import i4vec_uniform_ab
  from monomial_value import monomial_value

  m = 3
  n = 4192
  test_num = 20

  print ( '' )
  print ( 'BALL01_MONOMIAL_INTEGRAL_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  BALL01_MONOMIAL_INTEGRAL returns the integral of a monomial' )
  print ( '  over the interior of the unit ball in 3D.' )
  print ( '  Compate against Monte Carlo estimates.' )
#
#  Get sample points.
#
  seed = 123456789
  x, seed = ball01_sample ( n, seed )

  print ( '' )
  print ( '  Number of sample points used is %d' % ( n ) )
#
#  Randomly choose X,Y exponents between 0 and 8.
#
  print ( '' )
  print ( '  If any exponent is odd, the integral is zero.' )
  print ( '  We will restrict this test to randomly chosen even exponents.' )
  print ( '' )
  print ( '  Ex  Ey  Ez     MC-Estimate           Exact      Error' )
  print ( '' )

  for test in range ( 0, test_num ):

    e, seed = i4vec_uniform_ab ( m, 0, 4, seed )

    for i in range ( 0, 3 ):
      e[i] = e[i] * 2

    value = monomial_value ( m, n, e, x )

    result = ball01_volume ( ) * np.sum ( value ) / float ( n )
    exact = ball01_monomial_integral ( e )
    error = abs ( result - exact )

    print ( '  %2d  %2d  %2d  %14.6g  %14.6g  %10.2e' \
      % ( e[0], e[1], e[1], result, exact, error ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'BALL01_MONOMIAL_INTEGRAL_TEST' )
  print ( '  Normal end of execution.' )
  return

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