function value = r4_besk1e ( x )

%*****************************************************************************80
%
%% R4_BESK1E evaluates the exponentially scaled Bessel function K1(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 K1(X).
%
  persistent ak12cs
  persistent ak1cs
  persistent bk1cs
  persistent ntak1
  persistent ntak12
  persistent ntk1
  persistent xmin
  persistent xsml

  if ( isempty ( ntk1 ) )

    ak12cs = [ ...
     +0.06379308343739001E+00, ...
     +0.02832887813049721E+00, ...
     -0.00024753706739052E+00, ...
     +0.00000577197245160E+00, ...
     -0.00000020689392195E+00, ...
     +0.00000000973998344E+00, ...
     -0.00000000055853361E+00, ...
     +0.00000000003732996E+00, ...
     -0.00000000000282505E+00, ...
     +0.00000000000023720E+00, ...
     -0.00000000000002176E+00, ...
     +0.00000000000000215E+00, ...
     -0.00000000000000022E+00, ...
     +0.00000000000000002E+00 ]';
    ak1cs = [ ...
     +0.2744313406973883E+00, ...
     +0.07571989953199368E+00, ...
     -0.00144105155647540E+00, ...
     +0.00006650116955125E+00, ...
     -0.00000436998470952E+00, ...
     +0.00000035402774997E+00, ...
     -0.00000003311163779E+00, ...
     +0.00000000344597758E+00, ...
     -0.00000000038989323E+00, ...
     +0.00000000004720819E+00, ...
     -0.00000000000604783E+00, ...
     +0.00000000000081284E+00, ...
     -0.00000000000011386E+00, ...
     +0.00000000000001654E+00, ...
     -0.00000000000000248E+00, ...
     +0.00000000000000038E+00, ...
     -0.00000000000000006E+00 ]';
    bk1cs = [ ...
     +0.0253002273389477705E+00, ...
     -0.353155960776544876E+00, ...
     -0.122611180822657148E+00, ...
     -0.0069757238596398643E+00, ...
     -0.0001730288957513052E+00, ...
     -0.0000024334061415659E+00, ...
     -0.0000000221338763073E+00, ...
     -0.0000000001411488392E+00, ...
     -0.0000000000006666901E+00, ...
     -0.0000000000000024274E+00, ...
     -0.0000000000000000070E+00 ]';

    ntk1 = r4_inits ( bk1cs, 11, 0.1 * r4_mach ( 3 ) );
    ntak1 = r4_inits ( ak1cs, 17, 0.1 * r4_mach ( 3 ) );
    ntak12 = r4_inits ( ak12cs, 14, 0.1 * r4_mach ( 3 ) );
    xmin = exp ( max ( log ( r4_mach ( 1 ) ), ...
      - log ( r4_mach ( 2 ) ) ) + 0.01 );
    xsml = sqrt ( 4.0 * r4_mach ( 3 ) );

  end

  if ( x <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_BESK1E = Fatal error!\n' );
    fprintf ( 1, '  X <= 0.\n' );
    error ( 'R4_BESK1E = Fatal error!' )
  elseif ( x <= xsml )
    y = 0.0;
    value = exp ( x ) * ( log ( 0.5 * x ) * r4_besi1 ( x ) ...
      + ( 0.75 + r4_csevl ( 0.5 * y - 1.0, bk1cs, ntk1 ) ) / x );
  elseif ( x <= 2.0 )
    y = x * x;
    value = exp ( x ) * ( log ( 0.5 * x ) * r4_besi1 ( x ) ...
      + ( 0.75 + r4_csevl ( 0.5 * y - 1.0, bk1cs, ntk1 ) ) / x );
  elseif ( x <= 8.0 )
    value = ( 1.25 + r4_csevl ( ( 16.0 / x - 5.0 ) / 3.0, ...
      ak1cs, ntak1 ) ) / sqrt ( x );
  else
    value = ( 1.25 + r4_csevl ( 16. / x - 1.0, ak12cs, ntak12 ) ) ...
      / sqrt ( x );
  end

  return
end