#! /usr/bin/env python
#
def r8_besi0 ( x ):

#*****************************************************************************80
#
## R8_BESI0 evaluates the Bessel function I of order 0 of an R8 argument.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 April 2016
#
#  Author:
#
#    Original FORTRAN77 version by Wayne Fullerton.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Wayne Fullerton,
#    Portable Special Function Routines,
#    in Portability of Numerical Software,
#    edited by Wayne Cowell,
#    Lecture Notes in Computer Science, Volume 57,
#    Springer 1977,
#    ISBN: 978-3-540-08446-4,
#    LC: QA297.W65.
#
#  Parameters:
#
#    Input, real X, the argument.
#
#    Output, real VALUE, the Bessel function I of order 0 of X.
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach
  from sys import exit

  bi0cs = np.array ( [ \
      -0.7660547252839144951081894976243285E-01, \
      +0.1927337953993808269952408750881196E+01, \
      +0.2282644586920301338937029292330415, \
      +0.1304891466707290428079334210691888E-01, \
      +0.4344270900816487451378682681026107E-03, \
      +0.9422657686001934663923171744118766E-05, \
      +0.1434006289510691079962091878179957E-06, \
      +0.1613849069661749069915419719994611E-08, \
      +0.1396650044535669699495092708142522E-10, \
      +0.9579451725505445344627523171893333E-13, \
      +0.5333981859862502131015107744000000E-15, \
      +0.2458716088437470774696785919999999E-17, \
      +0.9535680890248770026944341333333333E-20, \
      +0.3154382039721427336789333333333333E-22, \
      +0.9004564101094637431466666666666666E-25, \
      +0.2240647369123670016000000000000000E-27, \
      +0.4903034603242837333333333333333333E-30, \
      +0.9508172606122666666666666666666666E-33 ] )

  nti0 = r8_inits ( bi0cs, 18, 0.1 * r8_mach ( 3 ) )
  xsml = np.sqrt ( 8.0 * r8_mach ( 3 ) )
  xmax = np.log ( r8_mach ( 2 ) )

  y = abs ( x )

  if ( y <= xsml ):
    value = 1.0
  elif ( y <= 3.0 ):
    value = 2.75 + r8_csevl ( y * y / 4.5 - 1.0, bi0cs, nti0 )
  elif ( y <= xmax ):
    value = np.exp ( y ) * r8_besi0e ( x )
  else:
    print ( '' )
    print ( 'R8_BESI0 - Fatal error!' )
    print ( '  |X| too large.' )
    exit ( 'R8_BESI0 - Fatal error!' )

  return value

def r8_besi0_test ( ):

#*****************************************************************************80
#
## R8_BESI0_TEST tests R8_BESI0.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 April 2016
#
#  Author:
#
#    John Burkardt
#
  import platform
  from bessel_i0_values import bessel_i0_values

  print ( '' )
  print ( 'R8_BESI0_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_BESI0 evaluates the Bessel I0(x) function.' )
  print ( '' )
  print ( '             X       BESI0(X)  R8_BESI0(X)         Diff' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = bessel_i0_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = r8_besi0 ( x )

    print ( '  %14.4f  %14.6g  %14.6g  %14.6g' % ( x, fx1, fx2, abs ( fx1 - fx2 ) ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8_BESI0_TEST:' )
  print ( '  Normal end of execution.' )
  return

def r8_besi0e ( x ):

#*****************************************************************************80
#
## R8_BESI0E evaluates the exponentially scaled Bessel function I0(X).
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 September 2011
#
#  Author:
#
#    Original FORTRAN77 version by Wayne Fullerton.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Wayne Fullerton,
#    Portable Special Function Routines,
#    in Portability of Numerical Software,
#    edited by Wayne Cowell,
#    Lecture Notes in Computer Science, Volume 57,
#    Springer 1977,
#    ISBN: 978-3-540-08446-4,
#    LC: QA297.W65.
#
#  Parameters:
#
#    Input, real X, the argument.
#
#    Output, real VALUE, the exponentially scaled Bessel function I0(X).
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach

  ai02cs = np.array ( [ \
      +0.5449041101410883160789609622680E-01, \
      +0.3369116478255694089897856629799E-02, \
      +0.6889758346916823984262639143011E-04, \
      +0.2891370520834756482966924023232E-05, \
      +0.2048918589469063741827605340931E-06, \
      +0.2266668990498178064593277431361E-07, \
      +0.3396232025708386345150843969523E-08, \
      +0.4940602388224969589104824497835E-09, \
      +0.1188914710784643834240845251963E-10, \
      -0.3149916527963241364538648629619E-10, \
      -0.1321581184044771311875407399267E-10, \
      -0.1794178531506806117779435740269E-11, \
      +0.7180124451383666233671064293469E-12, \
      +0.3852778382742142701140898017776E-12, \
      +0.1540086217521409826913258233397E-13, \
      -0.4150569347287222086626899720156E-13, \
      -0.9554846698828307648702144943125E-14, \
      +0.3811680669352622420746055355118E-14, \
      +0.1772560133056526383604932666758E-14, \
      -0.3425485619677219134619247903282E-15, \
      -0.2827623980516583484942055937594E-15, \
      +0.3461222867697461093097062508134E-16, \
      +0.4465621420296759999010420542843E-16, \
      -0.4830504485944182071255254037954E-17, \
      -0.7233180487874753954562272409245E-17, \
      +0.9921475412173698598880460939810E-18, \
      +0.1193650890845982085504399499242E-17, \
      -0.2488709837150807235720544916602E-18, \
      -0.1938426454160905928984697811326E-18, \
      +0.6444656697373443868783019493949E-19, \
      +0.2886051596289224326481713830734E-19, \
      -0.1601954907174971807061671562007E-19, \
      -0.3270815010592314720891935674859E-20, \
      +0.3686932283826409181146007239393E-20, \
      +0.1268297648030950153013595297109E-22, \
      -0.7549825019377273907696366644101E-21, \
      +0.1502133571377835349637127890534E-21, \
      +0.1265195883509648534932087992483E-21, \
      -0.6100998370083680708629408916002E-22, \
      -0.1268809629260128264368720959242E-22, \
      +0.1661016099890741457840384874905E-22, \
      -0.1585194335765885579379705048814E-23, \
      -0.3302645405968217800953817667556E-23, \
      +0.1313580902839239781740396231174E-23, \
      +0.3689040246671156793314256372804E-24, \
      -0.4210141910461689149219782472499E-24, \
      +0.4791954591082865780631714013730E-25, \
      +0.8459470390221821795299717074124E-25, \
      -0.4039800940872832493146079371810E-25, \
      -0.6434714653650431347301008504695E-26, \
      +0.1225743398875665990344647369905E-25, \
      -0.2934391316025708923198798211754E-26, \
      -0.1961311309194982926203712057289E-26, \
      +0.1503520374822193424162299003098E-26, \
      -0.9588720515744826552033863882069E-28, \
      -0.3483339380817045486394411085114E-27, \
      +0.1690903610263043673062449607256E-27, \
      +0.1982866538735603043894001157188E-28, \
      -0.5317498081491816214575830025284E-28, \
      +0.1803306629888392946235014503901E-28, \
      +0.6213093341454893175884053112422E-29, \
      -0.7692189292772161863200728066730E-29, \
      +0.1858252826111702542625560165963E-29, \
      +0.1237585142281395724899271545541E-29, \
      -0.1102259120409223803217794787792E-29, \
      +0.1886287118039704490077874479431E-30, \
      +0.2160196872243658913149031414060E-30, \
      -0.1605454124919743200584465949655E-30, \
      +0.1965352984594290603938848073318E-31 ] )

  ai0cs = np.array ( [ \
      +0.7575994494023795942729872037438E-01, \
      +0.7591380810823345507292978733204E-02, \
      +0.4153131338923750501863197491382E-03, \
      +0.1070076463439073073582429702170E-04, \
      -0.7901179979212894660750319485730E-05, \
      -0.7826143501438752269788989806909E-06, \
      +0.2783849942948870806381185389857E-06, \
      +0.8252472600612027191966829133198E-08, \
      -0.1204463945520199179054960891103E-07, \
      +0.1559648598506076443612287527928E-08, \
      +0.2292556367103316543477254802857E-09, \
      -0.1191622884279064603677774234478E-09, \
      +0.1757854916032409830218331247743E-10, \
      +0.1128224463218900517144411356824E-11, \
      -0.1146848625927298877729633876982E-11, \
      +0.2715592054803662872643651921606E-12, \
      -0.2415874666562687838442475720281E-13, \
      -0.6084469888255125064606099639224E-14, \
      +0.3145705077175477293708360267303E-14, \
      -0.7172212924871187717962175059176E-15, \
      +0.7874493403454103396083909603327E-16, \
      +0.1004802753009462402345244571839E-16, \
      -0.7566895365350534853428435888810E-17, \
      +0.2150380106876119887812051287845E-17, \
      -0.3754858341830874429151584452608E-18, \
      +0.2354065842226992576900757105322E-19, \
      +0.1114667612047928530226373355110E-19, \
      -0.5398891884396990378696779322709E-20, \
      +0.1439598792240752677042858404522E-20, \
      -0.2591916360111093406460818401962E-21, \
      +0.2238133183998583907434092298240E-22, \
      +0.5250672575364771172772216831999E-23, \
      -0.3249904138533230784173432285866E-23, \
      +0.9924214103205037927857284710400E-24, \
      -0.2164992254244669523146554299733E-24, \
      +0.3233609471943594083973332991999E-25, \
      -0.1184620207396742489824733866666E-26, \
      -0.1281671853950498650548338687999E-26, \
      +0.5827015182279390511605568853333E-27, \
      -0.1668222326026109719364501503999E-27, \
      +0.3625309510541569975700684800000E-28, \
      -0.5733627999055713589945958399999E-29, \
      +0.3736796722063098229642581333333E-30, \
      +0.1602073983156851963365512533333E-30, \
      -0.8700424864057229884522495999999E-31, \
      +0.2741320937937481145603413333333E-31 ] )

  bi0cs = np.array ( [ \
      -0.7660547252839144951081894976243285E-01, \
      +0.1927337953993808269952408750881196E+01, \
      +0.2282644586920301338937029292330415, \
      +0.1304891466707290428079334210691888E-01, \
      +0.4344270900816487451378682681026107E-03, \
      +0.9422657686001934663923171744118766E-05, \
      +0.1434006289510691079962091878179957E-06, \
      +0.1613849069661749069915419719994611E-08, \
      +0.1396650044535669699495092708142522E-10, \
      +0.9579451725505445344627523171893333E-13, \
      +0.5333981859862502131015107744000000E-15, \
      +0.2458716088437470774696785919999999E-17, \
      +0.9535680890248770026944341333333333E-20, \
      +0.3154382039721427336789333333333333E-22, \
      +0.9004564101094637431466666666666666E-25, \
      +0.2240647369123670016000000000000000E-27, \
      +0.4903034603242837333333333333333333E-30, \
      +0.9508172606122666666666666666666666E-33 ] )

  eta = 0.1 * r8_mach ( 3 )
  nti0 = r8_inits ( bi0cs, 18, eta )
  ntai0 = r8_inits ( ai0cs, 46, eta )
  ntai02 = r8_inits ( ai02cs, 69, eta )
  xsml = np.sqrt ( 8.0 * r8_mach ( 3 ) )

  y = abs ( x )

  if ( y <= xsml ):
    value = 1.0
  elif ( y <= 3.0 ):
    value = np.exp ( - y ) * ( 2.75 \
      + r8_csevl ( y * y / 4.5 - 1.0, bi0cs, nti0 ) )
  elif ( y <= 8.0 ):
    value = ( 0.375 + r8_csevl ( ( 48.0 / y - 11.0 ) / 5.0, \
      ai0cs, ntai0 ) ) / np.sqrt ( y )
  else:
    value = ( 0.375 + r8_csevl ( 16.0 / y - 1.0, ai02cs, ntai02 ) ) \
      / np.sqrt ( y )

  return value

def r8_besk0 ( x ):

#*****************************************************************************80
#
## R8_BESK0 evaluates the Bessel function K of order 0 of an R8 argument.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 April 2016
#
#  Author:
#
#    Original FORTRAN77 version by Wayne Fullerton.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Wayne Fullerton,
#    Portable Special Function Routines,
#    in Portability of Numerical Software,
#    edited by Wayne Cowell,
#    Lecture Notes in Computer Science, Volume 57,
#    Springer 1977,
#    ISBN: 978-3-540-08446-4,
#    LC: QA297.W65.
#
#  Parameters:
#
#    Input, real X, the argument.
#
#    Output, real VALUE, the Bessel function K of order 0 of X.
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach
  from sys import exit

  bk0cs = np.array ( [ \
      -0.353273932339027687201140060063153E-01, \
      +0.344289899924628486886344927529213, \
      +0.359799365153615016265721303687231E-01, \
      +0.126461541144692592338479508673447E-02, \
      +0.228621210311945178608269830297585E-04, \
      +0.253479107902614945730790013428354E-06, \
      +0.190451637722020885897214059381366E-08, \
      +0.103496952576336245851008317853089E-10, \
      +0.425981614279108257652445327170133E-13, \
      +0.137446543588075089694238325440000E-15, \
      +0.357089652850837359099688597333333E-18, \
      +0.763164366011643737667498666666666E-21, \
      +0.136542498844078185908053333333333E-23, \
      +0.207527526690666808319999999999999E-26, \
      +0.271281421807298560000000000000000E-29, \
      +0.308259388791466666666666666666666E-32 ] )

  ntk0 = r8_inits ( bk0cs, 16, 0.1 * r8_mach ( 3 ) )
  xsml = np.sqrt ( 4.0 * r8_mach ( 3 ) )
  xmax = - np.log ( r8_mach ( 1 ) )
  xmax = xmax - 0.5 * xmax * np.log ( xmax ) / ( xmax + 0.5 )

  if ( x <= 0.0 ):
    print ( '' )
    print ( 'R8_BESK0 = Fatal error!' )
    print ( '  X <= 0.' )
    exit ( 'R8_BESK0 = Fatal error!' )
  elif ( x <= xsml ):
    y = 0.0
    value = - np.log ( 0.5 * x ) * r8_besi0 ( x ) \
      - 0.25 + r8_csevl ( 0.5 * y - 1.0, bk0cs, ntk0 )
  elif ( x <= 2.0 ):
    y = x * x
    value = - np.log ( 0.5 * x ) * r8_besi0 ( x ) \
      - 0.25 + r8_csevl ( 0.5 * y - 1.0, bk0cs, ntk0 )
  elif ( x <= xmax ):
    value = np.exp ( - x ) * r8_besk0e ( x )
  else:
    value = 0.0

  return value

def r8_besk0_test ( ):

#*****************************************************************************80
#
## R8_BESK0_TEST tests R8_BESK0.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 April 2016
#
#  Author:
#
#    John Burkardt
#
  import platform
  from bessel_k0_values import bessel_k0_values

  print ( '' )
  print ( 'R8_BESK0_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_BESK0 evaluates Bessel functions K0(X)' )
  print ( '' )
  print ( '             X       BESK0(X)  R8_BESK0(X)         Diff' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = bessel_k0_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = r8_besk0 ( x )

    print ( '  %14.4f  %14.6g  %14.6g  %14.6g' % ( x, fx1, fx2, abs ( fx1 - fx2 ) ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8_BESK0_TEST:' )
  print ( '  Normal end of execution.' )
  return

def r8_besk0e ( x ):

#*****************************************************************************80
#
## R8_BESK0E evaluates the exponentially scaled Bessel function K0(X).
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 April 2016
#
#  Author:
#
#    Original FORTRAN77 version by Wayne Fullerton.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Wayne Fullerton,
#    Portable Special Function Routines,
#    in Portability of Numerical Software,
#    edited by Wayne Cowell,
#    Lecture Notes in Computer Science, Volume 57,
#    Springer 1977,
#    ISBN: 978-3-540-08446-4,
#    LC: QA297.W65.
#
#  Parameters:
#
#    Input, real X, the argument.
#
#    Output, real VALUE, the exponentially scaled Bessel function K0(X).
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach
  from sys import exit

  ak02cs = np.array ( [ \
      -0.1201869826307592239839346212452E-01, \
      -0.9174852691025695310652561075713E-02, \
      +0.1444550931775005821048843878057E-03, \
      -0.4013614175435709728671021077879E-05, \
      +0.1567831810852310672590348990333E-06, \
      -0.7770110438521737710315799754460E-08, \
      +0.4611182576179717882533130529586E-09, \
      -0.3158592997860565770526665803309E-10, \
      +0.2435018039365041127835887814329E-11, \
      -0.2074331387398347897709853373506E-12, \
      +0.1925787280589917084742736504693E-13, \
      -0.1927554805838956103600347182218E-14, \
      +0.2062198029197818278285237869644E-15, \
      -0.2341685117579242402603640195071E-16, \
      +0.2805902810643042246815178828458E-17, \
      -0.3530507631161807945815482463573E-18, \
      +0.4645295422935108267424216337066E-19, \
      -0.6368625941344266473922053461333E-20, \
      +0.9069521310986515567622348800000E-21, \
      -0.1337974785423690739845005311999E-21, \
      +0.2039836021859952315522088960000E-22, \
      -0.3207027481367840500060869973333E-23, \
      +0.5189744413662309963626359466666E-24, \
      -0.8629501497540572192964607999999E-25, \
      +0.1472161183102559855208038400000E-25, \
      -0.2573069023867011283812351999999E-26, \
      +0.4601774086643516587376640000000E-27, \
      -0.8411555324201093737130666666666E-28, \
      +0.1569806306635368939301546666666E-28, \
      -0.2988226453005757788979199999999E-29, \
      +0.5796831375216836520618666666666E-30, \
      -0.1145035994347681332155733333333E-30, \
      +0.2301266594249682802005333333333E-31 ] )

  ak0cs = np.array ( [ \
      -0.7643947903327941424082978270088E-01, \
      -0.2235652605699819052023095550791E-01, \
      +0.7734181154693858235300618174047E-03, \
      -0.4281006688886099464452146435416E-04, \
      +0.3081700173862974743650014826660E-05, \
      -0.2639367222009664974067448892723E-06, \
      +0.2563713036403469206294088265742E-07, \
      -0.2742705549900201263857211915244E-08, \
      +0.3169429658097499592080832873403E-09, \
      -0.3902353286962184141601065717962E-10, \
      +0.5068040698188575402050092127286E-11, \
      -0.6889574741007870679541713557984E-12, \
      +0.9744978497825917691388201336831E-13, \
      -0.1427332841884548505389855340122E-13, \
      +0.2156412571021463039558062976527E-14, \
      -0.3349654255149562772188782058530E-15, \
      +0.5335260216952911692145280392601E-16, \
      -0.8693669980890753807639622378837E-17, \
      +0.1446404347862212227887763442346E-17, \
      -0.2452889825500129682404678751573E-18, \
      +0.4233754526232171572821706342400E-19, \
      -0.7427946526454464195695341294933E-20, \
      +0.1323150529392666866277967462400E-20, \
      -0.2390587164739649451335981465599E-21, \
      +0.4376827585923226140165712554666E-22, \
      -0.8113700607345118059339011413333E-23, \
      +0.1521819913832172958310378154666E-23, \
      -0.2886041941483397770235958613333E-24, \
      +0.5530620667054717979992610133333E-25, \
      -0.1070377329249898728591633066666E-25, \
      +0.2091086893142384300296328533333E-26, \
      -0.4121713723646203827410261333333E-27, \
      +0.8193483971121307640135680000000E-28, \
      -0.1642000275459297726780757333333E-28, \
      +0.3316143281480227195890346666666E-29, \
      -0.6746863644145295941085866666666E-30, \
      +0.1382429146318424677635413333333E-30, \
      -0.2851874167359832570811733333333E-31 ] )

  bk0cs = np.array ( [ \
      -0.353273932339027687201140060063153E-01, \
      +0.344289899924628486886344927529213, \
      +0.359799365153615016265721303687231E-01, \
      +0.126461541144692592338479508673447E-02, \
      +0.228621210311945178608269830297585E-04, \
      +0.253479107902614945730790013428354E-06, \
      +0.190451637722020885897214059381366E-08, \
      +0.103496952576336245851008317853089E-10, \
      +0.425981614279108257652445327170133E-13, \
      +0.137446543588075089694238325440000E-15, \
      +0.357089652850837359099688597333333E-18, \
      +0.763164366011643737667498666666666E-21, \
      +0.136542498844078185908053333333333E-23, \
      +0.207527526690666808319999999999999E-26, \
      +0.271281421807298560000000000000000E-29, \
      +0.308259388791466666666666666666666E-32 ] )

  eta = 0.1 * r8_mach ( 3 )
  ntk0 = r8_inits ( bk0cs, 16, eta )
  ntak0 = r8_inits ( ak0cs, 38, eta )
  ntak02 = r8_inits ( ak02cs, 33, eta )
  xsml = np.sqrt ( 4.0 * r8_mach ( 3 ) )

  if ( x <= 0.0 ):
    print ( '' )
    print ( 'R8_BESK0E = Fatal error!' )
    print ( '  X <= 0.' )
    exit ( 'R8_BESK0E = Fatal error!' )
  elif ( x <= xsml ):
    y = 0.0
    value = np.exp ( x ) * ( - np.log ( 0.5 * x ) * r8_besi0 ( x ) - 0.25 \
      + r8_csevl ( 0.5 * y - 1.0, bk0cs, ntk0 ) )
  elif ( x <= 2.0 ):
    y = x * x
    value = np.exp ( x ) * ( - np.log ( 0.5 * x ) * r8_besi0 ( x ) - 0.25 \
      + r8_csevl ( 0.5 * y - 1.0, bk0cs, ntk0 ) )
  elif ( x <= 8.0 ):
    value = ( 1.25 \
      + r8_csevl ( ( 16.0 / x - 5.0 ) / 3.0, ak0cs, ntak0 ) ) / np.sqrt ( x )
  else:
    value = ( 1.25 + r8_csevl ( 16.0 / x - 1.0, ak02cs, ntak02 ) ) \
      / np.sqrt ( x )

  return value

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  r8_besi0_test ( )
  r8_besk0_test ( )
  timestamp ( )