function value = bessel_i0 ( arg )

%*****************************************************************************80
%
%% BESSEL_I0 evaluates the modified Bessel function of the first kind and order zero.
%
%  Discussion:
%
%    The main computation evaluates slightly modified forms of
%    minimax approximations generated by Blair and Edwards, Chalk
%    River (Atomic Energy of Canada Limited) Report AECL-4928,
%    October, 1974.  This transportable program is patterned after
%    the machine dependent FUNPACK packet NATSI0, but cannot match
%    that version for efficiency or accuracy.  This version uses
%    rational functions that theoretically approximate I-SUB-0(X)
%    to at least 18 significant decimal digits.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 October 2004
%
%  Author:
%
%    Original FORTRAN77 version by William Cody, Laura Stoltz.
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, real ARG, the argument.
%
%    Output, real VALUE, the value of the modified Bessel function
%    of the first kind.
%
  exp40 = 2.353852668370199854E+17;

  p = [ ...
    -5.2487866627945699800E-18, ...
    -1.5982226675653184646E-14, ...
    -2.6843448573468483278E-11, ...
    -3.0517226450451067446E-08, ...
    -2.5172644670688975051E-05, ...
    -1.5453977791786851041E-02, ...
    -7.0935347449210549190E+00, ...
    -2.4125195876041896775E+03, ...
    -5.9545626019847898221E+05, ...
    -1.0313066708737980747E+08, ...
    -1.1912746104985237192E+10, ...
    -8.4925101247114157499E+11, ...
    -3.2940087627407749166E+13, ...
    -5.5050369673018427753E+14, ...
    -2.2335582639474375249E+15 ];

  pp = [ ...
    -3.9843750000000000000E-01, ...
     2.9205384596336793945E+00, ...
    -2.4708469169133954315E+00, ...
     4.7914889422856814203E-01, ...
    -3.7384991926068969150E-03, ...
    -2.6801520353328635310E-03, ...
     9.9168777670983678974E-05, ...
    -2.1877128189032726730E-06 ];

  q = [ ...
    -3.7277560179962773046E+03, ...
     6.5158506418655165707E+06, ...
    -6.5626560740833869295E+09, ...
     3.7604188704092954661E+12, ...
    -9.7087946179594019126E+14 ];

  qq = [ ...
    -3.1446690275135491500E+01, ...
     8.5539563258012929600E+01, ...
    -6.0228002066743340583E+01, ...
     1.3982595353892851542E+01, ...
    -1.1151759188741312645E+00, ...
     3.2547697594819615062E-02, ...
    -5.5194330231005480228E-04 ];

  rec15 = 6.6666666666666666666E-02;
  xmax = 91.9E+00;
  xsmall = 2.98E-08;

  x = abs ( arg );

  if ( x < xsmall )
    value = 1.0;
  elseif ( x < 15.0 )
%
%  XSMALL <= ABS(ARG) < 15.0
%
    xx = x * x;
    sump = p(1);
    for i = 2 : 15
      sump = sump * xx + p(i);
    end

    xx = xx - 225.0;
    sumq = (((( ...
           xx + q(1) ) ...
         * xx + q(2) ) ...
         * xx + q(3) ) ...
         * xx + q(4) ) ...
         * xx + q(5);

    value = sump / sumq;

  elseif ( 15.0 <= x )

    if ( xmax < x )
      value = r8_huge ( );
    else
%
%  15.0 <= ABS(ARG)
%
      xx = 1.0 / x - rec15;

      sump = ((((((  ...
                  pp(1) ...
           * xx + pp(2) ) ...
           * xx + pp(3) ) ...
           * xx + pp(4) ) ...
           * xx + pp(5) ) ...
           * xx + pp(6) ) ...
           * xx + pp(7) ) ...
           * xx + pp(8);

      sumq = (((((( ...
             xx + qq(1) ) ...
           * xx + qq(2) ) ...
           * xx + qq(3) ) ...
           * xx + qq(4) ) ...
           * xx + qq(5) ) ...
           * xx + qq(6) ) ...
           * xx + qq(7);

      value = sump / sumq;
%
%  Calculation reformulated to avoid premature overflow.
%
      if ( x <= xmax - 15.0 )
        a = exp ( x );
        b = 1.0;
      else
        a = exp ( x - 40.0 );
        b = exp40;
      end

      value = ( ( value * a - pp(1) * a ) / sqrt ( x ) ) * b;

    end

  end

  return
end