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

#*****************************************************************************80
#
## U_MASS_MATRIX computes the mass matrix for the Chebyshev T polynomial.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    12 July 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  from u_polynomial import u_polynomial
  from u_quadrature_rule import u_quadrature_rule

  x, w = u_quadrature_rule ( n + 1 )

  phi = u_polynomial ( n + 1, n, x )

  phiw = np.zeros ( [ n + 1, n + 1 ] )

  for i in range ( 0, n + 1 ):
    for j in range ( 0, n + 1 ):
      phiw[j,i] = w[i] * phi[i,j]

  a = np.dot ( phiw, phi )

  return a

def u_mass_matrix_test ( ):

#*****************************************************************************80
#
## U_MASS_MATRIX_TEST tests U_MASS_MATRIX.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    15 July 2015
#
#  Author:
#
#    John Burkardt
#
  import platform
  from r8mat_print import r8mat_print

  print ( '' )
  print ( 'U_MASS_MATRIX_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  U_MASS_MATRIX computes the mass matrix for the' )
  print ( '  Chebyshev U polynomials U(i,x).' )
  print ( '  A(I,J) = integral ( -1 <=x <= +1 ) U(i,x) U(j,x) * sqrt ( 1 - x^2 ) dx' )
  print ( '  0    if i is not equal to j' )
  print ( '  pi/2 if i = j.' )

  n = 3

  a = u_mass_matrix ( n )

  r8mat_print ( n + 1, n + 1, a, '  U mass matrix:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'U_MASS_MATRIX_TEST:' )
  print ( '  Normal end of execution.' )
  return

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