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

#*****************************************************************************80
#
## V_MASS_MATRIX computes the mass matrix for the Chebyshev V polynomial.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    20 July 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  from v_polynomial import v_polynomial
  from v_quadrature_rule import v_quadrature_rule

  x, w = v_quadrature_rule ( n + 1 )

  phi = v_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 v_mass_matrix_test ( ):

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

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

  n = 3

  a = v_mass_matrix ( n )

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

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