#! /usr/bin/env python
#
def u_quadrature_rule ( n ):

#*****************************************************************************80
#
## U_QUADRATURE_RULE: quadrature rule for U(n,x).
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    15 July 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of the rule.
#
#    Output, real T(N,1), W(N,1), the points and weights of the rule.
#
  import numpy as np
  from imtqlx import imtqlx

  aj = np.zeros ( n )

  bj = np.ones ( n )
  for i in range ( 0, n ):
    bj[i] = 0.5 * bj[i]

  cj = np.zeros ( n )
  cj[0] = np.sqrt ( np.pi / 2.0 )

  t, w = imtqlx ( n, aj, bj, cj )

  for i in range ( 0, n ):
    w[i] = w[i] ** 2

  return t, w

def u_quadrature_rule_test ( ):

#*****************************************************************************80
#
## U_QUADRATURE_RULE_TEST tests U_QUADRATURE_RULE.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    15 July 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from r8vec2_print import r8vec2_print
  from u_moment import u_moment

  print ( '' )
  print ( 'U_QUADRATURE_RULE_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  U_QUADRATURE_RULE computes the quadrature rule' )
  print ( '  associated with U(n,x)' )

  n = 7
  x, w = u_quadrature_rule ( n )

  r8vec2_print ( n, x, w, '      X            W' )

  print ( '' )
  print ( '  Use the quadrature rule to estimate:' )
  print ( '' )
  print ( '    Q = Integral ( -1 <= X <= +1 ) X^E * sqrt ( 1-x^2) dx' )
  print ( '' )
  print ( '   E       Q_Estimate      Q_Exact' )
  print ( '' )

  f = np.zeros ( n )

  for e in range ( 0, 2 * n ):
    if ( e == 0 ):
      for i in range ( 0, n ):
        f[i] = 1.0
    else:
      for i in range ( 0, n ):
        f[i] = x[i] ** e

    q = np.dot ( w, f )
    q_exact = u_moment ( e )
    print ( '  %2d  %14g  %14g' % ( e, q, q_exact ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'U_QUADRATURE_RULE_TEST:' )
  print ( '  Normal end of execution.' )
  return

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