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