#! /usr/bin/env python
#
def pn_pair_product ( p ):

#*****************************************************************************80
#
## PN_PAIR_PRODUCT: pair products for normalized Legendre polynomial Pn(n,x).
#
#  Discussion:
#
#    Let Pn(n,x) represent the normalized Legendre polynomial of degree n.  
#
#    To check orthonormality, we compute
#
#      Tij = Integral ( -1.0 <= X <= +1.0 ) Pn(i,x) * Pn(j,x) dx
#
#    We will estimate these integrals using Gauss-Legendre quadrature.
#
#    The computed table should be the identity matrix.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    17 March 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer P, the maximum degree of the polyonomial 
#    factors.  0 <= P.
#
#    Output, real TABLE(1:P+1,1:P+1), the table of integrals.  
#
  import numpy as np
  from p_quadrature_rule import p_quadrature_rule
  from pn_polynomial_value import pn_polynomial_value

  table = np.zeros ( [ p + 1, p + 1 ] )
  xvec = np.zeros ( 1 )

  order = 2 * p + 1

  x_table, w_table = p_quadrature_rule ( order )

  for k in range ( 0, order ):

    x = x_table[k]
    xvec[0] = x

    h_table = pn_polynomial_value ( 1, p, xvec )

    for i in range ( 0, p + 1 ):
      for j in range ( 0, p + 1 ):
        table[i,j] = table[i,j] + w_table[k] * h_table[0,i] * h_table[0,j]

  return table

def pn_pair_product_test ( p ):

#*****************************************************************************80
#
## PN_PAIR_PRODUCT_TEST tests PN_PAIR_PRODUCT.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    17 March 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer P, the maximum degree of the polynomial 
#    factors.
#
  import platform
  from r8mat_print import r8mat_print

  print ( '' )
  print ( 'PN_PAIR_PRODUCT_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  PN_PAIR_PRODUCT computes a pair product table for Pn(n,x):' )
  print ( '' )
  print ( '  Tij = integral ( -1 <= x <= +1 ) Pn(i,x) Pn(j,x) dx' )
  print ( '' )
  print ( '  The Pn(n,x) polynomials are orthonormal,' )
  print ( '  so T should be the identity matrix.' )
  print ( '' )
  print ( '  Maximum degree P = %d' % ( p ) )

  table = pn_pair_product ( p )

  r8mat_print ( p + 1, p + 1, table, '  Pair product table:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'PN_PAIR_PRODUCT_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( ) 
  p = 5
  pn_pair_product_test ( p )
  timestamp ( )