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

#*****************************************************************************80
#
## HERMITE_GK24_SET sets a Genz-Keister 24 Hermite rule.
#
#  Discussion:
#
#    The integral:
#
#      integral ( -oo <= x <= +oo ) f(x) exp ( - x * x ) dx
#
#    The quadrature rule:
#
#      sum ( 1 <= i <= n ) w(i) * f ( x(i) )
#
#    A nested family of rules for the Hermite integration problem
#    was produced by Genz and Keister.  The structure of the nested
#    family was denoted by 1+2+6+10+24, that is, it comprised rules
#    of successive orders O = 1, 3, 9, 19, and 43.
#
#    The precisions of these rules are P = 1, 5, 15, 29, and 67.
#
#    Some of the data in this function was kindly supplied directly by
#    Alan Genz on 24 April 2011.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    20 June 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Alan Genz, Bradley Keister,
#    Fully symmetric interpolatory rules for multiple integrals
#    over infinite regions with Gaussian weight,
#    Journal of Computational and Applied Mathematics,
#    Volume 71, 1996, pages 299-309
#
#    Thomas Patterson,
#    The Optimal Addition of Points to Quadrature Formulae,
#    Mathematics of Computation,
#    Volume 22, Number 104, October 1968, pages 847-856.
#
#  Parameters:
#
#    Input, integer N, the order.
#    Legal values are 1, 3, 9, 19, and 43.
#
#    Output, real X(N,1), the abscissas.
#
#    Output, real W(N,1), the weights.
#
  import numpy as np
  from sys import exit

  if ( n == 1 ):

    x = np.array ( [ \
              0.0000000000000000E+00 ] )

    w = np.array ( [ \
              1.7724538509055159E+00 ] )

  elif ( n == 3 ):

    x = np.array ( [ \
             -1.2247448713915889E+00, \
              0.0000000000000000E+00, \
              1.2247448713915889E+00 ] )

    w = np.array ( [ \
              2.9540897515091930E-01, \
              1.1816359006036772E+00, \
              2.9540897515091930E-01 ] )

  elif ( n == 9 ):

    x = np.array ( [ \
             -2.9592107790638380E+00, \
             -2.0232301911005157E+00, \
             -1.2247448713915889E+00, \
             -5.2403354748695763E-01, \
              0.0000000000000000E+00, \
              5.2403354748695763E-01, \
              1.2247448713915889E+00, \
              2.0232301911005157E+00, \
              2.9592107790638380E+00 ] )

    w = np.array ( [ \
              1.6708826306882348E-04, \
              1.4173117873979098E-02, \
              1.6811892894767771E-01, \
              4.7869428549114124E-01, \
              4.5014700975378197E-01, \
              4.7869428549114124E-01, \
              1.6811892894767771E-01, \
              1.4173117873979098E-02, \
              1.6708826306882348E-04 ] )

  elif ( n == 19 ):

    x = np.array ( [ \
             -4.4995993983103881E+00, \
             -3.6677742159463378E+00, \
             -2.9592107790638380E+00, \
             -2.2665132620567876E+00, \
             -2.0232301911005157E+00, \
             -1.8357079751751868E+00, \
             -1.2247448713915889E+00, \
             -8.7004089535290285E-01, \
             -5.2403354748695763E-01, \
              0.0000000000000000E+00, \
              5.2403354748695763E-01, \
              8.7004089535290285E-01, \
              1.2247448713915889E+00, \
              1.8357079751751868E+00, \
              2.0232301911005157E+00, \
              2.2665132620567876E+00, \
              2.9592107790638380E+00, \
              3.6677742159463378E+00, \
              4.4995993983103881E+00 ] )

    w = np.array ( [ \
              1.5295717705322357E-09, \
              1.0802767206624762E-06, \
              1.0656589772852267E-04, \
              5.1133174390883855E-03, \
             -1.1232438489069229E-02, \
              3.2055243099445879E-02, \
              1.1360729895748269E-01, \
              1.0838861955003017E-01, \
              3.6924643368920851E-01, \
              5.3788160700510168E-01, \
              3.6924643368920851E-01, \
              1.0838861955003017E-01, \
              1.1360729895748269E-01, \
              3.2055243099445879E-02, \
             -1.1232438489069229E-02, \
              5.1133174390883855E-03, \
              1.0656589772852267E-04, \
              1.0802767206624762E-06, \
              1.5295717705322357E-09 ] )

  elif ( n == 43 ):

    x = np.array ( [ \
            -10.167574994881873, \
             -7.231746029072501, \
             -6.535398426382995, \
             -5.954781975039809, \
             -5.434053000365068, \
             -4.952329763008589, \
             -4.4995993983103881, \
             -4.071335874253583, \
             -3.6677742159463378, \
             -3.295265921534226, \
             -2.9592107790638380, \
             -2.633356763661946, \
             -2.2665132620567876, \
             -2.089340389294661, \
             -2.0232301911005157, \
             -1.8357079751751868, \
             -1.583643465293944, \
             -1.2247448713915889, \
             -0.87004089535290285, \
             -0.52403354748695763, \
             -0.196029453662011, \
              0.0000000000000000, \
              0.196029453662011, \
              0.52403354748695763, \
              0.87004089535290285, \
              1.2247448713915889, \
              1.583643465293944, \
              1.8357079751751868, \
              2.0232301911005157, \
              2.089340389294661, \
              2.2665132620567876, \
              2.633356763661946, \
              2.9592107790638380, \
              3.295265921534226, \
              3.6677742159463378, \
              4.071335874253583, \
              4.4995993983103881, \
              4.952329763008589, \
              5.434053000365068, \
              5.954781975039809, \
              6.535398426382995, \
              7.231746029072501, \
             10.167574994881873 ] )

    w = np.array ( [ \
              0.546191947478318097E-37, \
              0.87544909871323873E-23, \
              0.992619971560149097E-19, \
              0.122619614947864357E-15, \
              0.421921851448196032E-13, \
              0.586915885251734856E-11, \
              0.400030575425776948E-09, \
              0.148653643571796457E-07, \
              0.316018363221289247E-06, \
              0.383880761947398577E-05, \
              0.286802318064777813E-04, \
              0.184789465688357423E-03, \
              0.150909333211638847E-02, \
            - 0.38799558623877157E-02, \
              0.67354758901013295E-02, \
              0.139966252291568061E-02, \
              0.163616873493832402E-01, \
              0.450612329041864976E-01, \
              0.928711584442575456E-01, \
              0.145863292632147353E+00, \
              0.164880913687436689E+00, \
              0.579595986101181095E-01, \
              0.164880913687436689E+00, \
              0.145863292632147353E+00, \
              0.928711584442575456E-01, \
              0.450612329041864976E-01, \
              0.163616873493832402E-01, \
              0.139966252291568061E-02, \
              0.67354758901013295E-02, \
            - 0.38799558623877157E-02, \
              0.150909333211638847E-02, \
              0.184789465688357423E-03, \
              0.286802318064777813E-04, \
              0.383880761947398577E-05, \
              0.316018363221289247E-06, \
              0.148653643571796457E-07, \
              0.400030575425776948E-09, \
              0.586915885251734856E-11, \
              0.421921851448196032E-13, \
              0.122619614947864357E-15, \
              0.992619971560149097E-19, \
              0.87544909871323873E-23, \
              0.546191947478318097E-37 ] )

  else:

    print ( '' )
    print ( 'HERMITE_GK24_SET - Fatal error!' )
    print ( '  Illegal input value of N.' )
    print ( '  N must be 1, 3, 9, 19, or 43.' )
    exit ( 'HERMITE_GK24_SET - Fatal error!' )

  return x, w

def hermite_gk24_set_test ( ):

#*****************************************************************************80
#
## HERMITE_GK24_SET_TEST tests HERMITE_GK24_SET.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    20 June 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  l_max = 4

  n_list = np.array ( [ 1, 3, 9, 19, 43 ] )

  print ( '' )
  print ( 'HERMITE_GK24_SET_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  HERMITE_GK24_SET sets up a nested rule' )
  print ( '  for the Hermite integration problem.' )
  print ( '' )
  print ( '  Index       X             W' )

  for l in range ( 0, l_max + 1 ):

    n = n_list[l]

    x, w = hermite_gk24_set ( n )

    print ( '' )

    for i in range ( 0, n ):
      print ( '  %2d  %24.16g  %24.16g' % ( i, x[i], w[i] ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'HERMITE_GK24_SET_TEST:' )
  print ( '  Normal end of execution.' )
  return

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