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

#*****************************************************************************80
#
## PMNS_POLYNOMIAL_VALUE: sphere_normalized Legendre polynomial Pmn(n,m,x).
#
#  Discussion:
#
#    The unnormalized associated Legendre functions P_N^M(X) have
#    the property that
#
#      Integral ( -1 <= X <= 1 ) ( P_N^M(X) )^2 dX 
#      = 2 * ( N + M )! / ( ( 2 * N + 1 ) * ( N - M )! )
#
#    By dividing the function by the square root of this term,
#    the normalized associated Legendre functions have norm 1.
#
#    However, we plan to use these functions to build spherical
#    harmonics, so we use a slightly different normalization factor of
#
#      sqrt ( ( ( 2 * N + 1 ) * ( N - M )! ) / ( 4 * pi * ( N + M )! ) ) 
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    17 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.
#
#  Parameters:
#
#    Input, integer MM, the number of evaluation points.
#
#    Input, integer N, the maximum first index of the Legendre
#    function, which must be at least 0.
#
#    Input, integer M, the second index of the Legendre function,
#    which must be at least 0, and no greater than N.
#
#    Input, real X(MM,1), the evaluation points.
#
#    Output, real CX(MM,N+1), the function values.
#
  import numpy as np
  from r8_factorial import r8_factorial

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

  if ( m <= n ):
    for i in range ( 0, mm ):
      cx[i,m] = 1.0 
    factor = 1.0
    for j in range ( 0, m ):
      for i in range ( 0, mm ):
        cx[i,m] = - factor * cx[i,m] * np.sqrt ( 1.0 - x[i] ** 2 )
      factor = factor + 2.0

  if ( m + 1 <= n ):
    for i in range ( 0, mm ):
      cx[i,m+1] = ( 2 * m + 1 ) * x[i] * cx[i,m]

  for j in range ( m + 2, n + 1 ):
    for i in range ( 0, mm ):
      cx[i,j] = ( ( 2 * j     - 1 ) * x[i] * cx[i,j-1] \
                + (   - j - m + 1 ) *        cx[i,j-2] ) \
                / (     j - m     )
#
#  Normalization.
#
  for j in range ( m, n + 1 ):
    factor = np.sqrt ( ( ( 2 * j + 1 ) * r8_factorial ( j - m ) ) \
      / ( 4.0 * np.pi * r8_factorial ( j + m ) ) )
    for i in range ( 0, mm ):
      cx[i,j] = cx[i,j] * factor

  return cx

def pmns_polynomial_value_test ( ):

#*****************************************************************************80
#
## PMNS_POLYNOMIAL_VALUE_TEST tests PMNS_POLYNOMIAL_VALUE.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    17 March 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from pmns_polynomial_values import pmns_polynomial_values

  mm = 1
  xvec = np.zeros ( 1 )

  print ( '' )
  print ( 'PMNS_POLYNOMIAL_VALUE_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  PMNS_POLYNOMIAL_VALUE evaluates the Legendre polynomial Pmns(n,m,x).' )
  print ( '' )
  print ( '                                         Tabulated                 Computed' )
  print ( '     N     M             X               Pmns(N,M,X)               Pmns(N,M,X)         Error' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, n, m, x, fx1 = pmns_polynomial_values ( n_data )

    if ( n_data == 0 ):
      break

    xvec[0] = x
    v = pmns_polynomial_value ( mm, n, m, xvec )
    fx2 = v[0,n]

    e = fx1 - fx2

    print ( '  %4d  %4d  %12g  %24.16g  %24.16g  %8g' % ( n, m, x, fx1, fx2, e ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'PMNS_POLYNOMIAL_VALUE_TEST' )
  print ( '  Normal end of execution.' )
  return

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