function value = 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:
%
%    02 October 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB 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).
%
  persistent ak02cs
  persistent ak0cs
  persistent bk0cs
  persistent ntak0
  persistent ntak02
  persistent ntk0
  persistent xsml

  if ( isempty ( ntk0 ) )

    ak02cs = [ ...
      -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 = [ ...
      -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 = [ ...
      -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 = sqrt ( 4.0 * r8_mach ( 3 ) );

  end

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

  return
end
