#! /usr/bin/env python
#
def p_exponential_product ( p, b ):

#*****************************************************************************80
#
## P_EXPONENTIAL_PRODUCT: exponential products for 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 exp(B*X) with every possible pair
#    of basis functions.  That is, we'd like to form
#
#      Tij = Integral ( -1 <= X <= +1 ) exp(B*x) * P(i,x) * P(j,x) dx
#
#    Because of the exponential factor, the quadrature will not be exact.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    15 March 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer P, the maximum degree of the polyonomial factors.
#    0 <= P.
#
#    Input, real B, the coefficient of X in the exponential factor.
#
#    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

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

  order = ( ( 3 * p + 4 ) // 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.
#
    for i in range ( 0, p + 1 ):
      for j in range ( 0, p + 1 ):
        table[i,j] = table[i,j] \
          + w_table[k] * np.exp ( b * x ) * l_table[0,i] * l_table[0,j]

  return table

def p_exponential_product_test ( p, b ):

#*****************************************************************************80
#
## P_EXPONENTIAL_PRODUCT_TEST tests P_EXPONENTIAL_PRODUCT.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    15 March 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer P, the maximum degree of the polynomial 
#    factors.
#
#    Input, real B, the coefficient of X in the exponential factor.
#
  import platform
  from r8mat_print import r8mat_print

  print ( '' )
  print ( 'P_EXPONENTIAL_PRODUCT_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  P_EXPONENTIAL_PRODUCT computes an exponential product table for P(n,x):' )
  print ( '' )
  print ( '  Tij = integral ( -1 <= x <= +1 ) exp(b*x) P(i,x) P(j,x) dx' )

  print ( '' )
  print ( '  Maximum degree P = %d' % ( p ) )
  print ( '  Exponential argument coefficient B = %g' % ( b ) )

  table = p_exponential_product ( p, b )

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

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  p = 5
  b = 0.0
  p_exponential_product_test ( p, b )
  p = 5
  b = 1.0
  p_exponential_product_test ( p, b )
  timestamp ( )