#! /usr/bin/env python
#
def bessel_i0 ( arg ):

#*****************************************************************************80
#
## BESSEL_I0 evaluates the modified Bessel function of the first kind and order zero.
#
#  Discussion:
#
#    The main computation evaluates slightly modified forms of
#    minimax approximations generated by Blair and Edwards, Chalk
#    River (Atomic Energy of Canada Limited) Report AECL-4928,
#    October, 1974.  This transportable program is patterned after
#    the machine dependent FUNPACK packet NATSI0, but cannot match
#    that version for efficiency or accuracy.  This version uses
#    rational functions that theoretically approximate I-SUB-0(X)
#    to at least 18 significant decimal digits.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    02 March 2016
#
#  Author:
#
#    Original FORTRAN77 version by William Cody, Laura Stoltz.
#    Python version by John Burkardt.
#
#  Parameters:
#
#    Input, real ARG, the argument.
#
#    Output, real VALUE, the value of the modified Bessel function
#    of the first kind.
#
  import numpy as np
  from r8_huge import r8_huge

  exp40 = 2.353852668370199854E+17

  p = np.array ( [ \
    -5.2487866627945699800E-18, \
    -1.5982226675653184646E-14, \
    -2.6843448573468483278E-11, \
    -3.0517226450451067446E-08, \
    -2.5172644670688975051E-05, \
    -1.5453977791786851041E-02, \
    -7.0935347449210549190E+00, \
    -2.4125195876041896775E+03, \
    -5.9545626019847898221E+05, \
    -1.0313066708737980747E+08, \
    -1.1912746104985237192E+10, \
    -8.4925101247114157499E+11, \
    -3.2940087627407749166E+13, \
    -5.5050369673018427753E+14, \
    -2.2335582639474375249E+15 ] )

  pp = np.array ( [ \
    -3.9843750000000000000E-01, \
     2.9205384596336793945E+00, \
    -2.4708469169133954315E+00, \
     4.7914889422856814203E-01, \
    -3.7384991926068969150E-03, \
    -2.6801520353328635310E-03, \
     9.9168777670983678974E-05, \
    -2.1877128189032726730E-06 ] )

  q = np.array ( [ \
    -3.7277560179962773046E+03, \
     6.5158506418655165707E+06, \
    -6.5626560740833869295E+09, \
     3.7604188704092954661E+12, \
    -9.7087946179594019126E+14 ] )

  qq = np.array ( [ \
    -3.1446690275135491500E+01, \
     8.5539563258012929600E+01, \
    -6.0228002066743340583E+01, \
     1.3982595353892851542E+01, \
    -1.1151759188741312645E+00, \
     3.2547697594819615062E-02, \
    -5.5194330231005480228E-04 ] )

  rec15 = 6.6666666666666666666E-02
  xmax = 91.9E+00
  xsmall = 2.98E-08

  x = abs ( arg )

  if ( x < xsmall ):
    value = 1.0
  elif ( x < 15.0 ):
#
#  XSMALL <= ABS(ARG) < 15.0
#
    xx = x * x
    sump = p[0]
    for i in range ( 1, 15 ):
      sump = sump * xx + p[i]

    xx = xx - 225.0
    sumq = (((( \
           xx + q[0] ) \
         * xx + q[1] ) \
         * xx + q[2] ) \
         * xx + q[3] ) \
         * xx + q[4]

    value = sump / sumq

  elif ( 15.0 <= x ):

    if ( xmax < x ):
      value = r8_huge ( )
    else:
#
#  15.0 <= ABS(ARG)
#
      xx = 1.0 / x - rec15

      sump = ((((((  \
                  pp[0] \
           * xx + pp[1] ) \
           * xx + pp[2] ) \
           * xx + pp[3] ) \
           * xx + pp[4] ) \
           * xx + pp[5] ) \
           * xx + pp[6] ) \
           * xx + pp[7]

      sumq = (((((( \
             xx + qq[0] ) \
           * xx + qq[1] ) \
           * xx + qq[2] ) \
           * xx + qq[3] ) \
           * xx + qq[4] ) \
           * xx + qq[5] ) \
           * xx + qq[6]

      value = sump / sumq
#
#  Calculation reformulated to avoid premature overflow.
#
      if ( x <= xmax - 15.0 ):
        a = np.exp ( x )
        b = 1.0
      else:
        a = np.exp ( x - 40.0 )
        b = exp40

      value = ( ( value * a - pp[0] * a ) / np.sqrt ( x ) ) * b

  return value

def bessel_i0_test ( ):

#*****************************************************************************80
#
## BESSEL_I0_TEST tests BESSEL_I0.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    02 March 2016
#
#  Author:
#
#    John Burkardt
#
  import platform
  from bessel_i0_values import bessel_i0_values

  print ( '' )
  print ( 'BESSEL_I0_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  BESSEL_I0 evaluates the Bessel function I0(X)' )
  print ( '' )
  print ( '       X                 Exact F                 I0(X)' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx = bessel_i0_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = bessel_i0 ( x )

    print ( '  %8g  %24.16g  %24.16g' % ( x, fx, fx2 ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'BESSEL_I0_TEST' )
  print ( '  Normal end of execution.' )
  return

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