function value = r4_besj0 ( x )

%*****************************************************************************80
%
%% R4_BESJ0 evaluates the Bessel function J of order 0 of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 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 Bessel function J of order 0 of X.
%
  persistent bj0cs
  persistent bm0cs
  persistent bth0cs
  persistent ntj0
  persistent ntm0
  persistent ntth0
  persistent xmax
  persistent xsml

  if ( isempty ( ntj0 ) )

    bj0cs = [ ...
      0.100254161968939137E+00, ...
     -0.665223007764405132E+00, ...
      0.248983703498281314E+00, ...
     -0.0332527231700357697E+00, ...
      0.0023114179304694015E+00, ...
     -0.0000991127741995080E+00, ...
      0.0000028916708643998E+00, ...
     -0.0000000612108586630E+00, ...
      0.0000000009838650793E+00, ...
     -0.0000000000124235515E+00, ...
      0.0000000000001265433E+00, ...
     -0.0000000000000010619E+00, ...
      0.0000000000000000074E+00 ]';
    bm0cs = [ ...
      0.09284961637381644E+00, ...
     -0.00142987707403484E+00, ...
      0.00002830579271257E+00, ...
     -0.00000143300611424E+00, ...
      0.00000012028628046E+00, ...
     -0.00000001397113013E+00, ...
      0.00000000204076188E+00, ...
     -0.00000000035399669E+00, ...
      0.00000000007024759E+00, ...
     -0.00000000001554107E+00, ...
      0.00000000000376226E+00, ...
     -0.00000000000098282E+00, ...
      0.00000000000027408E+00, ...
     -0.00000000000008091E+00, ...
      0.00000000000002511E+00, ...
     -0.00000000000000814E+00, ...
      0.00000000000000275E+00, ...
     -0.00000000000000096E+00, ...
      0.00000000000000034E+00, ...
     -0.00000000000000012E+00, ...
      0.00000000000000004E+00 ]';
    bth0cs = [ ...
     -0.24639163774300119E+00, ...
      0.001737098307508963E+00, ...
     -0.000062183633402968E+00, ...
      0.000004368050165742E+00, ...
     -0.000000456093019869E+00, ...
      0.000000062197400101E+00, ...
     -0.000000010300442889E+00, ...
      0.000000001979526776E+00, ...
     -0.000000000428198396E+00, ...
      0.000000000102035840E+00, ...
     -0.000000000026363898E+00, ...
      0.000000000007297935E+00, ...
     -0.000000000002144188E+00, ...
      0.000000000000663693E+00, ...
     -0.000000000000215126E+00, ...
      0.000000000000072659E+00, ...
     -0.000000000000025465E+00, ...
      0.000000000000009229E+00, ...
     -0.000000000000003448E+00, ...
      0.000000000000001325E+00, ...
     -0.000000000000000522E+00, ...
      0.000000000000000210E+00, ...
     -0.000000000000000087E+00, ...
      0.000000000000000036E+00 ]';

    ntj0 = r4_inits ( bj0cs, 13, 0.1 * r4_mach ( 3 ) );
    ntm0 = r4_inits ( bm0cs, 21, 0.1 * r4_mach ( 3 ) );
    ntth0 = r4_inits ( bth0cs, 24, 0.1 * r4_mach ( 3 ) );
    xsml = sqrt ( 4.0 * r4_mach ( 3 ) );
    xmax = 1.0 / r4_mach ( 4 );

  end

  y = abs ( x );

  if ( y <= xsml )
    value = 1.0;
  elseif ( y <= 4.0 )
    value = r4_csevl ( 0.125 * y * y - 1.0, bj0cs, ntj0 );
  else
    z = 32.0 / y / y - 1.0;
    ampl = ( 0.75 + r4_csevl ( z, bm0cs, ntm0 ) ) / sqrt ( y );
    theta = y - 0.25 * pi + r4_csevl ( z, bth0cs, ntth0 ) / y;
    value = ampl * cos ( theta );
  end

  return
end
