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

#*****************************************************************************80
#
## TETRAHEDRON01_MONOMIAL_INTEGRAL: integrals in the unit tetrahedron in 3D.
#
#  Discussion:
#
#    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:
#
#    23 June 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer E(3), the exponents.  
#    Each exponent must be nonnegative.
#
#    Output, real INTEGRAL, the integral.
#
  from sys import exit

  m = 3

  for i in range ( 0, m ):
    if ( e[i] < 0 ):
      print ( '' )
      print ( 'TETRAHEDRON01_MONOMIAL_INTEGRAL - Fatal error!' )
      print ( '  All exponents must be nonnegative.' )
      error ( 'TETRAHEDRON01_MONOMIAL_INTEGRAL - Fatal error!\n' )

  k = 0
  integral = 1.0

  for i in range ( 0, m ):

    for j in range ( 1, e[i] + 1 ):
      k = k + 1
      integral = integral * float ( j ) / float ( k )

  for i in range ( 0, m ):
    k = k + 1
    integral = integral / float ( k )

  return integral

def tetrahedron01_monomial_integral_test ( ):

#*****************************************************************************80
#
## TETRAHEDRON_MONOMIAL_INTEGRAL_TEST01 tests TETRAHEDRON_MONOMIAL_INTEGRAL.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 June 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from monomial_value import monomial_value
  from tetrahedron01_sample import tetrahedron01_sample
  from tetrahedron01_volume import tetrahedron01_volume

  m = 3
  n = 4192
  test_num = 20

  print ( '' )
  print ( 'TETRAHEDRON_MONOMIAL_INTEGRAL_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  TETRAHEDRON_MONOMIAL_INTEGRAL returns the integral of a monomial' )
  print ( '  over the interior of the unit tetrahedron in 3D.' )
  print ( '  Compare with a Monte Carlo estimate.' )
#
#  Get sample points.
#
  seed = 123456789
  x, seed = tetrahedron01_sample ( n, seed )

  print ( '' )
  print ( '  Number of sample points used is %d' % ( n ) )
#
#  Run through the exponents.
#
  print ( '' )
  print ( '  Ex  Ey  Ez     MC-Estimate      Exact           Error' )
  print ( '' )

  e = np.zeros ( m, dtype = np.int32 )

  for i in range ( 0, m + 1 ):
    e[0] = i
    for j in range ( 0, m + 1 ):
      e[1] = j
      for k in range ( 0, m + 1 ):
        e[2] = k

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

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

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

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