function value = r8_bi ( x )

%*****************************************************************************80
%
%% R8_BI evaluates the Airy function Bi of an R8 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 Airy function Bi of X.
%
  persistent bifcs
  persistent bif2cs
  persistent bigcs
  persistent big2cs
  persistent nbif
  persistent nbif2
  persistent nbig
  persistent nbig2
  persistent x3sml
  persistent xmax

  if ( isempty ( nbif ) )

  bifcs = [ ...
    -0.16730216471986649483537423928176E-01, ...
    +0.10252335834249445611426362777757, ...
    +0.17083092507381516539429650242013E-02, ...
    +0.11862545467744681179216459210040E-04, ...
    +0.44932907017792133694531887927242E-07, ...
    +0.10698207143387889067567767663628E-09, ...
    +0.17480643399771824706010517628573E-12, ...
    +0.20810231071761711025881891834399E-15, ...
    +0.18849814695665416509927971733333E-18, ...
    +0.13425779173097804625882666666666E-21, ...
    +0.77159593429658887893333333333333E-25, ...
    +0.36533879617478566399999999999999E-28, ...
    +0.14497565927953066666666666666666E-31 ]';
  bif2cs = [ ...
    +0.0998457269381604104468284257993, ...
    +0.47862497786300553772211467318231, ...
    +0.25155211960433011771324415436675E-01, ...
    +0.58206938852326456396515697872216E-03, ...
    +0.74997659644377865943861457378217E-05, ...
    +0.61346028703493836681403010356474E-07, ...
    +0.34627538851480632900434268733359E-09, ...
    +0.14288910080270254287770846748931E-11, ...
    +0.44962704298334641895056472179200E-14, ...
    +0.11142323065833011708428300106666E-16, ...
    +0.22304791066175002081517866666666E-19, ...
    +0.36815778736393142842922666666666E-22, ...
    +0.50960868449338261333333333333333E-25, ...
    +0.60003386926288554666666666666666E-28, ...
    +0.60827497446570666666666666666666E-31 ]';
  bigcs = [ ...
    +0.22466223248574522283468220139024E-01, ...
    +0.37364775453019545441727561666752E-01, ...
    +0.44476218957212285696215294326639E-03, ...
    +0.24708075636329384245494591948882E-05, ...
    +0.79191353395149635134862426285596E-08, ...
    +0.16498079851827779880887872402706E-10, ...
    +0.24119906664835455909247501122841E-13, ...
    +0.26103736236091436985184781269333E-16, ...
    +0.21753082977160323853123792000000E-19, ...
    +0.14386946400390433219483733333333E-22, ...
    +0.77349125612083468629333333333333E-26, ...
    +0.34469292033849002666666666666666E-29, ...
    +0.12938919273216000000000000000000E-32 ]';
  big2cs = [ ...
    +0.033305662145514340465176188111647, ...
    +0.161309215123197067613287532084943, ...
    +0.631900730961342869121615634921173E-02, ...
    +0.118790456816251736389780192304567E-03, ...
    +0.130453458862002656147116485012843E-05, ...
    +0.937412599553521729546809615508936E-08, ...
    +0.474580188674725153788510169834595E-10, ...
    +0.178310726509481399800065667560946E-12, ...
    +0.516759192784958180374276356640000E-15, ...
    +0.119004508386827125129496251733333E-17, ...
    +0.222982880666403517277063466666666E-20, ...
    +0.346551923027689419722666666666666E-23, ...
    +0.453926336320504514133333333333333E-26, ...
    +0.507884996513522346666666666666666E-29, ...
    +0.491020674696533333333333333333333E-32 ]';

    eta = 0.1 * r8_mach ( 3 );
    nbif = r8_inits ( bifcs, 13, eta );
    nbig = r8_inits ( bigcs, 13, eta );
    nbif2 = r8_inits ( bif2cs, 15, eta );
    nbig2 = r8_inits ( big2cs, 15, eta );
    x3sml = eta^0.3333;
    xmax = ( 1.5 * log ( r8_mach ( 2 ) ) )^0.6666;

  end

  if ( x < - 1.0 )
    [ xm, theta ] = r8_aimp ( x );
    value = xm * sin ( theta );
  elseif ( abs ( x ) <= x3sml )
    z = 0.0;
    value = 0.625 + r8_csevl ( z, bifcs, nbif ) ...
      + x * ( 0.4375 + r8_csevl ( z, bigcs, nbig ) );
  elseif ( x <= 1.0 )
    z = x * x * x;
    value = 0.625 + r8_csevl ( z, bifcs, nbif ) ...
      + x * ( 0.4375 + r8_csevl ( z, bigcs, nbig ) );
  elseif ( x <= 2.0 )
    z = ( 2.0 * x * x * x - 9.0 ) / 7.0;
    value = 1.125 + r8_csevl ( z, bif2cs, nbif2 ) ...
      + x * ( 0.625 + r8_csevl ( z, big2cs, nbig2 ) );
  else
    value = r8_bie ( x ) * exp ( 2.0 * x * sqrt ( x ) / 3.0 );
  end

  return
end
