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

#*****************************************************************************80
#
#% BESSEL_I1 evaluates the Bessel I function of order I.
#
#  Discussion:
#
#    The main computation evaluates slightly modified forms of
#    minimax approximations generated by Blair and Edwards.
#    This transportable program is patterned after the machine-dependent
#    FUNPACK packet NATSI1, but cannot match that version for efficiency
#    or accuracy.  This version uses rational functions that theoretically
#    approximate I-SUB-1(X) to at least 18 significant decimal digits.
#    The accuracy achieved depends on the arithmetic system, the compiler,
#    the intrinsic functions, and proper selection of the machine-dependent
#    constants.
#
#  Machine-dependent constants:
#
#    beta   = Radix for the floating-point system.
#    maxexp = Smallest power of beta that overflows.
#    XMAX =   Largest argument acceptable to BESI1  Solution to
#             equation:
#               EXP(X) * (1-3/(8*X)) / SQRT(2*PI*X) = beta^maxexp
#
#
#    Approximate values for some important machines are:
#
#                            beta       maxexp    XMAX
#
#    CRAY-1        (S.P.)       2         8191    5682.810
#    Cyber 180/855
#      under NOS   (S.P.)       2         1070     745.894
#    IEEE (IBM/XT,
#      SUN, etc.)  (S.P.)       2          128      91.906
#    IEEE (IBM/XT,
#      SUN, etc.)  (D.P.)       2         1024     713.987
#    IBM 3033      (D.P.)      16           63     178.185
#    VAX           (S.P.)       2          127      91.209
#    VAX D-Format  (D.P.)       2          127      91.209
#    VAX G-Format  (D.P.)       2         1023     713.293
#
#  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.
#
#  Reference:
#
#    Blair and Edwards,
#    Chalk River Report AECL-4928,
#    Atomic Energy of Canada, Limited,
#    October, 1974.
#
#  Parameters:
#
#    Input, real ARG, the argument.
#
#    Output, real VALUE, the value of the Bessel I1 function.
#
  import numpy as np
  from r8_epsilon import r8_epsilon
  from r8_huge import r8_huge

  exp40 = 2.353852668370199854E+17
  forty = 40.0E+00
  half = 0.5E+00
  one = 1.0E+00
  one5 = 15.0E+00

  p = np.array ( [ \
    -1.9705291802535139930E-19, \
    -6.5245515583151902910E-16, \
    -1.1928788903603238754E-12, \
    -1.4831904935994647675E-09, \
    -1.3466829827635152875E-06, \
    -9.1746443287817501309E-04, \
    -4.7207090827310162436E-01, \
    -1.8225946631657315931E+02, \
    -5.1894091982308017540E+04, \
    -1.0588550724769347106E+07, \
    -1.4828267606612366099E+09, \
    -1.3357437682275493024E+11, \
    -6.9876779648010090070E+12, \
    -1.7732037840791591320E+14, \
    -1.4577180278143463643E+15 ] )

  pbar = 3.98437500E-01

  pp = np.array ( [ \
    -6.0437159056137600000E-02, \
     4.5748122901933459000E-01, \
    -4.2843766903304806403E-01, \
     9.7356000150886612134E-02, \
    -3.2457723974465568321E-03, \
    -3.6395264712121795296E-04, \
     1.6258661867440836395E-05, \
    -3.6347578404608223492E-07 ] )

  q = np.array ( [ \
    -4.0076864679904189921E+03, \
     7.4810580356655069138E+06, \
    -8.0059518998619764991E+09, \
     4.8544714258273622913E+12, \
    -1.3218168307321442305E+15 ] )

  qq = np.array ( [ \
    -3.8806586721556593450E+00, \
     3.2593714889036996297E+00, \
    -8.5017476463217924408E-01, \
     7.4212010813186530069E-02, \
    -2.2835624489492512649E-03, \
     3.7510433111922824643E-05 ] )

  rec15 = 6.6666666666666666666E-02
  two25 = 225.0E+00
  xmax = 713.987E+00
  zero = 0.0E+00

  x = abs ( arg )
#
#  ABS(ARG) < EPSILON ( ARG )
#
  if ( x < r8_epsilon ( ) ):

    value = half * x
#
#  EPSILON ( ARG ) <= ABS(ARG) < 15.0
#
  elif ( x < one5 ):

    xx = x * x
    sump = p[0]
    for j in range ( 1, 15 ):
      sump = sump * xx + p[j]

    xx = xx - two25

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

    value = ( sump / sumq ) * x

  elif ( xmax < x ):

    value = r8_huge ( )
#
#  15.0 <= ABS(ARG)
#
  else:

    xx = one / 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]

    value = sump / sumq

    if ( xmax - one5 < x ):
      a = np.exp ( x - forty )
      b = exp40
    else:
      a = np.exp ( x )
      b = one

    value = ( ( value * a + pbar * a ) / np.sqrt ( x ) ) * b

  if ( arg < zero ):
    value = -value

  return value

def bessel_i1_test ( ):

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

  print ( '' )
  print ( 'BESSEL_I1_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  BESSEL_I1 evaluates the Bessel function I1(X)' )
  print ( '' )
  print ( '       X                 Exact F                 I1(X)' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx = bessel_i1_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = bessel_i1 ( x )

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

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