#! /usr/bin/env python
#
def annulus_sample ( pc, r1, r2, n, seed ):

#*****************************************************************************80
#
## ANNULUS_SAMPLE samples a circular annulus.
#
#  Discussion:
#
#    A circular annulus with center PC, inner radius R1 and
#    outer radius R2, is the set of points P so that
#
#      R1^2 <= (P(1)-PC(1))^2 + (P(2)-PC(2))^2 <= R2^2
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    05 July 2018
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Peter Shirley,
#    Nonuniform Random Point Sets Via Warping,
#    Graphics Gems, Volume III,
#    edited by David Kirk,
#    AP Professional, 1992, 
#    ISBN: 0122861663,
#    LC: T385.G6973.
#
#  Parameters:
#
#    Input, real PC(2), the center.
#
#    Input, real R1, R2, the inner and outer radii.
#    0.0 <= R1 <= R2.
#
#    Input, integer N, the number of points to generate.
#
#    Input/output, integer SEED, a seed for the random number generator.
#
#    Output, real P(2,N), sample points.
#
  import numpy as np
  from r8vec_uniform_01 import r8vec_uniform_01
  from sys import exit

  if ( r1 < 0.0 ):
    print ( '' )
    print ( 'ANNULUS_SAMPLE - Fatal error!' )
    print ( '  Inner radius R1 < 0.0.' )
    exit ( 'ANNULUS_SAMPLE - Fatal error!' )

  if ( r2 < r1 ):
    print ( '' )
    print ( 'ANNULUS_SAMPLE - Fatal error!' )
    print ( '  Outer radius R1 < R1 = inner radius.' )
    exit ( 'ANNULUS_SAMPLE - Fatal error!' )

  u, seed = r8vec_uniform_01 ( n, seed )

  theta = u[:] * 2.0 * np.pi

  v, seed = r8vec_uniform_01 ( n, seed )

  r = np.sqrt ( ( 1.0 - v[:] ) * r1 * r1 \
    +                   v[:]   * r2 * r2 )

  p = np.zeros ( [ 2, n ] )

  p[0,:] = pc[0] + r[:] * np.cos ( theta[:] )
  p[1,:] = pc[1] + r[:] * np.sin ( theta[:] )

  return p, seed

def annulus_sample_test ( center, r1, r2 ):

#*****************************************************************************80
#
## ANNULUS_SAMPLE_TEST uses ANNULUS_SAMPLE to estimate integrals.
#
#  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:
#
#    05 July 2018
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  from annulus_area import annulus_area
  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_SAMPLE_TEST' )
  print ( '  ANNULUS_SAMPLE can sample an annulus uniformly.' )
  print ( '  Use it to estimate integrals in the annulus' )
  print ( '  centered at (%g,%g) with R1 = %g, R2 = %g' \
    % ( center[0], center[1], r1, r2 ) )

  seed = 123456789

  print ( '' )
  print ( '         N        1              X^2             Y^2             X^4             X^2Y^2           Y^4             X^6' )
  print ( '' )

  n = 1

  while ( n <= 65536 ):

    x, seed = annulus_sample ( center, r1, r2, n, seed )

    print ( '  %8d' % ( n ), end = '' )
    for i in range ( 0, 7 ):
      e = e_test[i,:]
      value = monomial_value ( 2, n, e, x )
      result = annulus_area ( center, r1, r2 ) * np.sum ( value[:] ) / n
      print ( '  %14.6g' % ( result ), end = '' )
    print ( '' )

    n = 2 * n

  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 ( '' )

  return

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