function [ ampl, theta ] = r4_aimp ( x )

%*****************************************************************************80
%
%% R4_AIMP evaluates the modulus and phase of the Airy function.
%
%  Description:
%
%    This function must only be called when X <= -1.0.
%
%  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 AMPL, PHI, the modulus and phase of the Airy function.
%
  persistent am21cs
  persistent am22cs
  persistent ath1cs
  persistent ath2cs
  persistent nam21
  persistent nam22
  persistent nath1
  persistent nath2
  persistent xsml

  if ( isempty ( nam21 ) )

    am21cs = [ ...
     +0.0065809191761485E+00, ...
     +0.0023675984685722E+00, ...
     +0.0001324741670371E+00, ...
     +0.0000157600904043E+00, ...
     +0.0000027529702663E+00, ...
     +0.0000006102679017E+00, ...
     +0.0000001595088468E+00, ...
     +0.0000000471033947E+00, ...
     +0.0000000152933871E+00, ...
     +0.0000000053590722E+00, ...
     +0.0000000020000910E+00, ...
     +0.0000000007872292E+00, ...
     +0.0000000003243103E+00, ...
     +0.0000000001390106E+00, ...
     +0.0000000000617011E+00, ...
     +0.0000000000282491E+00, ...
     +0.0000000000132979E+00, ...
     +0.0000000000064188E+00, ...
     +0.0000000000031697E+00, ...
     +0.0000000000015981E+00, ...
     +0.0000000000008213E+00, ...
     +0.0000000000004296E+00, ...
     +0.0000000000002284E+00, ...
     +0.0000000000001232E+00, ...
     +0.0000000000000675E+00, ...
     +0.0000000000000374E+00, ...
     +0.0000000000000210E+00, ...
     +0.0000000000000119E+00, ...
     +0.0000000000000068E+00, ...
     +0.0000000000000039E+00, ...
     +0.0000000000000023E+00, ...
     +0.0000000000000013E+00, ...
     +0.0000000000000008E+00, ...
     +0.0000000000000005E+00, ...
     +0.0000000000000003E+00, ...
     +0.0000000000000001E+00, ...
     +0.0000000000000001E+00, ...
     +0.0000000000000000E+00, ...
     +0.0000000000000000E+00, ...
     +0.0000000000000000E+00 ]';
    am22cs = [ ...
     -0.01562844480625341E+00, ...
     +0.00778336445239681E+00, ...
     +0.00086705777047718E+00, ...
     +0.00015696627315611E+00, ...
     +0.00003563962571432E+00, ...
     +0.00000924598335425E+00, ...
     +0.00000262110161850E+00, ...
     +0.00000079188221651E+00, ...
     +0.00000025104152792E+00, ...
     +0.00000008265223206E+00, ...
     +0.00000002805711662E+00, ...
     +0.00000000976821090E+00, ...
     +0.00000000347407923E+00, ...
     +0.00000000125828132E+00, ...
     +0.00000000046298826E+00, ...
     +0.00000000017272825E+00, ...
     +0.00000000006523192E+00, ...
     +0.00000000002490471E+00, ...
     +0.00000000000960156E+00, ...
     +0.00000000000373448E+00, ...
     +0.00000000000146417E+00, ...
     +0.00000000000057826E+00, ...
     +0.00000000000022991E+00, ...
     +0.00000000000009197E+00, ...
     +0.00000000000003700E+00, ...
     +0.00000000000001496E+00, ...
     +0.00000000000000608E+00, ...
     +0.00000000000000248E+00, ...
     +0.00000000000000101E+00, ...
     +0.00000000000000041E+00, ...
     +0.00000000000000017E+00, ...
     +0.00000000000000007E+00, ...
     +0.00000000000000002E+00 ]';
    ath1cs = [ ...
     -0.07125837815669365E+00, ...
     -0.00590471979831451E+00, ...
     -0.00012114544069499E+00, ...
     -0.00000988608542270E+00, ...
     -0.00000138084097352E+00, ...
     -0.00000026142640172E+00, ...
     -0.00000006050432589E+00, ...
     -0.00000001618436223E+00, ...
     -0.00000000483464911E+00, ...
     -0.00000000157655272E+00, ...
     -0.00000000055231518E+00, ...
     -0.00000000020545441E+00, ...
     -0.00000000008043412E+00, ...
     -0.00000000003291252E+00, ...
     -0.00000000001399875E+00, ...
     -0.00000000000616151E+00, ...
     -0.00000000000279614E+00, ...
     -0.00000000000130428E+00, ...
     -0.00000000000062373E+00, ...
     -0.00000000000030512E+00, ...
     -0.00000000000015239E+00, ...
     -0.00000000000007758E+00, ...
     -0.00000000000004020E+00, ...
     -0.00000000000002117E+00, ...
     -0.00000000000001132E+00, ...
     -0.00000000000000614E+00, ...
     -0.00000000000000337E+00, ...
     -0.00000000000000188E+00, ...
     -0.00000000000000105E+00, ...
     -0.00000000000000060E+00, ...
     -0.00000000000000034E+00, ...
     -0.00000000000000020E+00, ...
     -0.00000000000000011E+00, ...
     -0.00000000000000007E+00, ...
     -0.00000000000000004E+00, ...
     -0.00000000000000002E+00 ]';
    ath2cs = [ ...
     +0.00440527345871877E+00, ...
     -0.03042919452318455E+00, ...
     -0.00138565328377179E+00, ...
     -0.00018044439089549E+00, ...
     -0.00003380847108327E+00, ...
     -0.00000767818353522E+00, ...
     -0.00000196783944371E+00, ...
     -0.00000054837271158E+00, ...
     -0.00000016254615505E+00, ...
     -0.00000005053049981E+00, ...
     -0.00000001631580701E+00, ...
     -0.00000000543420411E+00, ...
     -0.00000000185739855E+00, ...
     -0.00000000064895120E+00, ...
     -0.00000000023105948E+00, ...
     -0.00000000008363282E+00, ...
     -0.00000000003071196E+00, ...
     -0.00000000001142367E+00, ...
     -0.00000000000429811E+00, ...
     -0.00000000000163389E+00, ...
     -0.00000000000062693E+00, ...
     -0.00000000000024260E+00, ...
     -0.00000000000009461E+00, ...
     -0.00000000000003716E+00, ...
     -0.00000000000001469E+00, ...
     -0.00000000000000584E+00, ...
     -0.00000000000000233E+00, ...
     -0.00000000000000093E+00, ...
     -0.00000000000000037E+00, ...
     -0.00000000000000015E+00, ...
     -0.00000000000000006E+00, ...
     -0.00000000000000002E+00 ]';

    eta = 0.1 * r4_mach ( 3 );
    nam21 = r4_inits ( am21cs, 40, eta );
    nath1 = r4_inits ( ath1cs, 36, eta );
    nam22 = r4_inits ( am22cs, 33, eta );
    nath2 = r4_inits ( ath2cs, 32, eta );
    xsml = - ( 16.0 / r4_mach ( 3 ) )^0.3333;

  end

  if ( x <= xsml )
    z = 1.0;
    ampl = 0.3125 + r4_csevl ( z, am21cs, nam21 );
    theta = - 0.625 + r4_csevl ( z, ath1cs, nath1 );
  elseif ( x < - 2.0 )
    z = 16.0 / x / x / x + 1.0;
    ampl = 0.3125 + r4_csevl ( z, am21cs, nam21 );
    theta = - 0.625 + r4_csevl ( z, ath1cs, nath1 );
  elseif ( x <= - 1.0 )
    z = ( 16.0 / x / x / x + 9.0 ) / 7.0;
    ampl = 0.3125 + r4_csevl ( z, am22cs, nam22 );
    theta = - 0.625 + r4_csevl ( z, ath2cs, nath2 );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_AIMP - Fatal error!\n' );
    fprintf ( 1, '  - 1.0 < X.\n' );
    error ( 'R4_AIMP - Fatal error!' )
  end

  sqrtx = sqrt ( - x );
  ampl = sqrt ( ampl / sqrtx );
  theta = 0.25 * pi - x * sqrtx * theta;

  return
end
