#! /usr/bin/env python
#
def p_power_product ( p, e ):

#*****************************************************************************80
#
#% P_POWER_PRODUCT: power products for Legendre polynomial P(n,x).
#
#  Discussion:
#
#    Let P(n,x) represent the Legendre polynomial of degree i.  
#
#    For polynomial chaos applications, it is of interest to know the
#    value of the integrals of products of powers of X with every possible pair
#    of basis functions.  That is, we'd like to form
#
#      Tij = Integral ( -1 <= X <= +1 ) X^E * P(i,x) * P(j,x) dx
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    16 March 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer P, the maximum degree of the polyonomial factors.
#    0 <= P.
#
#    Input, integer E, the exponent of X in the integrand.
#    0 <= E.
#
#    Output, real TABLE(P+1,P+1), the table of integrals.
#
  import numpy as np
  from p_polynomial_value import p_polynomial_value
  from p_quadrature_rule import p_quadrature_rule

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

  order = p + 1 + ( ( e + 1 ) // 2 )

  x_table, w_table = p_quadrature_rule ( order )

  for k in range ( 0, order ):

    x = x_table[k]
    xvec[0] = x
    l_table = p_polynomial_value ( 1, p, xvec )
#
#  The following formula is an outer product in L_TABLE.
#
    if ( e == 0 ):
      for i in range ( 0, p + 1 ):
        for j in range ( 0, p + 1 ):
          table[i,j] = table[i,j] + w_table[k] * l_table[0,i] * l_table[0,j]
    else:
      for i in range ( 0, p + 1 ):
        for j in range ( 0, p + 1 ):
          table[i,j] = table[i,j] + w_table[k] * x ** e * l_table[0,i] * l_table[0,j]

  return table

def p_power_product_test ( p, e ):

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

  print ( '' )
  print ( 'P_POWER_PRODUCT_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  P_POWER_PRODUCT computes a power product table for P(n,x):' )
  print ( '' )
  print ( '  Tij = integral ( -1 <= x <= +1 ) x^e P(i,x) P(j,x) dx' )

  print ( '' )
  print ( '  Maximum degree P = %d' % ( p ) )
  print ( '  Exponent of X, E = %d' % ( e ) )

  table = p_power_product ( p, e )

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

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( ) 
  p = 5
  e = 0
  p_power_product_test ( p, e )
  p = 5
  e = 1
  p_power_product_test ( p, e )
  timestamp ( )