function value = 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:
%
%    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.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 = [ ...
      +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 = [ ...
      +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 = exp ( max ( log ( r8_mach ( 1 ) ), ...
      - log ( r8_mach ( 2 ) ) ) + 0.01 );
    xsml = sqrt ( 4.0 * r8_mach ( 3 ) );

  end

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

  return
end