#! /usr/bin/env python
#
def annulus_rule_monomial_test ( center, r1, r2 ):

#*****************************************************************************80
#
## ANNULUS_RULE_MONOMIAL_TEST estimates monomial integrals using quadrature.
#
#  Discussion:
#
#    If CENTER=(0,0) and R1 = 0 and R2 = 1, then we can compare exact values.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 July 2018
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, real CENTER(2), the coordinates of the center.
#
#    Input, real R1, R2, the inner and outer radii of the annulus.
#    0.0 <= R1 <= R2.
#
  import numpy as np
  import platform
  from annulus_rule_compute import annulus_rule_compute
  from disk01_monomial_integral import disk01_monomial_integral
  from monomial_value import monomial_value

  e_test = np.array ( [ \
    [ 0, 0 ], \
    [ 2, 0 ], \
    [ 0, 2 ], \
    [ 4, 0 ], \
    [ 2, 2 ], \
    [ 0, 4 ], \
    [ 6, 0 ] ] )

  print ( '' )
  print ( 'ANNULUS_RULE_MONOMIAL_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  ANNULUS_RULE_COMPUTE can supply a quadrature rule for' )
  print ( '  the annulus centered at (%g,%g) with R1 = %g, R2 = %g' \
    % ( center[0], center[1], r1, r2 ) )
  print ( '  Apply this rule to a variety of monomials.' )

  print ( '' )
  print ( '    NR    NT           1              X^2             Y^2             X^4               X^2Y^2           Y^4             X^6' )
  print ( '' )

  nr = 4
  
  while ( nr <= 64 ):

    nt = 4 * nr
    n = nr * nt

    w, x, y = annulus_rule_compute ( center, r1, r2, nr, nt )

    xy = np.zeros ( [ 2, n ] )
    xy[0,:] = x[:]
    xy[1,:] = y[:]

    print ( '  %4d  %4d' % ( nr, nt ), end = '' )
 
    for i in range ( 0, 7 ):
      e = e_test[i,:]
      value = monomial_value ( 2, n, e, xy )
      result = np.sum ( w[:] * value[:] )
      print ( '  %14.6g' % ( result ), end = '' )
    print ( '' )

    nr = 2 * nr

  if ( \
    center[0] == 0.0 and \
    center[1] == 0.0 and \
    r1 == 0.0 and \
    r2 == 1.0 ):
    print ( '' )
    print ( '     Exact  ', end = '' )
    for i in range ( 0, 7 ):
      e = e_test[i,:]
      result = disk01_monomial_integral ( e )
      print ( '  %14.6g' % ( result ), end = '' )
    print ( '' )
#
#  Terminate.
#
  print ( '' )
  print ( 'ANNULUS_RULE_MONOMIAL_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  import numpy as np
  timestamp ( )
  center = np.array ( [ 0.0, 0.0 ] )
  r1 = 0.0
  r2 = 1.0
  annulus_rule_monomial_test ( center, r1, r2 )
  timestamp ( )