function value = r4_aide ( x ) %*****************************************************************************80 % %% R4_AIDE: exponentially scaled derivative, Airy function Ai of an R4 argument. % % Discussion: % % if X <= 0, % R4_AIDE ( X ) = R4_AID ( X ) % else % R4_AIDE ( X ) = R4_AID ( X ) * exp ( 2/3 * X^(3/2) ) % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 26 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 exponentially scaled derivative of % the Airy function Ai of X. % persistent aifcs persistent aigcs persistent aip1cs persistent aip2cs persistent naif persistent naig persistent naip1 persistent naip2 persistent x2sml persistent x32sml persistent x3sml persistent xbig if ( isempty ( naif ) ) aifcs = [ ... 0.10527461226531408809E+00, ... 0.01183613628152997844E+00, ... 0.00012328104173225664E+00, ... 0.00000062261225638140E+00, ... 0.00000000185298887844E+00, ... 0.00000000000363328873E+00, ... 0.00000000000000504622E+00, ... 0.00000000000000000522E+00 ]'; aigcs = [ ... 0.021233878150918666852E+00, ... 0.086315930335214406752E+00, ... 0.001797594720383231358E+00, ... 0.000014265499875550693E+00, ... 0.000000059437995283683E+00, ... 0.000000000152403366479E+00, ... 0.000000000000264587660E+00, ... 0.000000000000000331562E+00, ... 0.000000000000000000314E+00 ]'; aip1cs = [ ... 0.0358865097808301538E+00, ... 0.0114668575627764899E+00, ... -0.0007592073583861400E+00, ... 0.0000869517610893841E+00, ... -0.0000128237294298592E+00, ... 0.0000022062695681038E+00, ... -0.0000004222295185921E+00, ... 0.0000000874686415726E+00, ... -0.0000000192773588418E+00, ... 0.0000000044668460054E+00, ... -0.0000000010790108052E+00, ... 0.0000000002700029447E+00, ... -0.0000000000696480108E+00, ... 0.0000000000184489907E+00, ... -0.0000000000050027817E+00, ... 0.0000000000013852243E+00, ... -0.0000000000003908218E+00, ... 0.0000000000001121536E+00, ... -0.0000000000000326862E+00, ... 0.0000000000000096619E+00, ... -0.0000000000000028935E+00, ... 0.0000000000000008770E+00, ... -0.0000000000000002688E+00, ... 0.0000000000000000832E+00, ... -0.0000000000000000260E+00 ]'; aip2cs = [ ... 0.0065457691989713757E+00, ... 0.0023833724120774592E+00, ... -0.0000430700770220586E+00, ... 0.0000015629125858629E+00, ... -0.0000000815417186163E+00, ... 0.0000000054103738057E+00, ... -0.0000000004284130883E+00, ... 0.0000000000389497963E+00, ... -0.0000000000039623161E+00, ... 0.0000000000004428184E+00, ... -0.0000000000000536297E+00, ... 0.0000000000000069650E+00, ... -0.0000000000000009620E+00, ... 0.0000000000000001403E+00, ... -0.0000000000000000215E+00 ]'; eta = 0.1 * r4_mach ( 3 ); naif = r4_inits ( aifcs, 8, eta ); naig = r4_inits ( aigcs, 9, eta ); naip1 = r4_inits ( aip1cs, 25, eta ); naip2 = r4_inits ( aip2cs, 15, eta ); x2sml = sqrt ( eta ); x3sml = eta^0.3333; x32sml = 1.3104 * x3sml * x3sml; xbig = r4_mach ( 2 )^0.6666; end if ( x < - 1.0 ) [ xn, phi ] = r4_admp ( x ); value = xn * cos ( phi ); elseif ( abs ( x ) <= x2sml ) x2 = 0.0; x3 = 0.0; value = ( x2 * ( 0.125 ... + r4_csevl ( x3, aifcs, naif ) ) ... - r4_csevl ( x3, aigcs, naig ) ) - 0.25; elseif ( abs ( x ) <= x3sml ) x2 = x * x; x3 = 0.0; value = ( x2 * ( 0.125 ... + r4_csevl ( x3, aifcs, naif ) ) ... - r4_csevl ( x3, aigcs, naig ) ) - 0.25; elseif ( abs ( x ) <= x32sml ) x2 = x * x; x3 = x * x * x; value = ( x2 * ( 0.125 ... + r4_csevl ( x3, aifcs, naif ) ) ... - r4_csevl ( x3, aigcs, naig ) ) - 0.25; elseif ( x <= 1.0 ) x2 = x * x; x3 = x * x * x; value = ( x2 * ( 0.125 ... + r4_csevl ( x3, aifcs, naif ) ) ... - r4_csevl ( x3, aigcs, naig ) ) - 0.25; value = value * exp ( 2.0 * x * sqrt ( x ) / 3.0 ); elseif ( x <= 4.0 ) sqrtx = sqrt ( x ); z = ( 16.0 / ( x * sqrtx ) - 9.0 ) / 7.0; value = ( - 0.28125 - r4_csevl ( z, aip1cs, naip1 ) ) * sqrt ( sqrtx ); elseif ( x < xbig ) sqrtx = sqrt ( x ); z = 16.0 / ( x * sqrtx ) - 1.0; value = ( - 0.28125 - r4_csevl ( z, aip2cs, naip2 ) ) * sqrt ( sqrtx ); else sqrtx = sqrt ( x ); z = - 1.0; value = ( - 0.28125 - r4_csevl ( z, aip2cs, naip2 ) ) * sqrt ( sqrtx ); end return end