function value = r4_besj1 ( x )

%*****************************************************************************80
%
%% R4_BESJ1 evaluates the Bessel function J of order 1 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 1 of X.
%
  persistent bj1cs
  persistent bm1cs
  persistent bth1cs
  persistent ntj1
  persistent ntm1
  persistent ntth1
  persistent xmax
  persistent xmin
  persistent xsml

  if ( isempty ( ntj1 ) )

    bj1cs = [ ...
     -0.11726141513332787E+00, ...
     -0.25361521830790640E+00, ...
     +0.050127080984469569E+00, ...
     -0.004631514809625081E+00, ...
     +0.000247996229415914E+00, ...
     -0.000008678948686278E+00, ...
     +0.000000214293917143E+00, ...
     -0.000000003936093079E+00, ...
     +0.000000000055911823E+00, ...
     -0.000000000000632761E+00, ...
     +0.000000000000005840E+00, ...
     -0.000000000000000044E+00 ]';
    bm1cs = [ ...
     +0.1047362510931285E+00, ...
     +0.00442443893702345E+00, ...
     -0.00005661639504035E+00, ...
     +0.00000231349417339E+00, ...
     -0.00000017377182007E+00, ...
     +0.00000001893209930E+00, ...
     -0.00000000265416023E+00, ...
     +0.00000000044740209E+00, ...
     -0.00000000008691795E+00, ...
     +0.00000000001891492E+00, ...
     -0.00000000000451884E+00, ...
     +0.00000000000116765E+00, ...
     -0.00000000000032265E+00, ...
     +0.00000000000009450E+00, ...
     -0.00000000000002913E+00, ...
     +0.00000000000000939E+00, ...
     -0.00000000000000315E+00, ...
     +0.00000000000000109E+00, ...
     -0.00000000000000039E+00, ...
     +0.00000000000000014E+00, ...
     -0.00000000000000005E+00 ]';
    bth1cs = [ ...
     +0.74060141026313850E+00, ...
     -0.004571755659637690E+00, ...
     +0.000119818510964326E+00, ...
     -0.000006964561891648E+00, ...
     +0.000000655495621447E+00, ...
     -0.000000084066228945E+00, ...
     +0.000000013376886564E+00, ...
     -0.000000002499565654E+00, ...
     +0.000000000529495100E+00, ...
     -0.000000000124135944E+00, ...
     +0.000000000031656485E+00, ...
     -0.000000000008668640E+00, ...
     +0.000000000002523758E+00, ...
     -0.000000000000775085E+00, ...
     +0.000000000000249527E+00, ...
     -0.000000000000083773E+00, ...
     +0.000000000000029205E+00, ...
     -0.000000000000010534E+00, ...
     +0.000000000000003919E+00, ...
     -0.000000000000001500E+00, ...
     +0.000000000000000589E+00, ...
     -0.000000000000000237E+00, ...
     +0.000000000000000097E+00, ...
     -0.000000000000000040E+00 ]';

    ntj1 = r4_inits ( bj1cs, 12, 0.1 * r4_mach ( 3 ) );
    ntm1 = r4_inits ( bm1cs, 21, 0.1 * r4_mach ( 3 ) );
    ntth1 = r4_inits ( bth1cs, 24, 0.1 * r4_mach ( 3 ) );
    xsml = sqrt ( 8.0 * r4_mach ( 3 ) );
    xmin = 2.0 * r4_mach ( 1 );
    xmax = 1.0 / r4_mach ( 4 );

  end

  y = abs ( x );

  if ( y <= xmin )
    value = 0.0;
  elseif ( y <= xsml )
    value = 0.5 * x;
  elseif ( y <= 4.0 )
    value = x * ( 0.25 + r4_csevl ( 0.125 * y * y - 1.0, bj1cs, ntj1 ) );
  else
    z = 32.0 / y / y - 1.0;
    ampl = ( 0.75 + r4_csevl ( z, bm1cs, ntm1 ) ) / sqrt ( y );
    theta = y - 0.75 * pi + r4_csevl ( z, bth1cs, ntth1 ) / y;
    if ( x < 0.0 )
      value = - ampl * cos ( theta );
    else
      value = + ampl * cos ( theta );
    end
  end

  return
end