#! /usr/bin/env python3
#
def wedge01_monomial_integral ( e ):

#*****************************************************************************80
#
## WEDGE01_MONOMIAL_INTEGRAL: integral of a monomial in the unit wedge in 3D.
#
#  Discussion:
#
#    This routine returns the integral of
#
#      product ( 1 <= I <= 3 ) X(I)^E(I)
#
#    over the unit wedge.
#
#    The integration region is:
#
#      0 <= X
#      0 <= Y
#      X + Y <= 1
#      -1 <= Z <= 1.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 June 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Arthur Stroud,
#    Approximate Calculation of Multiple Integrals,
#    Prentice Hall, 1971,
#    ISBN: 0130438936,
#    LC: QA311.S85.
#
#  Parameters:
#
#    Input, integer E(3), the exponents.
#
#    Output, real VALUE, the integral of the monomial.
#
  from sys import exit

  value = 1.0

  k = e[0]

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

  k = k + 1
  value = value / float ( k )

  k = k + 1
  value = value / float ( k )
#
#  Now account for integration in Z.
#
  if ( e[2] == - 1 ):
    print ( '' )
    print ( 'WEDGE01_MONOMIAL_INTEGRAL - Fatal error!' )
    print ( '  E(3) = -1 is not a legal input.' )
    exit ( 'WEDGE01_MONOMIAL_INTEGRAL - Fatal error!' )
  elif ( ( e[2] % 2 ) == 1 ):
    value = 0.0
  else:
    value = value * 2.0 / float ( e[2] + 1 )

  return value

def wedge01_monomial_integral_test ( ):

#*****************************************************************************80
#
## WEDGE01_MONOMIAL_INTEGRAL_TEST tests WEDGE01_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 wedge01_sample import wedge01_sample
  from wedge01_volume import wedge01_volume

  m = 3
  n = 500000
  e_max = 6

  print ( '' )
  print ( 'WEDGE01_MONOMIAL_INTEGRAL_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  WEDGE01_MONOMIAL_INTEGRAL computes the integral of a monomial' )
  print ( '  over the interior of the unit wedge in 3D.' )
  print ( '  Compare with a Monte Carlo estimate.' )

  seed = 123456789
  x, seed = wedge01_sample ( n, seed )

  print ( '' )
  print ( '  Number of sample points used is %d' % ( n ) )
  print ( '' )
  print ( '   E1  E2  E3     MC-Estimate      Exact           Error' )
  print ( '' )
#
#  Check all monomials up to total degree E_MAX.
#
  e = np.zeros ( 3, dtype = np.int32 )

  for e3 in range ( 0, e_max + 1 ):
    e[2] = e3
    for e2 in range ( 1, e_max - e3 + 1 ):
      e[1] = e2
      for e1 in range ( 0, e_max - e3 - e2 + 1 ):
        e[0] = e1

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

        q = wedge01_volume ( ) * np.sum ( value ) / float ( n )
        exact = wedge01_monomial_integral ( e )
        error = abs ( q - exact )

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

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