#! /usr/bin/env python
#
def annulus_rule_compute ( center, r1, r2, nr, nt ):

#*****************************************************************************80
#
## ANNULUS_RULE_COMPUTE computes a quadrature rule for an annulus.
#
#  Discussion:
#
#    The integration region is points (X,Y) such that
#
#      R1^2 <= ( X - CENTER(1) )^2 + ( Y - CENTER(2) )^2 <= R2^2
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    08 July 2018
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    William Peirce,
#    Numerical Integration Over the Planar Annulus,
#    Journal of the Society for Industrial and Applied Mathematics,
#    Volume 5, Issue 2, June 1957, pages 66-73.
#
#  Parameters:
#
#    Input, real CENTER(2), the center of the annulus.
#
#    Input, real R1, R2, the inner and outer radius.
#
#    Input, integer NR, the number of points in the radial rule.
#
#    Input, integer NT, the number of angles to use.
#    The value NT=4*NR is recommended.
#
#    Output, real W(NR*NT), the weights for the rule.
#
#    Output, real X(NR*NT), Y(NR*NT), the points for the rule.
#
  import numpy as np
  from annulus_area import annulus_area
  from legendre_ek_compute import legendre_ek_compute
  from rule_adjust import rule_adjust
#
#  Get the Legendre rule for [-1,+1].
#
  ra, rw = legendre_ek_compute ( nr )
#
#  Adjust the rule from [-1,+1] to [r1^2,r2^2].
#
  a = -1.0
  b = +1.0
  c = r1 * r1
  d = r2 * r2
  ra, rw = rule_adjust ( a, b, c, d, nr, ra, rw )
#
#  Convert from R^2 to R.
#
  ra = np.sqrt ( ra )
  rw = rw / ( r2 + r1 ) / ( r2 - r1 )
#
#  Set the angular weight.
#
  tw = 1.0 / nt
#
#  Get area of annulus.
#
  area = annulus_area ( center, r1, r2 )
#
#  Form the abscissa and weight vectors.
#
  x = np.zeros ( nr * nt )
  y = np.zeros ( nr * nt )
  w = np.zeros ( nr * nt )

  k = 0
  for i in range ( 0, nt ):
    t = 2.0 * np.pi * float ( i ) / float ( nt )
    for j in range ( 0, nr ):
      x[k] = center[0] + ra[j] * np.cos ( t )
      y[k] = center[1] + ra[j] * np.sin ( t )
      w[k] = area * tw * rw[j]
      k = k + 1

  return w, x, y

def annulus_rule_compute_test ( ):

#*****************************************************************************80
#
## ANNULUS_RULE_COMPUTE_TEST tests ANNULUS_RULE_COMPUTE.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 July 2018
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from r8vec3_print import r8vec3_print

  print ( '' )
  print ( 'ANNULUS_RULE_COMPUTE_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Test ANNULUS_RULE_COMPUTE.' )

  center = np.array ( [ 0.0, 0.0 ] )
  r1 = 0.5
  r2 = 1.0
  nr = 3
  nt = 12

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

  r8vec3_print ( nr * nt, w, x, y, '  W, X, Y for annulus quadrature:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'ANNULUS_RULE_COMPUTE_TEST:' )
  print ( '  Normal end of execution.' )
  return

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