function value = r4_besy0 ( x )

%*****************************************************************************80
%
%% R4_BESY0 evaluates the Bessel function Y of order 0 of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 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 Y of order 0 of X.
%
  persistent alnhaf
  persistent bm0cs
  persistent bth0cs
  persistent by0cs
  persistent ntm0
  persistent ntth0
  persistent nty0
  persistent pi4
  persistent twodpi
  persistent xmax
  persistent xsml

  twodpi = 0.63661977236758134;
  pi4 = 0.78539816339744831;
  alnhaf = -0.693147180559945309;

  if ( isempty ( nty0 ) )

    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 ]';
    by0cs = [ ...
     -0.011277839392865573E+00, ...
     -0.12834523756042035E+00, ...
     -0.10437884799794249E+00, ...
     +0.023662749183969695E+00, ...
     -0.002090391647700486E+00, ...
     +0.000103975453939057E+00, ...
     -0.000003369747162423E+00, ...
     +0.000000077293842676E+00, ...
     -0.000000001324976772E+00, ...
     +0.000000000017648232E+00, ...
     -0.000000000000188105E+00, ...
     +0.000000000000001641E+00, ...
     -0.000000000000000011E+00 ]';

    nty0 = r4_inits ( by0cs, 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

  if ( x <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_BESY0 - Fatal error!\n' );
    fprintf ( 1, '  X <= 0.\n' );
    error ( 'R4_BESY0 - Fatal error!' )
  elseif ( x <= xsml )
    y = 0.0;
    value = twodpi * ( alnhaf + log ( x ) ) ...
      * r4_besj0 ( x ) + 0.375 ...
      + r4_csevl ( 0.125 * y - 1.0, by0cs, nty0 );
  elseif ( x <= 4.0 )
    y = x * x;
    value = twodpi * ( alnhaf + log ( x ) ) ...
      * r4_besj0 ( x ) + 0.375 ...
      + r4_csevl ( 0.125 * y - 1.0, by0cs, nty0 );
  else
    z = 32.0 / x / x - 1.0;
    ampl = ( 0.75 + r4_csevl ( z, bm0cs, ntm0 ) ) / sqrt ( x );
    theta = x - pi4 + r4_csevl ( z, bth0cs, ntth0 ) / x;
    value = ampl * sin ( theta );
  end

  return
end