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

#*****************************************************************************80
#
## R8_BESI1 evaluates the Bessel function I of order 1 of an R8 argument.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    24 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 1 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

  bi1cs = np.array ( [ \
      -0.19717132610998597316138503218149E-02, \
      +0.40734887667546480608155393652014, \
      +0.34838994299959455866245037783787E-01, \
      +0.15453945563001236038598401058489E-02, \
      +0.41888521098377784129458832004120E-04, \
      +0.76490267648362114741959703966069E-06, \
      +0.10042493924741178689179808037238E-07, \
      +0.99322077919238106481371298054863E-10, \
      +0.76638017918447637275200171681349E-12, \
      +0.47414189238167394980388091948160E-14, \
      +0.24041144040745181799863172032000E-16, \
      +0.10171505007093713649121100799999E-18, \
      +0.36450935657866949458491733333333E-21, \
      +0.11205749502562039344810666666666E-23, \
      +0.29875441934468088832000000000000E-26, \
      +0.69732310939194709333333333333333E-29, \
      +0.14367948220620800000000000000000E-31 ] )

  nti1 = r8_inits ( bi1cs, 17, 0.1 * r8_mach ( 3 ) )
  xmin = 2.0 * r8_mach ( 1 )
  xsml = np.sqrt ( 8.0 * r8_mach ( 3 ) )
  xmax = np.log ( r8_mach ( 2 ) )

  y = abs ( x )

  if ( y <= xmin ):
    value = 0.0
  elif ( y <= xsml ):
    value = 0.5 * x
  elif ( y <= 3.0 ):
    value = x * ( 0.875 + r8_csevl ( y * y / 4.5 - 1.0, bi1cs, nti1 ) )
  elif ( y <= xmax ):
    value = np.exp ( y ) * r8_besi1e ( x )
  else:
    print ( '' )
    print ( 'R8_BESI1 - Fatal error!' )
    print ( '  Result overflows.' )
    exit ( 'R8_BESI1 - Fatal error!' )

  return value

def r8_besi1_test ( ):

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

  print ( '' )
  print ( 'R8_BESI1_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_BESI1 evaluates the Bessel I1(x) function' )
  print ( '' )
  print ( '             X       BESI1(X)  R8_BESI1(X)         Diff' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = bessel_i1_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = r8_besi1 ( x )

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

def r8_besi1e ( x ):

#*****************************************************************************80
#
## R8_BESI1E evaluates the exponentially scaled Bessel function I1(X).
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    24 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 I1(X).
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach

  ai12cs = np.array ( [ \
    +0.2857623501828012047449845948469E-01, \
    -0.9761097491361468407765164457302E-02, \
    -0.1105889387626237162912569212775E-03, \
    -0.3882564808877690393456544776274E-05, \
    -0.2512236237870208925294520022121E-06, \
    -0.2631468846889519506837052365232E-07, \
    -0.3835380385964237022045006787968E-08, \
    -0.5589743462196583806868112522229E-09, \
    -0.1897495812350541234498925033238E-10, \
    +0.3252603583015488238555080679949E-10, \
    +0.1412580743661378133163366332846E-10, \
    +0.2035628544147089507224526136840E-11, \
    -0.7198551776245908512092589890446E-12, \
    -0.4083551111092197318228499639691E-12, \
    -0.2101541842772664313019845727462E-13, \
    +0.4272440016711951354297788336997E-13, \
    +0.1042027698412880276417414499948E-13, \
    -0.3814403072437007804767072535396E-14, \
    -0.1880354775510782448512734533963E-14, \
    +0.3308202310920928282731903352405E-15, \
    +0.2962628997645950139068546542052E-15, \
    -0.3209525921993423958778373532887E-16, \
    -0.4650305368489358325571282818979E-16, \
    +0.4414348323071707949946113759641E-17, \
    +0.7517296310842104805425458080295E-17, \
    -0.9314178867326883375684847845157E-18, \
    -0.1242193275194890956116784488697E-17, \
    +0.2414276719454848469005153902176E-18, \
    +0.2026944384053285178971922860692E-18, \
    -0.6394267188269097787043919886811E-19, \
    -0.3049812452373095896084884503571E-19, \
    +0.1612841851651480225134622307691E-19, \
    +0.3560913964309925054510270904620E-20, \
    -0.3752017947936439079666828003246E-20, \
    -0.5787037427074799345951982310741E-22, \
    +0.7759997511648161961982369632092E-21, \
    -0.1452790897202233394064459874085E-21, \
    -0.1318225286739036702121922753374E-21, \
    +0.6116654862903070701879991331717E-22, \
    +0.1376279762427126427730243383634E-22, \
    -0.1690837689959347884919839382306E-22, \
    +0.1430596088595433153987201085385E-23, \
    +0.3409557828090594020405367729902E-23, \
    -0.1309457666270760227845738726424E-23, \
    -0.3940706411240257436093521417557E-24, \
    +0.4277137426980876580806166797352E-24, \
    -0.4424634830982606881900283123029E-25, \
    -0.8734113196230714972115309788747E-25, \
    +0.4045401335683533392143404142428E-25, \
    +0.7067100658094689465651607717806E-26, \
    -0.1249463344565105223002864518605E-25, \
    +0.2867392244403437032979483391426E-26, \
    +0.2044292892504292670281779574210E-26, \
    -0.1518636633820462568371346802911E-26, \
    +0.8110181098187575886132279107037E-28, \
    +0.3580379354773586091127173703270E-27, \
    -0.1692929018927902509593057175448E-27, \
    -0.2222902499702427639067758527774E-28, \
    +0.5424535127145969655048600401128E-28, \
    -0.1787068401578018688764912993304E-28, \
    -0.6565479068722814938823929437880E-29, \
    +0.7807013165061145280922067706839E-29, \
    -0.1816595260668979717379333152221E-29, \
    -0.1287704952660084820376875598959E-29, \
    +0.1114548172988164547413709273694E-29, \
    -0.1808343145039336939159368876687E-30, \
    -0.2231677718203771952232448228939E-30, \
    +0.1619029596080341510617909803614E-30, \
    -0.1834079908804941413901308439210E-31 ] )

  ai1cs = np.array ( [ \
    -0.2846744181881478674100372468307E-01, \
    -0.1922953231443220651044448774979E-01, \
    -0.6115185857943788982256249917785E-03, \
    -0.2069971253350227708882823777979E-04, \
    +0.8585619145810725565536944673138E-05, \
    +0.1049498246711590862517453997860E-05, \
    -0.2918338918447902202093432326697E-06, \
    -0.1559378146631739000160680969077E-07, \
    +0.1318012367144944705525302873909E-07, \
    -0.1448423418183078317639134467815E-08, \
    -0.2908512243993142094825040993010E-09, \
    +0.1266388917875382387311159690403E-09, \
    -0.1664947772919220670624178398580E-10, \
    -0.1666653644609432976095937154999E-11, \
    +0.1242602414290768265232168472017E-11, \
    -0.2731549379672432397251461428633E-12, \
    +0.2023947881645803780700262688981E-13, \
    +0.7307950018116883636198698126123E-14, \
    -0.3332905634404674943813778617133E-14, \
    +0.7175346558512953743542254665670E-15, \
    -0.6982530324796256355850629223656E-16, \
    -0.1299944201562760760060446080587E-16, \
    +0.8120942864242798892054678342860E-17, \
    -0.2194016207410736898156266643783E-17, \
    +0.3630516170029654848279860932334E-18, \
    -0.1695139772439104166306866790399E-19, \
    -0.1288184829897907807116882538222E-19, \
    +0.5694428604967052780109991073109E-20, \
    -0.1459597009090480056545509900287E-20, \
    +0.2514546010675717314084691334485E-21, \
    -0.1844758883139124818160400029013E-22, \
    -0.6339760596227948641928609791999E-23, \
    +0.3461441102031011111108146626560E-23, \
    -0.1017062335371393547596541023573E-23, \
    +0.2149877147090431445962500778666E-24, \
    -0.3045252425238676401746206173866E-25, \
    +0.5238082144721285982177634986666E-27, \
    +0.1443583107089382446416789503999E-26, \
    -0.6121302074890042733200670719999E-27, \
    +0.1700011117467818418349189802666E-27, \
    -0.3596589107984244158535215786666E-28, \
    +0.5448178578948418576650513066666E-29, \
    -0.2731831789689084989162564266666E-30, \
    -0.1858905021708600715771903999999E-30, \
    +0.9212682974513933441127765333333E-31, \
    -0.2813835155653561106370833066666E-31 ] )

  bi1cs = np.array ( [ \
    -0.19717132610998597316138503218149E-02, \
    +0.40734887667546480608155393652014E+00, \
    +0.34838994299959455866245037783787E-01, \
    +0.15453945563001236038598401058489E-02, \
    +0.41888521098377784129458832004120E-04, \
    +0.76490267648362114741959703966069E-06, \
    +0.10042493924741178689179808037238E-07, \
    +0.99322077919238106481371298054863E-10, \
    +0.76638017918447637275200171681349E-12, \
    +0.47414189238167394980388091948160E-14, \
    +0.24041144040745181799863172032000E-16, \
    +0.10171505007093713649121100799999E-18, \
    +0.36450935657866949458491733333333E-21, \
    +0.11205749502562039344810666666666E-23, \
    +0.29875441934468088832000000000000E-26, \
    +0.69732310939194709333333333333333E-29, \
    +0.14367948220620800000000000000000E-31 ] )

  eta = 0.1 * r8_mach ( 3 )
  nti1 = r8_inits ( bi1cs, 17, eta )
  ntai1 = r8_inits ( ai1cs, 46, eta )
  ntai12 = r8_inits ( ai12cs, 69, eta )
  xmin = 2.0 * r8_mach ( 1 )
  xsml = np.sqrt ( 8.0 * r8_mach ( 3 ) )

  y = abs ( x )

  if ( y <= xmin ):
    value = 0.0
  elif ( y <= xsml ):
    value = 0.5 * x * np.exp ( - y )
  elif ( y <= 3.0 ):
    value = x * ( 0.875 + r8_csevl ( y * y / 4.5 - 1.0, bi1cs, nti1 ) ) \
      * np.exp ( - y )
  elif ( y <= 8.0 ):
    value = ( 0.375 + r8_csevl ( ( 48.0 / y - 11.0) / 5.0, \
      ai1cs, ntai1 ) ) / np.sqrt ( y )
    if ( x < 0.0 ):
      value = - value
  else:
    value = ( 0.375 + r8_csevl ( 16.0 / y - 1.0, ai12cs, ntai12 ) ) \
      / np.sqrt ( y )
    if ( x < 0.0 ):
      value = - value

  return value

def r8_besk1 ( x ):

#*****************************************************************************80
#
## R8_BESK1 evaluates the Bessel function K of order 1 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 1 of X.
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach

  bk1cs = np.array ( [ \
      +0.25300227338947770532531120868533E-01, \
      -0.35315596077654487566723831691801, \
      -0.12261118082265714823479067930042, \
      -0.69757238596398643501812920296083E-02, \
      -0.17302889575130520630176507368979E-03, \
      -0.24334061415659682349600735030164E-05, \
      -0.22133876307347258558315252545126E-07, \
      -0.14114883926335277610958330212608E-09, \
      -0.66669016941993290060853751264373E-12, \
      -0.24274498505193659339263196864853E-14, \
      -0.70238634793862875971783797120000E-17, \
      -0.16543275155100994675491029333333E-19, \
      -0.32338347459944491991893333333333E-22, \
      -0.53312750529265274999466666666666E-25, \
      -0.75130407162157226666666666666666E-28, \
      -0.91550857176541866666666666666666E-31 ] )

  ntk1 = r8_inits ( bk1cs, 16, 0.1 * r8_mach ( 3 ) )
  xmin = np.exp ( max ( np.log ( r8_mach ( 1 ) ), \
    - np.log ( r8_mach ( 2 ) ) ) + 0.01 )
  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 ) - 0.01

  if ( x <= 0.0 ):
    print ( '' )
    print ( 'R8_BESK1 = Fatal error!' )
    print ( '  X <= 0.' )
    exit ( 'R8_BESK1 - Fatal error!' )
  elif ( x <= xsml ):
    y = 0.0
    value = np.log ( 0.5 * x ) * r8_besi1 ( x ) + ( 0.75 \
      + r8_csevl ( 0.5 * y - 1.0, bk1cs, ntk1 ) ) / x
  elif ( x <= 2.0 ):
    y = x * x
    value = np.log ( 0.5 * x ) * r8_besi1 ( x ) + ( 0.75 \
      + r8_csevl ( 0.5 * y - 1.0, bk1cs, ntk1 ) ) / x
  elif ( x <= xmax ):
    value = np.exp ( - x ) * r8_besk1e ( x )
  else:
    value = 0.0

  return value

def r8_besk1_test ( ):

#*****************************************************************************80
#
#% R8_BESK1_TEST tests R8_BESK1.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 April 2016
#
#  Author:
#
#    John Burkardt
#
  import platform
  from bessel_k1_values import bessel_k1_values

  print ( '' )
  print ( 'R8_BESK1_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_BESK1 evaluates Bessel functions K1(x)' )
  print ( '' )
  print ( '             X       BESK1(X)  R8_BESK1(X)         Diff' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = bessel_k1_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = r8_besk1 ( x )

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

def r8_besk1e ( x ):

#*****************************************************************************80
#
## R8_BESK1E evaluates the exponentially scaled Bessel function K1(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 K1(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

  ak12cs = np.array ( [ \
      +0.6379308343739001036600488534102E-01, \
      +0.2832887813049720935835030284708E-01, \
      -0.2475370673905250345414545566732E-03, \
      +0.5771972451607248820470976625763E-05, \
      -0.2068939219536548302745533196552E-06, \
      +0.9739983441381804180309213097887E-08, \
      -0.5585336140380624984688895511129E-09, \
      +0.3732996634046185240221212854731E-10, \
      -0.2825051961023225445135065754928E-11, \
      +0.2372019002484144173643496955486E-12, \
      -0.2176677387991753979268301667938E-13, \
      +0.2157914161616032453939562689706E-14, \
      -0.2290196930718269275991551338154E-15, \
      +0.2582885729823274961919939565226E-16, \
      -0.3076752641268463187621098173440E-17, \
      +0.3851487721280491597094896844799E-18, \
      -0.5044794897641528977117282508800E-19, \
      +0.6888673850418544237018292223999E-20, \
      -0.9775041541950118303002132480000E-21, \
      +0.1437416218523836461001659733333E-21, \
      -0.2185059497344347373499733333333E-22, \
      +0.3426245621809220631645388800000E-23, \
      -0.5531064394246408232501248000000E-24, \
      +0.9176601505685995403782826666666E-25, \
      -0.1562287203618024911448746666666E-25, \
      +0.2725419375484333132349439999999E-26, \
      -0.4865674910074827992378026666666E-27, \
      +0.8879388552723502587357866666666E-28, \
      -0.1654585918039257548936533333333E-28, \
      +0.3145111321357848674303999999999E-29, \
      -0.6092998312193127612416000000000E-30, \
      +0.1202021939369815834623999999999E-30, \
      -0.2412930801459408841386666666666E-31 ] )

  ak1cs = np.array ( [ \
      +0.27443134069738829695257666227266, \
      +0.75719899531993678170892378149290E-01, \
      -0.14410515564754061229853116175625E-02, \
      +0.66501169551257479394251385477036E-04, \
      -0.43699847095201407660580845089167E-05, \
      +0.35402774997630526799417139008534E-06, \
      -0.33111637792932920208982688245704E-07, \
      +0.34459775819010534532311499770992E-08, \
      -0.38989323474754271048981937492758E-09, \
      +0.47208197504658356400947449339005E-10, \
      -0.60478356628753562345373591562890E-11, \
      +0.81284948748658747888193837985663E-12, \
      -0.11386945747147891428923915951042E-12, \
      +0.16540358408462282325972948205090E-13, \
      -0.24809025677068848221516010440533E-14, \
      +0.38292378907024096948429227299157E-15, \
      -0.60647341040012418187768210377386E-16, \
      +0.98324256232648616038194004650666E-17, \
      -0.16284168738284380035666620115626E-17, \
      +0.27501536496752623718284120337066E-18, \
      -0.47289666463953250924281069568000E-19, \
      +0.82681500028109932722392050346666E-20, \
      -0.14681405136624956337193964885333E-20, \
      +0.26447639269208245978085894826666E-21, \
      -0.48290157564856387897969868800000E-22, \
      +0.89293020743610130180656332799999E-23, \
      -0.16708397168972517176997751466666E-23, \
      +0.31616456034040694931368618666666E-24, \
      -0.60462055312274989106506410666666E-25, \
      +0.11678798942042732700718421333333E-25, \
      -0.22773741582653996232867840000000E-26, \
      +0.44811097300773675795305813333333E-27, \
      -0.88932884769020194062336000000000E-28, \
      +0.17794680018850275131392000000000E-28, \
      -0.35884555967329095821994666666666E-29, \
      +0.72906290492694257991679999999999E-30, \
      -0.14918449845546227073024000000000E-30, \
      +0.30736573872934276300799999999999E-31 ] )

  bk1cs = np.array ( [ \
      +0.25300227338947770532531120868533E-01, \
      -0.35315596077654487566723831691801, \
      -0.12261118082265714823479067930042, \
      -0.69757238596398643501812920296083E-02, \
      -0.17302889575130520630176507368979E-03, \
      -0.24334061415659682349600735030164E-05, \
      -0.22133876307347258558315252545126E-07, \
      -0.14114883926335277610958330212608E-09, \
      -0.66669016941993290060853751264373E-12, \
      -0.24274498505193659339263196864853E-14, \
      -0.70238634793862875971783797120000E-17, \
      -0.16543275155100994675491029333333E-19, \
      -0.32338347459944491991893333333333E-22, \
      -0.53312750529265274999466666666666E-25, \
      -0.75130407162157226666666666666666E-28, \
      -0.91550857176541866666666666666666E-31 ] )

  eta = 0.1 * r8_mach ( 3 )
  ntk1 = r8_inits ( bk1cs, 16, eta )
  ntak1 = r8_inits ( ak1cs, 38, eta )
  ntak12 = r8_inits ( ak12cs, 33, eta )
  xmin = np.exp ( max ( np.log ( r8_mach ( 1 ) ), \
    - np.log ( r8_mach ( 2 ) ) ) + 0.01 )
  xsml = np.sqrt ( 4.0 * r8_mach ( 3 ) )

  if ( x <= 0.0 ):
    print ( '' )
    print ( 'R8_BESK1E = Fatal error!' )
    print ( '  X <= 0.' )
    exit ( 'R8_BESK1E - Fatal error!' )
  elif ( x <= xsml ):
    y = 0.0
    value = np.exp ( x ) * ( np.log ( 0.5 * x ) * r8_besi1 ( x ) \
      + ( 0.75 + r8_csevl ( 0.5 * y - 1.0, bk1cs, ntk1 ) ) / x )
  elif ( x <= 2.0 ):
    y = x * x
    value = np.exp ( x ) * ( np.log ( 0.5 * x ) * r8_besi1 ( x ) \
      + ( 0.75 + r8_csevl ( 0.5 * y - 1.0, bk1cs, ntk1 ) ) / x )
  elif ( x <= 8.0 ):
    value = ( 1.25 + r8_csevl ( ( 16.0 / x - 5.0 ) / 3.0, ak1cs, \
      ntak1 ) ) / np.sqrt ( x )
  else:
    value = ( 1.25 + r8_csevl ( 16.0 / x - 1.0, ak12cs, ntak12 ) ) / np.sqrt ( x )

  return value

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  r8_besi1_test ( )
  r8_besk1_test ( )
  timestamp ( )