function value = r4_besi0e ( x )

%*****************************************************************************80
%
%% R4_BESI0E evaluates the exponentially scaled Bessel function I0(X).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 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 I0(X).
%
  persistent ai02cs
  persistent ai0cs
  persistent bi0cs
  persistent ntai0
  persistent ntai02
  persistent nti0
  persistent xsml

  if ( isempty ( nti0 ) )
    ai02cs = [ ...
      0.05449041101410882E+00, ...
      0.00336911647825569E+00, ...
      0.00006889758346918E+00, ...
      0.00000289137052082E+00, ...
      0.00000020489185893E+00, ...
      0.00000002266668991E+00, ...
      0.00000000339623203E+00, ...
      0.00000000049406022E+00, ...
      0.00000000001188914E+00, ...
     -0.00000000003149915E+00, ...
     -0.00000000001321580E+00, ...
     -0.00000000000179419E+00, ...
      0.00000000000071801E+00, ...
      0.00000000000038529E+00, ...
      0.00000000000001539E+00, ...
     -0.00000000000004151E+00, ...
     -0.00000000000000954E+00, ...
      0.00000000000000382E+00, ...
      0.00000000000000176E+00, ...
     -0.00000000000000034E+00, ...
     -0.00000000000000027E+00, ...
      0.00000000000000003E+00 ]';
    ai0cs = [ ...
      0.07575994494023796E+00, ...
      0.00759138081082334E+00, ...
      0.00041531313389237E+00, ...
      0.00001070076463439E+00, ...
     -0.00000790117997921E+00, ...
     -0.00000078261435014E+00, ...
      0.00000027838499429E+00, ...
      0.00000000825247260E+00, ...
     -0.00000001204463945E+00, ...
      0.00000000155964859E+00, ...
      0.00000000022925563E+00, ...
     -0.00000000011916228E+00, ...
      0.00000000001757854E+00, ...
      0.00000000000112822E+00, ...
     -0.00000000000114684E+00, ...
      0.00000000000027155E+00, ...
     -0.00000000000002415E+00, ...
     -0.00000000000000608E+00, ...
      0.00000000000000314E+00, ...
     -0.00000000000000071E+00, ...
      0.00000000000000007E+00 ]';
    bi0cs = [ ...
     -0.07660547252839144951E+00, ...
      1.927337953993808270E+00, ...
      0.2282644586920301339E+00, ...
      0.01304891466707290428E+00, ...
      0.00043442709008164874E+00, ...
      0.00000942265768600193E+00, ...
      0.00000014340062895106E+00, ...
      0.00000000161384906966E+00, ...
      0.00000000001396650044E+00, ...
      0.00000000000009579451E+00, ...
      0.00000000000000053339E+00, ...
      0.00000000000000000245E+00 ]';
    nti0 = r4_inits ( bi0cs, 12, 0.1 * r4_mach ( 3 ) );
    ntai0 = r4_inits ( ai0cs, 21, 0.1 * r4_mach ( 3 ) );
    ntai02 = r4_inits ( ai02cs, 22, 0.1 * r4_mach ( 3 ) );
    xsml = sqrt ( 4.0 * r4_mach ( 3 ) );
  end

  y = abs ( x );

  if ( y <= xsml )
    value = 1.0;
  elseif ( y <= 3.0 )
    value = exp ( - y ) * ( 2.75 + ...
      r4_csevl ( y * y / 4.5 - 1.0, bi0cs, nti0 ) );
  elseif ( y <= 8.0 )
    value = ( 0.375 + r4_csevl ...
      ( ( 48.0 / y - 11.0 ) / 5.0, ai0cs, ntai0 ) ) / sqrt ( y );
  else
    value = ( 0.375 + r4_csevl ...
      ( 16.0 / y - 1.0, ai02cs, ntai02 ) ) / sqrt ( y );
  end

  return
end