function value = r4_besi1e ( x )

%*****************************************************************************80
%
%% R4_BESI1E: exponentially scaled Bessel function I of order 1 of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 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 I
%    of order 1 of X.
%
  persistent ai1cs
  persistent ai12cs
  persistent bi1cs
  persistent ntai1
  persistent ntai12
  persistent nti1
  persistent xmin
  persistent xsml

  if ( isempty ( nti1 ) )

    ai1cs = [ ...
     -0.02846744181881479E+00, ...
     -0.01922953231443221E+00, ...
     -0.00061151858579437E+00, ...
     -0.00002069971253350E+00, ...
      0.00000858561914581E+00, ...
      0.00000104949824671E+00, ...
     -0.00000029183389184E+00, ...
     -0.00000001559378146E+00, ...
      0.00000001318012367E+00, ...
     -0.00000000144842341E+00, ...
     -0.00000000029085122E+00, ...
      0.00000000012663889E+00, ...
     -0.00000000001664947E+00, ...
     -0.00000000000166665E+00, ...
      0.00000000000124260E+00, ...
     -0.00000000000027315E+00, ...
      0.00000000000002023E+00, ...
      0.00000000000000730E+00, ...
     -0.00000000000000333E+00, ...
      0.00000000000000071E+00, ...
     -0.00000000000000006E+00 ]';
    ai12cs = [ ...
      0.02857623501828014E+00, ...
     -0.00976109749136147E+00, ...
     -0.00011058893876263E+00, ...
     -0.00000388256480887E+00, ...
     -0.00000025122362377E+00, ...
     -0.00000002631468847E+00, ...
     -0.00000000383538039E+00, ...
     -0.00000000055897433E+00, ...
     -0.00000000001897495E+00, ...
      0.00000000003252602E+00, ...
      0.00000000001412580E+00, ...
      0.00000000000203564E+00, ...
     -0.00000000000071985E+00, ...
     -0.00000000000040836E+00, ...
     -0.00000000000002101E+00, ...
      0.00000000000004273E+00, ...
      0.00000000000001041E+00, ...
     -0.00000000000000382E+00, ...
     -0.00000000000000186E+00, ...
      0.00000000000000033E+00, ...
      0.00000000000000028E+00, ...
     -0.00000000000000003E+00 ]';
    bi1cs = [ ...
     -0.001971713261099859E+00, ...
      0.40734887667546481E+00, ...
      0.034838994299959456E+00, ...
      0.001545394556300123E+00, ...
      0.000041888521098377E+00, ...
      0.000000764902676483E+00, ...
      0.000000010042493924E+00, ...
      0.000000000099322077E+00, ...
      0.000000000000766380E+00, ...
      0.000000000000004741E+00, ...
      0.000000000000000024E+00 ]';

    nti1 = r4_inits ( bi1cs, 11, 0.1 * r4_mach ( 3 ) );
    ntai1 = r4_inits ( ai1cs, 21, 0.1 * r4_mach ( 3 ) );
    ntai12 = r4_inits ( ai12cs, 22, 0.1 * r4_mach ( 3 ) );
    xmin = 2.0 * r4_mach ( 1 );
    xsml = sqrt ( 8.0 * r4_mach ( 3 ) );

  end

  y = abs ( x );

  if ( x == 0.0 )
    value = 0.0;
  elseif ( y <= xmin )
    value = 0.0;
  elseif ( y <= xsml )
    value = 0.5 * x;
    value = exp ( - y ) * value;
  elseif ( y <= 3.0 )
    value = x * ( 0.875 ...
      + r4_csevl ( y * y / 4.5 - 1.0, bi1cs, nti1 ) );
    value = exp ( - y ) * value;
  elseif ( y <= 8.0 )
    value = ( 0.375 + r4_csevl ( ( 48.0 / y - 11.0 ) / 5.0, ...
      ai1cs, ntai1) ) / sqrt ( y );
    if ( x < 0.0 )
      value = - value;
    end
  else
    value = ( 0.375 + r4_csevl ( 16.0 / y - 1.0, ai12cs, ntai12 ) ) ...
      / sqrt ( y );
    if ( x < 0.0 )
      value = - value;
    end
  end

  return
end