#! /usr/bin/env python
#
def p_polynomial_prime2 ( m, n, x ):

#*****************************************************************************80
#
## P_POLYNOMIAL_PRIME2: second derivative of Legendre polynomials P(n,x).
#
#  Discussion:
#
#    P(0,X) = 1
#    P(1,X) = X
#    P(N,X) = ( (2*N-1)*X*P(N-1,X)-(N-1)*P(N-2,X) ) / N
#
#    P'(0,X) = 0
#    P'(1,X) = 1
#    P'(N,X) = ( (2*N-1)*(P(N-1,X)+X*P'(N-1,X)-(N-1)*P'(N-2,X) ) / N
#
#    P"(0,X) = 0
#    P"(1,X) = 0
#    P"(N,X) = ( (2*N-1)*(2*P'(N-1,X)+X*P"(N-1,X)-(N-1)*P'(N-2,X) ) / N
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    16 March 2016
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Milton Abramowitz, Irene Stegun,
#    Handbook of Mathematical Functions,
#    National Bureau of Standards, 1964,
#    ISBN: 0-486-61272-4,
#    LC: QA47.A34.
#
#    Daniel Zwillinger, editor,
#    CRC Standard Mathematical Tables and Formulae,
#    30th Edition,
#    CRC Press, 1996.
#
#  Parameters:
#
#    Input, integer M, the number of evaluation points.
#
#    Input, integer N, the highest order polynomial to evaluate.
#    Note that polynomials 0 through N will be evaluated.
#
#    Input, real X(M,1), the evaluation points.
#
#    Output, real VPP(M,N+1), the second derivatives of the
#    Legendre polynomials of order 0 through N at the point X.
#
  import numpy as np

  v = np.zeros ( [ m, n + 1 ] )
  vp = np.zeros ( [ m, n + 1 ] )
  vpp = np.zeros ( [ m, n + 1 ] )

  for i in range ( 0, m ):
    v[i,0] = 1.0
    vp[i,0] = 0.0
    vpp[i,0] = 0.0

  if ( n < 1 ):
    return vpp

  for i in range ( 0, m ):
    v[i,1] = x[i]
    vp[i,1] = 1.0
    vpp[i,1] = 0.0
 
  for j in range ( 2, n + 1 ):
 
    for i in range ( 0, m ):

      v[i,j] = ( ( 2 * j - 1 ) * x[i] * v[i,j-1]   \
               - (     j - 1 ) *        v[i,j-2] ) \
               / (     j     )
 
      vp[i,j] = ( ( 2 * j - 1 ) * ( v[i,j-1] + x[i] * vp[i,j-1] ) \
                - (     j - 1 ) *                     vp[i,j-2] ) \
                / (     j     )

      vpp[i,j] = ( ( 2 * j - 1 ) * ( 2.0 * vp[i,j-1] + x[i] * vpp[i,j-1] ) \
                 - (     j - 1 ) *                            vpp[i,j-2] ) \
                 / (     j     )
  
  return vpp

def p_polynomial_prime2_test ( ):

#*****************************************************************************80
#
## P_POLYNOMIAL_PRIME2_TEST tests P_POLYNOMIAL_PRIME2.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    16 March 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'P_POLYNOMIAL_PRIME2_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  P_POLYNOMIAL_PRIME2 evaluates the second derivative of' )
  print ( '  the Legendre polynomial P(N,X).' )
  print ( '' )
  print ( '                                      Computed' )
  print ( '     N        X                         P"(N,X)' )

  m = 11
  n = 5
  x = np.linspace ( -1.0, +1.0, m )
  vpp = p_polynomial_prime2 ( m, n, x )

  for i in range ( 0, m ):
    print ( '' )
    for j in range ( 0, n + 1 ):
      print ( '  %4d  %12g  %24g' % ( j, x[i], vpp[i,j] ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'P_POLYNOMIAL_PRIME2_TEST' )
  print ( '  Normal end of execution.' )
  return

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