function value = r4_besy1 ( x )

%*****************************************************************************80
%
%% R4_BESY1 evaluates the Bessel function Y of order 1 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 1 of X.
%
  persistent bm1cs
  persistent bth1cs
  persistent by1cs
  persistent ntm1
  persistent ntth1
  persistent nty1
  persistent pi4
  persistent twodpi
  persistent xmax
  persistent xmin
  persistent xsml

  pi4 = 0.78539816339744831;
  twodpi = 0.63661977236758134;

  if ( isempty ( nty1 ) )

    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 ]';
    by1cs = [ ...
     +0.03208047100611908629E+00, ...
     +1.262707897433500450E+00, ...
     +0.00649996189992317500E+00, ...
     -0.08936164528860504117E+00, ...
     +0.01325088122175709545E+00, ...
     -0.00089790591196483523E+00, ...
     +0.00003647361487958306E+00, ...
     -0.00000100137438166600E+00, ...
     +0.00000001994539657390E+00, ...
     -0.00000000030230656018E+00, ...
     +0.00000000000360987815E+00, ...
     -0.00000000000003487488E+00, ...
     +0.00000000000000027838E+00, ...
     -0.00000000000000000186E+00 ]';

    nty1 = r4_inits ( by1cs, 14, 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 ) );
    xmin = 1.571 * exp ( max ( log ( r4_mach ( 1 ) ), ...
      - log ( r4_mach ( 2 ) ) ) + 0.01 );
    xsml = sqrt ( 4.0 * r4_mach ( 3 ) );
    xmax = 1.0 / r4_mach ( 4 );

  end

  if ( x <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_BESY1 - Fatal error!\n' );
    fprintf ( 1, '  X <= 0.\n' );
    error ( 'R4_BESY1 - Fatal error!' )
  elseif ( x <= xsml )
    y = 0.0;
    value = twodpi * log ( 0.5 * x ) * r4_besj1 ( x ) ...
      + ( 0.5 + r4_csevl ( 0.125 * y - 1.0, by1cs, nty1 ) ) / x;
  elseif ( x <= 4.0 )
    y = x * x;
    value = twodpi * log ( 0.5 * x ) * r4_besj1 ( x ) ...
      + ( 0.5 + r4_csevl ( 0.125 * y - 1.0, by1cs, nty1 ) ) / x;
  else
    z = 32.0 / x / x - 1.0;
    ampl = ( 0.75 + r4_csevl ( z, bm1cs, ntm1 ) ) / sqrt ( x );
    theta = x - 3.0 * pi4 + r4_csevl ( z, bth1cs, ntth1 ) / x;
    value = ampl * sin ( theta );
  end

  return
end
