function value = r4_bie ( x )

%*****************************************************************************80
%
%% R4_BIE evaluates the exponentially scaled Airy function Bi of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 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 exponentially scaled Airy
%    function Bi of X.
%
  persistent bif2cs
  persistent bifcs
  persistent big2cs
  persistent bigcs
  persistent bip2cs
  persistent bipcs
  persistent nbif
  persistent nbif2
  persistent nbig
  persistent nbig2
  persistent nbip
  persistent nbip2
  persistent x32sml
  persistent x3sml
  persistent xbig

  atr = 8.7506905708484345;
  btr = -2.093836321356054;

  if ( isempty ( nbif ) )

    bif2cs = [ ...
      0.09984572693816041E+00, ...
      0.478624977863005538E+00, ...
      0.0251552119604330118E+00, ...
      0.0005820693885232645E+00, ...
      0.0000074997659644377E+00, ...
      0.0000000613460287034E+00, ...
      0.0000000003462753885E+00, ...
      0.0000000000014288910E+00, ...
      0.0000000000000044962E+00, ...
      0.0000000000000000111E+00 ]';
    bifcs = [ ...
     -0.01673021647198664948E+00, ...
      0.1025233583424944561E+00, ...
      0.00170830925073815165E+00, ...
      0.00001186254546774468E+00, ...
      0.00000004493290701779E+00, ...
      0.00000000010698207143E+00, ...
      0.00000000000017480643E+00, ...
      0.00000000000000020810E+00, ...
      0.00000000000000000018E+00 ]';
    big2cs = [ ...
      0.033305662145514340E+00, ...
      0.161309215123197068E+00, ...
      0.0063190073096134286E+00, ...
      0.0001187904568162517E+00, ...
      0.0000013045345886200E+00, ...
      0.0000000093741259955E+00, ...
      0.0000000000474580188E+00, ...
      0.0000000000001783107E+00, ...
      0.0000000000000005167E+00, ...
      0.0000000000000000011E+00 ]';
    bigcs = [ ...
      0.02246622324857452E+00, ...
      0.03736477545301955E+00, ...
      0.00044476218957212E+00, ...
      0.00000247080756363E+00, ...
      0.00000000791913533E+00, ...
      0.00000000001649807E+00, ...
      0.00000000000002411E+00, ...
      0.00000000000000002E+00 ]';
    bip2cs = [ ...
     -0.113596737585988679E+00, ...
      0.0041381473947881595E+00, ...
      0.0001353470622119332E+00, ...
      0.0000104273166530153E+00, ...
      0.0000013474954767849E+00, ...
      0.0000001696537405438E+00, ...
     -0.0000000100965008656E+00, ...
     -0.0000000167291194937E+00, ...
     -0.0000000045815364485E+00, ...
      0.0000000003736681366E+00, ...
      0.0000000005766930320E+00, ...
      0.0000000000621812650E+00, ...
     -0.0000000000632941202E+00, ...
     -0.0000000000149150479E+00, ...
      0.0000000000078896213E+00, ...
      0.0000000000024960513E+00, ...
     -0.0000000000012130075E+00, ...
     -0.0000000000003740493E+00, ...
      0.0000000000002237727E+00, ...
      0.0000000000000474902E+00, ...
     -0.0000000000000452616E+00, ...
     -0.0000000000000030172E+00, ...
      0.0000000000000091058E+00, ...
     -0.0000000000000009814E+00, ...
     -0.0000000000000016429E+00, ...
      0.0000000000000005533E+00, ...
      0.0000000000000002175E+00, ...
     -0.0000000000000001737E+00, ...
     -0.0000000000000000010E+00 ]';
    bipcs = [ ...
     -0.08322047477943447E+00, ...
      0.01146118927371174E+00, ...
      0.00042896440718911E+00, ...
     -0.00014906639379950E+00, ...
     -0.00001307659726787E+00, ...
      0.00000632759839610E+00, ...
     -0.00000042226696982E+00, ...
     -0.00000019147186298E+00, ...
      0.00000006453106284E+00, ...
     -0.00000000784485467E+00, ...
     -0.00000000096077216E+00, ...
      0.00000000070004713E+00, ...
     -0.00000000017731789E+00, ...
      0.00000000002272089E+00, ...
      0.00000000000165404E+00, ...
     -0.00000000000185171E+00, ...
      0.00000000000059576E+00, ...
     -0.00000000000012194E+00, ...
      0.00000000000001334E+00, ...
      0.00000000000000172E+00, ...
     -0.00000000000000145E+00, ...
      0.00000000000000049E+00, ...
     -0.00000000000000011E+00, ...
      0.00000000000000001E+00 ]';

    eta = 0.1 * r4_mach ( 3 );
    nbif = r4_inits ( bifcs, 9, eta );
    nbig = r4_inits ( bigcs, 8, eta );
    nbif2 = r4_inits ( bif2cs, 10, eta );
    nbig2 = r4_inits ( big2cs, 10, eta );
    nbip  = r4_inits ( bipcs, 24, eta );
    nbip2 = r4_inits ( bip2cs, 29, eta );
    x3sml = eta^0.3333;
    x32sml = 1.3104 * x3sml * x3sml;
    xbig = r4_mach ( 2 )^0.6666;

  end

  if ( x < -1.0 )
    [ xm, theta ] = r4_aimp ( x );
    value = xm * sin ( theta );
  elseif ( abs ( x ) <= x32sml )
    z = 0.0;
    value = 0.625 + r4_csevl ( z, bifcs, nbif ) ...
      + x * ( 0.4375 + r4_csevl ( z, bigcs, nbig ) );
  elseif ( abs ( x ) <= x3sml )
    z = 0.0;
    value = 0.625 + r4_csevl ( z, bifcs, nbif ) ...
      + x * ( 0.4375 + r4_csevl ( z, bigcs, nbig ) );
    value = value * exp ( - 2.0 * x * sqrt ( x ) / 3.0 );
  elseif ( x <= 1.0 )
    z = x * x * x;
    value = 0.625 + r4_csevl ( z, bifcs, nbif ) ...
      + x * ( 0.4375 + r4_csevl ( z, bigcs, nbig ) );
    value = value * exp ( - 2.0 * x * sqrt ( x ) / 3.0 );
  elseif ( x <= 2.0 )
    z = ( 2.0 * x * x * x - 9.0 ) / 7.0;
    value = exp ( - 2.0 * x * sqrt ( x ) / 3.0 ) ...
      * ( 1.125 + r4_csevl ( z, bif2cs, nbif2 )...
      + x * ( 0.625 + r4_csevl ( z, big2cs, nbig2 ) ) );
  elseif ( x <= 4.0 )
    sqrtx = sqrt ( x );
    z = atr / ( x * sqrtx ) + btr;
    value = ( 0.625 + r4_csevl ( z, bipcs, nbip ) ) / sqrt ( sqrtx );
  elseif ( x <= xbig )
    sqrtx = sqrt ( x );
    z = 16.0 / ( x * sqrtx ) - 1.0;
    value = ( 0.625 + r4_csevl ( z, bip2cs, nbip2 ) ) / sqrt ( sqrtx );
  else
    sqrtx = sqrt ( x );
    z = - 1.0;
    value = ( 0.625 + r4_csevl ( z, bip2cs, nbip2 ) ) / sqrt ( sqrtx );
  end

  return
end
