function value = r4_besk0e ( x )

%*****************************************************************************80
%
%% R4_BESK0E evaluates the exponentially scaled Bessel function K0(X).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 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.01201869826307592E+00, ...
     -0.00917485269102569E+00, ...
     +0.00014445509317750E+00, ...
     -0.00000401361417543E+00, ...
     +0.00000015678318108E+00, ...
     -0.00000000777011043E+00, ...
     +0.00000000046111825E+00, ...
     -0.00000000003158592E+00, ...
     +0.00000000000243501E+00, ...
     -0.00000000000020743E+00, ...
     +0.00000000000001925E+00, ...
     -0.00000000000000192E+00, ...
     +0.00000000000000020E+00, ...
     -0.00000000000000002E+00 ]';
    ak0cs = [ ...
     -0.07643947903327941E+00, ...
     -0.02235652605699819E+00, ...
     +0.00077341811546938E+00, ...
     -0.00004281006688886E+00, ...
     +0.00000308170017386E+00, ...
     -0.00000026393672220E+00, ...
     +0.00000002563713036E+00, ...
     -0.00000000274270554E+00, ...
     +0.00000000031694296E+00, ...
     -0.00000000003902353E+00, ...
     +0.00000000000506804E+00, ...
     -0.00000000000068895E+00, ...
     +0.00000000000009744E+00, ...
     -0.00000000000001427E+00, ...
     +0.00000000000000215E+00, ...
     -0.00000000000000033E+00, ...
     +0.00000000000000005E+00 ]';
    bk0cs = [ ...
     -0.03532739323390276872E+00, ...
     +0.3442898999246284869E+00, ...
     +0.03597993651536150163E+00, ...
     +0.00126461541144692592E+00, ...
     +0.00002286212103119451E+00, ...
     +0.00000025347910790261E+00, ...
     +0.00000000190451637722E+00, ...
     +0.00000000001034969525E+00, ...
     +0.00000000000004259816E+00, ...
     +0.00000000000000013744E+00, ...
     +0.00000000000000000035E+00 ]';

    ntk0 = r4_inits ( bk0cs, 11, 0.1 * r4_mach ( 3 ) );
    ntak0 = r4_inits ( ak0cs, 17, 0.1 * r4_mach ( 3 ) );
    ntak02 = r4_inits ( ak02cs, 14, 0.1 * r4_mach ( 3 ) );
    xsml = sqrt ( 4.0 * r4_mach ( 3 ) );

  end

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

  return
end