#! /usr/bin/env python
#
def p_quadrature_rule ( nt ):

#*****************************************************************************80
#
## P_QUADRATURE_RULE: quadrature for Legendre function P(n,x).
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    16 March 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer NT, the order of the rule.
#
#    Output, real T(NT), W(NT), the points and weights
#    of the rule.
#
  import numpy as np
  from imtqlx import imtqlx

  a = np.zeros ( nt )

  b = np.zeros ( nt )

  for i in range ( 0, nt ):
    ip1 = i + 1
    b[i] = ip1 / np.sqrt ( 4 * ip1 * ip1 - 1 )

  c = np.zeros ( nt )
  c[0] = np.sqrt ( 2.0 )

  t, w = imtqlx ( nt, a, b, c )

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

  return t, w

def p_quadrature_rule_test ( ):

#*****************************************************************************80
#
## P_QUADRATURE_RULE_TEST tests P_QUADRATURE_RULE.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    16 March 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from p_integral import p_integral
  from r8vec2_print import r8vec2_print

  print ( '' )
  print ( 'P_QUADRATURE_RULE_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  P_QUADRATURE_RULE computes the quadrature rule' )
  print ( '  associated with P(n,x)' )

  n = 5
  x, w = p_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 dx' )
  print ( '' )
  print ( '   E       Q_Estimate      Q_Exact' )
  print ( '' )

  for e in range ( 0, 2 * n ):
    if ( e == 0 ):
      f = np.ones ( n )
    else:
      f = x ** e
    q = np.dot ( w, f )
    q_exact = p_integral ( e )
    print ( '  %2d  %14g  %14g' % ( e, q, q_exact ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'P_QUADRATURE_RULE_TEST' )
  print ( '  Normal end of execution.' )
  return

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