function main(varargin)
%*******************************************************************************
%
%! TOMS757_PRB tests the routines in TOMS757.
%
%  Discussion:
%
%    This program tests the 37 functions in the MISCFUN package.
%    It is a fairly simple code with each function being tested
%    at 20 different arguments.  The code compares the value
%    from the function with a pre-computed value, and produces
%    the absolute and relative errors.
%
%  Author:
%
%    Allan McLeod,
%    Department of Mathematics and Statistics,
%    Paisley University, High Street, Paisley, Scotland, PA12BE
%    macl_ms0@paisley.ac.uk
%
clear global; clear functions;
global GlobInArgs nargs
GlobInArgs={mfilename,varargin{:}}; nargs=nargin+1;

timestamp( );
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TOMS757_PRB:');
writef(1,['%s','\n'], '  Test the uncommon special function routines in TOMS757.');
test01;
test02;
test03;
test04;
test05;
test06;
test07;
test08;
test09;
test10;
test11;
test12;
test13;
test14;
test15;
test16;
test17;
test18;
test19;
test20;
test21;
test22;
test23;
test24;
test25;
test26;
test27;
test28;
test29;
test30;
test31;
test32;
test33;
test34;
test35;
test36;
test37;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TOMS757_PRB:');
writef(1,['%s','\n'], '  Normal end of execution.');
writef(1,['%s','\n'], ' ');
timestamp( );
%stop
end
function test01(varargin)
%*******************************************************************************
%
%! TEST01 tests ABRAM0.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST01');
writef(1,['%s','\n'], '  Testing function ABRAM0');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=abram0_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=abram0( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test02(varargin)
%*******************************************************************************
%
%! TEST02 tests ABRAM1.
%
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST02');
writef(1,['%s','\n'], '  Testing function ABRAM1');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=abram1_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=abram1( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test03(varargin)
%*******************************************************************************
%
%! TEST03 tests ABRAM2.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST03');
writef(1,['%s','\n'], '  Testing function ABRAM2');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=abram2_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=abram2( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test04(varargin)
%*******************************************************************************
%
%! TEST04 tests AIRY_AI_INT.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST04');
writef(1,['%s','\n'], '  Testing function AIRY_AI_INT');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=airy_ai_int_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=airy_ai_int( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test05(varargin)
%*******************************************************************************
%
%! TEST05 tests AIRY_BI_INT.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST05');
writef(1,['%s','\n'], '  Testing function AIRY_BI_INT');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=airy_bi_int_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=airy_bi_int( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test06(varargin)
%*******************************************************************************
%
%! TEST06 tests AIRY_GI.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST06');
writef(1,['%s','\n'], '  Testing function AIRY_GI');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=airy_gi_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=airy_gi( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test07(varargin)
%*******************************************************************************
%
%! TEST07 tests AIRY_HI.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST07');
writef(1,['%s','\n'], '  Testing function AIRY_HI');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=airy_hi_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=airy_hi( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test08(varargin)
%*******************************************************************************
%
%! TEST08 tests ARCTAN_INT.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST08');
writef(1,['%s','\n'], '  Testing function ARCTAN_INT');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=arctan_int_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=arctan_int( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test09(varargin)
%*******************************************************************************
%
%! TEST09 tests BESSEL_I0_INT.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST09');
writef(1,['%s','\n'], '  Testing function BESSEL_I0_INT');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=bessel_i0_int_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=bessel_i0_int( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test10(varargin)
%*******************************************************************************
%
%! TEST10 tests BESSEL_J0_INT.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST10');
writef(1,['%s','\n'], '  Testing function BESSEL_J0_INT');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=bessel_j0_int_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=bessel_j0_int( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test11(varargin)
%*******************************************************************************
%
%! TEST11 tests BESSEL_K0_INT.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST11');
writef(1,['%s','\n'], '  Testing function BESSEL_K0_INT');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=bessel_k0_int_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=bessel_k0_int( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test12(varargin)
%*******************************************************************************
%
%! TEST12 tests BESSEL_Y0_INT.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST12');
writef(1,['%s','\n'], '  Testing function BESSEL_Y0_INT');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=bessel_y0_int_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=bessel_y0_int( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test13(varargin)
%*******************************************************************************
%
%! TEST13 tests CLAUSEN.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST13');
writef(1,['%s','\n'], '  Testing function CLAUSEN');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=clausen_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=clausen( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test14(varargin)
%*******************************************************************************
%
%! TEST14 tests DEBYE1.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST14');
writef(1,['%s','\n'], '  Testing function DEBYE1');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=debye1_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=debye1( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test15(varargin)
%*******************************************************************************
%
%! TEST15 tests DEBYE2.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST15');
writef(1,['%s','\n'], '  Testing function DEBYE2');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=debye2_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=debye2( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test16(varargin)
%*******************************************************************************
%
%! TEST16 tests DEBYE3.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST16');
writef(1,['%s','\n'], '  Testing function DEBYE3');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=debye3_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=debye3( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test17(varargin)
%*******************************************************************************
%
%! TEST17 tests DEBYE4.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST17');
writef(1,['%s','\n'], '  Testing function DEBYE4');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=debye4_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=debye4( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test18(varargin)
%*******************************************************************************
%
%! TEST18 tests EXP3_INT.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST18');
writef(1,['%s','\n'], '  Testing function EXP3_INT');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=exp3_int_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=exp3_int( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test19(varargin)
%*******************************************************************************
%
%! TEST19 tests GOODWIN.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST19');
writef(1,['%s','\n'], '  Testing function GOODWIN');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=goodwin_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=goodwin( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test20(varargin)
%*******************************************************************************
%
%! TEST20 tests I0ML0.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST20');
writef(1,['%s','\n'], '  Testing function I0ML0');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=i0ml0_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=i0ml0( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test21(varargin)
%*******************************************************************************
%
%! TEST21 tests I1ML1.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST21');
writef(1,['%s','\n'], '  Testing function I1ML1');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=i1ml1_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=i1ml1( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test22(varargin)
%*******************************************************************************
%
%! TEST22 tests LOBACHEVSKY.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST22');
writef(1,['%s','\n'], '  Testing function LOBACHEVSKY');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=lobachevsky_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=lobachevsky( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test23(varargin)
%*******************************************************************************
%
%! TEST23 tests STROMGEN.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST23');
writef(1,['%s','\n'], '  Testing function STROMGEN');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=stromgen_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=stromgen( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test24(varargin)
%*******************************************************************************
%
%! TEST24 tests STRUVE_H0.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST24');
writef(1,['%s','\n'], '  Testing function STRUVE_H0');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=struve_h0_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=struve_h0( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test25(varargin)
%*******************************************************************************
%
%! TEST25 tests STRUVE_H1.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST25');
writef(1,['%s','\n'], '  Testing function STRUVE_H1');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=struve_h1_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=struve_h1( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test26(varargin)
%*******************************************************************************
%
%! TEST26 tests STRUVE_L0.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST26');
writef(1,['%s','\n'], '  Testing function STRUVE_L0');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=struve_l0_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=struve_l0( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test27(varargin)
%*******************************************************************************
%
%! TEST27 tests STRUVE_L1.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST27');
writef(1,['%s','\n'], '  Testing function STRUVE_L1');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=struve_l1_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=struve_l1( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test28(varargin)
%*******************************************************************************
%
%! TEST28 tests SYNCH1.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST28');
writef(1,['%s','\n'], '  Testing function SYNCH1');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=synch1_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=synch1( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test29(varargin)
%*******************************************************************************
%
%! TEST29 tests SYNCH2.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST29');
writef(1,['%s','\n'], '  Testing function SYNCH2');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=synch2_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=synch2( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test30(varargin)
%*******************************************************************************
%
%! TEST30 tests TRAN02.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST30');
writef(1,['%s','\n'], '  Testing function TRAN02');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=tran02_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=tran02( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test31(varargin)
%*******************************************************************************
%
%! TEST31 tests TRAN03.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST31');
writef(1,['%s','\n'], '  Testing function TRAN03');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=tran03_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=tran03( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test32(varargin)
%*******************************************************************************
%
%! TEST32 tests TRAN04.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST32');
writef(1,['%s','\n'], '  Testing function TRAN04');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=tran04_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=tran04( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test33(varargin)
%*******************************************************************************
%
%! TEST33 tests TRAN05.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST33');
writef(1,['%s','\n'], '  Testing function TRAN05');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=tran05_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=tran05( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test34(varargin)
%*******************************************************************************
%
%! TEST34 tests TRAN06.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST34');
writef(1,['%s','\n'], '  Testing function TRAN06');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=tran06_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=tran06( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test35(varargin)
%*******************************************************************************
%
%! TEST35 tests TRAN07.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST35');
writef(1,['%s','\n'], '  Testing function TRAN07');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=tran07_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=tran07( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test36(varargin)
%*******************************************************************************
%
%! TEST36 tests TRAN08.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST36');
writef(1,['%s','\n'], '  Testing function TRAN08');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=tran08_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=tran08( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function test37(varargin)
%*******************************************************************************
%
%! TEST37 tests TRAN09.
%
persistent abserr comp fx n_data relerr x 

if isempty(abserr), abserr=0; end;
if isempty(comp), comp=0; end;
if isempty(fx), fx=0; end;
if isempty(n_data), n_data=0; end;
if isempty(relerr), relerr=0; end;
if isempty(x), x=0; end;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'TEST37');
writef(1,['%s','\n'], '  Testing function TRAN09');
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], '     Argument             Abs. error             Rel. error');
writef(1,['%s','\n'], ' ');
n_data = 0;
while (1);
[ n_data, x, fx ]=tran09_values( n_data, x, fx );
if( n_data <= 0 ) ;
break;
end;
[comp , x ]=tran09( x );
abserr = abs( fx - comp );
relerr = abserr ./ abs( fx );
writef(1,[repmat(' ',1,2),'%15.10f',repmat(' ',1,2),'%15.5f',repmat(' ',1,8),'%15.5f','\n'], x, abserr, relerr);
end;
return;
end
function [abram0Result, xvalue ]=abram0( xvalue ,varargin);
%*******************************************************************************
%
%! ABRAM0 evaluates the Abramowitz function of order 0.
%
%  Discussion:
%
%    The function is defined by:
%
%      ABRAM0(x) = Integral ( 0 <= t < infinity ) exp ( -t^2 - x / t ) dt
%
%    The code uses Chebyshev expansions with the coefficients
%    given to an accuracy of 20 decimal places.
%
%    This subroutine is set up to work on IEEE machines.
%
%  Modified:
%
%    07 August 2004
%
%  Author:
%
%    Allan McLeod,
%    Department of Mathematics and Statistics,
%    Paisley University, High Street, Paisley, Scotland, PA12BE
%    macl_ms0@paisley.ac.uk
%
%  Reference:
%
%    Allan Mcleod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input, real ( kind = 8 ) XVALUE, the argument of the function.
%
%    Output, real ( kind = 8 ) ABRAM0, the value of the function.
%
abram0Result=[];
persistent ab0as ab0f ab0g ab0h asln asval fval gval gval0 half hval lnxmin nterma ntermf ntermg ntermh onerpi rt3bpi rtpib2 six t three two v x xlow1 zero 

if isempty(ab0f), ab0f([0:8]+1) =[real( -0.68121927093549469816d0),real( -0.78867919816149252495d0),real(  0.5121581776818819543d-1),real( -0.71092352894541296d-3),real(  0.368681808504287d-5),real( -0.917832337237d-8),real(  0.1270202563d-10),real( -0.1076888d-13),real(  0.599d-17) ]; end;
if isempty(ab0g), ab0g([0:8]+1) =[real( -0.60506039430868273190d0),real( -0.41950398163201779803d0),real(  0.1703265125190370333d-1),real( -0.16938917842491397d-3),real(  0.67638089519710d-6),real( -0.135723636255d-8),real(  0.156297065d-11),real( -0.112887d-14),real(  0.55d-18) ]; end;
if isempty(ab0h), ab0h([0:8]+1) =[real(  1.38202655230574989705d0),real( -0.30097929073974904355d0),real(  0.794288809364887241d-2),real( -0.6431910276847563d-4),real(  0.22549830684374d-6),real( -0.41220966195d-9),real(  0.44185282d-12),real( -0.30123d-15),real(  0.14d-18)  ]; end;
if isempty(ab0as), ab0as([0:27]+1) =[real(  1.97755499723693067407d+0),real( -0.1046024792004819485d-1),real(  0.69680790253625366d-3),real( -0.5898298299996599d-4),real(  0.577164455305320d-5),real( -0.61523013365756d-6),real(  0.6785396884767d-7),real( -0.723062537907d-8),real(  0.63306627365d-9),real( -0.989453793d-11),real( -0.1681980530d-10),real(  0.673799551d-11),real( -0.200997939d-11),real(  0.54055903d-12),real( -0.13816679d-12),real(  0.3422205d-13),real( -0.826686d-14),real(  0.194566d-14),real( -0.44268d-15),real(  0.9562d-16),real( -0.1883d-16),real(  0.301d-17),real( -0.19d-18),real( -0.14d-18),real(  0.11d-18),real( -0.4d-19),real(  0.2d-19),real( -0.1d-19) ]; end;
if isempty(asln), asln=0; end;
if isempty(asval), asval=0; end;
if isempty(fval), fval=0; end;
if isempty(gval), gval=0; end;
if isempty(gval0), gval0 = 0.13417650264770070909d+00; end;
if isempty(half), half = 0.5d+00; end;
if isempty(hval), hval=0; end;
if isempty(lnxmin), lnxmin = -708.3964d+00; end;
if isempty(nterma), nterma = 22; end;
if isempty(ntermf), ntermf = 8; end;
if isempty(ntermg), ntermg = 8; end;
if isempty(ntermh), ntermh = 8; end;
if isempty(onerpi), onerpi = 0.56418958354775628695d+00; end;
if isempty(rt3bpi), rt3bpi = 0.97720502380583984317d+00; end;
if isempty(rtpib2), rtpib2 = 0.88622692545275801365d+00; end;
if isempty(six), six = 6.0d+00; end;
if isempty(t), t=0; end;
if isempty(three), three = 3.0d+00; end;
if isempty(two), two = 2.0d+00; end;
if isempty(v), v=0; end;
if isempty(x), x=0; end;
if isempty(xlow1), xlow1 = 1.490116d-08; end;
if isempty(zero), zero = 0.0d+00; end;
x = xvalue;
if( x < zero ) ;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'ABRAM0 - Fatal error!');
writef(1,['%s','\n'], '  Argument X < 0.');
abram0Result = zero;
elseif( x == zero ) ;
abram0Result = rtpib2;
elseif( x < xlow1 ) ;
abram0Result = rtpib2 + x .*( log( x ) - gval0 );
elseif( x <= two ) ;
t =( x .* x ./ two - half ) - half;
[fval , ntermf, ab0f, t ]=cheval( ntermf, ab0f, t );
[gval , ntermg, ab0g, t ]=cheval( ntermg, ab0g, t );
[hval , ntermh, ab0h, t ]=cheval( ntermh, ab0h, t );
abram0Result = fval ./ onerpi + x .*( log( x ) .* hval - gval );
else;
v = three .*(( x ./ two ) .^( two ./ three ) );
t =( six ./ v - half ) - half;
[asval , nterma, ab0as, t ]=cheval( nterma, ab0as, t );
asln = log( asval ./ rt3bpi ) - v;
if( asln < lnxmin ) ;
abram0Result = zero;
else;
abram0Result = exp( asln );
end;
end;
if ~isempty(inputname(1)), assignin('caller','FUntemp',xvalue); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end
function [n_data, x, fx]=abram0_values( n_data, x, fx ,varargin);
%*******************************************************************************
%
%! ABRAM0_VALUES returns some values of the Abramowitz0 function.
%
%  Discussion:
%
%    The function is defined by:
%
%      ABRAM0(x) = Integral ( 0 <= t < infinity ) exp ( -t^2 - x / t ) dt
%
%    The data was reported by McLeod.
%
%  Modified:
%
%    21 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real ( kind = 8 ) X, the argument of the function.
%
%    Output, real ( kind = 8 ) FX, the value of the function.
%
persistent fx_vec n_max x_vec 

if isempty(n_max), n_max = 20; end;
if isempty(fx_vec), fx_vec([1:n_max]) =[real( 0.87377726306985360531d+00),real( 0.84721859650456925922d+00),real( 0.77288934483988301615d+00),real( 0.59684345853450151603d+00),real( 0.29871735283675888392d+00),real( 0.15004596450516388138d+00),real( 0.11114662419157955096d+00),real( 0.83909567153151897766d-01),real( 0.56552321717943417515d-01),real( 0.49876496603033790206d-01),real( 0.44100889219762791328d-01),real( 0.19738535180254062496d-01),real( 0.86193088287161479900d-02),real( 0.40224788162540127227d-02),real( 0.19718658458164884826d-02),real( 0.10045868340133538505d-02),real( 0.15726917263304498649d-03),real( 0.10352666912350263437d-04),real( 0.91229759190956745069e-06),real( 0.25628287737952698742e-09) ]; end;
if isempty(x_vec), x_vec([1:n_max]) =[real( 0.0019531250d+00),real( 0.0078125000d+00),real( 0.0312500000d+00),real( 0.1250000000d+00),real( 0.5000000000d+00),real( 1.0000000000d+00),real( 1.2500000000d+00),real( 1.5000000000d+00),real( 1.8750000000d+00),real( 2.0000000000d+00),real( 2.1250000000d+00),real( 3.0000000000d+00),real( 4.0000000000d+00),real( 5.0000000000d+00),real( 6.0000000000d+00),real( 7.0000000000d+00),real( 10.0000000000d+00),real( 15.0000000000d+00),real( 20.0000000000d+00),real( 40.0000000000d+00) ]; end;
if( n_data < 0 ) ;
n_data = 0;
end;
n_data = n_data + 1;
if( n_max < n_data ) ;
n_data = 0;
x = 0.0d+00;
fx = 0.0d+00;
else;
x  = x_vec(n_data);
fx = fx_vec(n_data);
end;
return;
end
function [abram1Result, xvalue ]=abram1( xvalue ,varargin);
%*******************************************************************************
%
%! ABRAM1 evaluates the Abramowitz function of order 1.
%
%  Discussion:
%
%    The function is defined by:
%
%      ABRAM1(x) = Integral ( 0 <= t < infinity ) t * exp ( -t^2 - x / t ) dt
%
%    The code uses Chebyshev expansions with the coefficients
%    given to an accuracy of 20 decimal places.
%
%    This subroutine is set up to work on IEEE machines.
%
%  Modified:
%
%    07 August 2004
%
%  Author:
%
%    Allan McLeod,
%    Department of Mathematics and Statistics,
%    Paisley University, High Street, Paisley, Scotland, PA12BE
%    macl_ms0@paisley.ac.uk
%
%  Reference:
%
%    Allan Mcleod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input, real ( kind = 8 ) XVALUE, the argument of the function.
%
%    Output, real ( kind = 8 ) ABRAM1, the value of the function.
%
abram1Result=[];
persistent ab1as ab1f ab1g ab1h asln asval fval gval half hval lnxmin nterma ntermf ntermg ntermh one onerpi rt3bpi six t three two v x xlow xlow1 zero 

if isempty(asln), asln=0; end;
if isempty(asval), asval=0; end;
if isempty(fval), fval=0; end;
if isempty(gval), gval=0; end;
if isempty(half), half = 0.5d+00; end;
if isempty(hval), hval=0; end;
if isempty(nterma), nterma = 23; end;
if isempty(ntermf), ntermf = 9; end;
if isempty(ntermg), ntermg = 8; end;
if isempty(ntermh), ntermh = 8; end;
if isempty(one), one = 1.0d+00; end;
if isempty(rt3bpi), rt3bpi = 0.97720502380583984317d+00; end;
if isempty(six), six = 6.0d+00; end;
if isempty(t), t=0; end;
if isempty(three), three = 3.0d+00; end;
if isempty(two), two = 2.0d+00; end;
if isempty(v), v=0; end;
if isempty(x), x=0; end;
if isempty(zero), zero = 0.0d+00; end;
if isempty(ab1f),   ab1f=[1.47285192577978807369d0,0.10903497570168956257d0,-0.12430675360056569753d0,0.306197946853493315d-2,-0.2218410323076511d-4,0.6989978834451d-7,-0.11597076444d-9,0.11389776d-12,-0.7173d-16,0.3d-19];  end;
if isempty(ab1g),   ab1g=[0.39791277949054503528d0,-0.29045285226454720849d0,0.1048784695465363504d-1,-0.10249869522691336d-3,0.41150279399110d-6,-0.83652638940d-9,0.97862595d-12,-0.71868d-15,0.35d-18];  end;
if isempty(ab1h),   ab1h=[0.84150292152274947030d0,-0.7790050698774143395d-1,0.133992455878390993d-2,-0.808503907152788d-5,0.2261858281728d-7,-0.3441395838d-10,0.3159858d-13,-0.1884d-16,0.1d-19];  end;
try;ab1as(0+1);catch;   ab1as(0+1)=[2.13013643429065549448d0];  end;
try;ab1as(1+1);catch;   ab1as(1+1)=[0.6371526795218539933d-1];  end;
try;ab1as(2+1);catch;   ab1as(2+1)=[-0.129334917477510647d-2];  end;
try;ab1as(3+1);catch;   ab1as(3+1)=[0.5678328753228265d-4];  end;
try;ab1as(4+1);catch;   ab1as(4+1)=[-0.279434939177646d-5];  end;
try;ab1as(5+1);catch;   ab1as(5+1)=[0.5600214736787d-7];  end;
try;ab1as(6+1);catch;   ab1as(6+1)=[0.2392009242798d-7];  end;
try;ab1as(7+1);catch;   ab1as(7+1)=[-0.750984865009d-8];  end;
try;ab1as(8+1);catch;   ab1as(8+1)=[0.173015330776d-8];  end;
try;ab1as(9+1);catch;   ab1as(9+1)=[-0.36648877955d-9];  end;
try;ab1as(10+1);catch;   ab1as(10+1)=[0.7520758307d-10];  end;
try;ab1as(11+1);catch;   ab1as(11+1)=[-0.1517990208d-10];  end;
try;ab1as(12+1);catch;   ab1as(12+1)=[0.301713710d-11];  end;
try;ab1as(13+1);catch;   ab1as(13+1)=[-0.58596718d-12];  end;
try;ab1as(14+1);catch;   ab1as(14+1)=[0.10914455d-12];  end;
try;ab1as(15+1);catch;   ab1as(15+1)=[-0.1870536d-13];  end;
try;ab1as(16+1);catch;   ab1as(16+1)=[0.262542d-14];  end;
try;ab1as(17+1);catch;   ab1as(17+1)=[-0.14627d-15];  end;
try;ab1as(18+1);catch;   ab1as(18+1)=[-0.9500d-16];  end;
try;ab1as(19+1);catch;   ab1as(19+1)=[0.5873d-16];  end;
try;ab1as(20+1);catch;   ab1as(20+1)=[-0.2420d-16];  end;
try;ab1as(21+1);catch;   ab1as(21+1)=[0.868d-17];  end;
try;ab1as(22+1);catch;   ab1as(22+1)=[-0.290d-17];  end;
try;ab1as(23+1);catch;   ab1as(23+1)=[0.93d-18];  end;
try;ab1as(24+1);catch;   ab1as(24+1)=[-0.29d-18];  end;
try;ab1as(25+1);catch;   ab1as(25+1)=[0.9d-19];  end;
try;ab1as(26+1);catch;   ab1as(26+1)=[-0.3d-19];  end;
try;ab1as(27+1);catch;   ab1as(27+1)=[0.1d-19];  end;
if isempty(onerpi),   onerpi=[0.56418958354775628695d0];  end;
%
%  Machine-dependent constants (suitable for IEEE machines)
%
if isempty(xlow),   xlow=[1.11023d-16];  end;
if isempty(xlow1), xlow1=[1.490116d-8];  end;
if isempty(lnxmin), lnxmin=[-708.3964d0];  end;
x = xvalue;
if( x < zero ) ;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'ABRAM1 - Fatal error!');
writef(1,['%s','\n'], '  Argument X < 0.');
abram1Result = zero;
elseif( x == zero ) ;
abram1Result = half;
elseif( x < xlow ) ;
abram1Result = half;
elseif( x < xlow1 ) ;
abram1Result =( one - x ./ onerpi - x .* x .* log( x ) ) .* half;
elseif( x <= two ) ;
t =( x .* x ./ two - half ) - half;
[fval , ntermf, ab1f, t ]=cheval( ntermf, ab1f, t );
[gval , ntermg, ab1g, t ]=cheval( ntermg, ab1g, t );
[hval , ntermh, ab1h, t ]=cheval( ntermh, ab1h, t );
abram1Result = fval - x .*( gval ./ onerpi + x .* log( x ) .* hval );
else;
v = three .*(( x ./ two ) .^( two ./ three ) );
t =( six ./ v - half ) - half;
[asval , nterma, ab1as, t ]=cheval( nterma, ab1as, t );
asln = log( asval .* sqrt( v ./ three ) ./ rt3bpi ) - v;
if( asln < lnxmin ) ;
abram1Result = zero;
else;
abram1Result = exp( asln );
end;
end;
if ~isempty(inputname(1)), assignin('caller','FUntemp',xvalue); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end
function [n_data, x, fx]=abram1_values( n_data, x, fx ,varargin);
%*******************************************************************************
%
%! ABRAM1_VALUES returns some values of the Abramowitz1 function.
%
%  Discussion:
%
%    The function is defined by:
%
%      ABRAM1(x) = Integral ( 0 <= t < infinity ) t * exp ( -t^2 - x / t ) dt
%
%    The data was reported by McLeod.
%
%  Modified:
%
%    21 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real ( kind = 8 ) X, the argument of the function.
%
%    Output, real ( kind = 8 ) FX, the value of the function.
%
persistent fx_vec n_max x_vec 

if isempty(n_max), n_max = 20; end;
if isempty(fx_vec), fx_vec([1:n_max]) =[real( 0.49828219848799921792d+00),real( 0.49324391773047288556d+00),real( 0.47431612784691234649d+00),real( 0.41095983258760410149d+00),real( 0.25317617388227035867d+00),real( 0.14656338138597777543d+00),real( 0.11421547056018366587d+00),real( 0.90026307383483764795d-01),real( 0.64088214170742303375d-01),real( 0.57446614314166191085d-01),real( 0.51581624564800730959d-01),real( 0.25263719555776416016d-01),real( 0.11930803330196594536d-01),real( 0.59270542280915272465d-02),real( 0.30609215358017829567d-02),real( 0.16307382136979552833d-02),real( 0.28371851916959455295d-03),real( 0.21122150121323238154d-04),real( 0.20344578892601627337d-05),real( 0.71116517236209642290e-09) ]; end;
if isempty(x_vec), x_vec([1:n_max]) =[real( 0.0019531250d+00),real( 0.0078125000d+00),real( 0.0312500000d+00),real( 0.1250000000d+00),real( 0.5000000000d+00),real( 1.0000000000d+00),real( 1.2500000000d+00),real( 1.5000000000d+00),real( 1.8750000000d+00),real( 2.0000000000d+00),real( 2.1250000000d+00),real( 3.0000000000d+00),real( 4.0000000000d+00),real( 5.0000000000d+00),real( 6.0000000000d+00),real( 7.0000000000d+00),real( 10.0000000000d+00),real( 15.0000000000d+00),real( 20.0000000000d+00),real( 40.0000000000d+00) ]; end;
if( n_data < 0 ) ;
n_data = 0;
end;
n_data = n_data + 1;
if( n_max < n_data ) ;
n_data = 0;
x = 0.0d+00;
fx = 0.0d+00;
else;
x  = x_vec(n_data);
fx = fx_vec(n_data);
end;
return;
end
function [abram2Result, xvalue ]=abram2( xvalue ,varargin);
%*******************************************************************************
%
%! ABRAM2 evaluates the Abramowitz function of order 2.
%
%  Discussion:
%
%    The function is defined by:
%
%      ABRAM2(x) = Integral ( 0 <= t < infinity ) t^2 * exp ( -t^2 - x / t ) dt
%
%    The code uses Chebyshev expansions with the coefficients
%    given to an accuracy of 20 decimal places.
%
%    This subroutine is set up to work on IEEE machines.
%
%  Modified:
%
%    07 August 2004
%
%  Author:
%
%    Allan McLeod,
%    Department of Mathematics and Statistics,
%    Paisley University, High Street, Paisley, Scotland, PA12BE
%    macl_ms0@paisley.ac.uk
%
%  Reference:
%
%    Allan Mcleod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input, real ( kind = 8 ) XVALUE, the argument of the function.
%
%    Output, real ( kind = 8 ) ABRAM2, the value of the function.
%
global realmlv
abram2Result=[];
persistent ab2as ab2f ab2g ab2h asln asval fval gval half hval lnxmin nterma ntermf ntermg ntermh onerpi rt3bpi rtpib4 six t three two v x xlow xlow1 zero 

if isempty(half), half = 0.5d+00; end;
if isempty(nterma), nterma = 23; end;
if isempty(ntermf), ntermf = 9; end;
if isempty(ntermg), ntermg = 8; end;
if isempty(ntermh), ntermh = 7; end;
if isempty(six), six = real( 6.0d+00); end;
if isempty(three), three = 3.0d+00; end;
if isempty(two), two = 2.0d+00; end;
if isempty(x), x=0; end;
if isempty(zero), zero = 0.0d+00; end;
if isempty(asln), asln=0; end;
if isempty(asval), asval=0; end;
if isempty(fval), fval=0; end;
if isempty(gval), gval=0; end;
if isempty(hval), hval=0; end;
if isempty(t), t=0; end;
if isempty(v), v=0; end;
if isempty(ab2f),   ab2f=[1.03612162804243713846d0,0.19371246626794570012d0,-0.7258758839233007378d-1,0.174790590864327399d-2,-0.1281223233756549d-4,0.4115018153651d-7,-0.6971047256d-10,0.6990183d-13,-0.4492d-16,0.2d-19];  end;
if isempty(ab2g),   ab2g=[1.46290157198630741150d0,0.20189466883154014317d0,-0.2908292087997129022d-1,0.47061049035270050d-3,-0.257922080359333d-5,0.656133712946d-8,-0.914110203d-11,0.774276d-14,-0.429d-17];  end;
if isempty(ab2h),   ab2h=[0.30117225010910488881d0,-0.1588667818317623783d-1,0.19295936935584526d-3,-0.90199587849300d-6,0.206105041837d-8,-0.265111806d-11,0.210864d-14,-0.111d-17];  end;
try;ab2as(0+1);catch;   ab2as(0+1)=[2.46492325304334856893d0];  end;
try;ab2as(1+1);catch;   ab2as(1+1)=[0.23142797422248905432d0];  end;
try;ab2as(2+1);catch;   ab2as(2+1)=[-0.94068173010085773d-3];  end;
try;ab2as(3+1);catch;   ab2as(3+1)=[0.8290270038089733d-4];  end;
try;ab2as(4+1);catch;   ab2as(4+1)=[-0.883894704245866d-5];  end;
try;ab2as(5+1);catch;   ab2as(5+1)=[0.106638543567985d-5];  end;
try;ab2as(6+1);catch;   ab2as(6+1)=[-0.13991128538529d-6];  end;
try;ab2as(7+1);catch;   ab2as(7+1)=[0.1939793208445d-7];  end;
try;ab2as(8+1);catch;   ab2as(8+1)=[-0.277049938375d-8];  end;
try;ab2as(9+1);catch;   ab2as(9+1)=[0.39590687186d-9];  end;
try;ab2as(10+1);catch;   ab2as(10+1)=[-0.5408354342d-10];  end;
try;ab2as(11+1);catch;   ab2as(11+1)=[0.635546076d-11];  end;
try;ab2as(12+1);catch;   ab2as(12+1)=[-0.38461613d-12];  end;
try;ab2as(13+1);catch;   ab2as(13+1)=[-0.11696067d-12];  end;
try;ab2as(14+1);catch;   ab2as(14+1)=[0.6896671d-13];  end;
try;ab2as(15+1);catch;   ab2as(15+1)=[-0.2503113d-13];  end;
try;ab2as(16+1);catch;   ab2as(16+1)=[0.785586d-14];  end;
try;ab2as(17+1);catch;   ab2as(17+1)=[-0.230334d-14];  end;
try;ab2as(18+1);catch;   ab2as(18+1)=[0.64914d-15];  end;
try;ab2as(19+1);catch;   ab2as(19+1)=[-0.17797d-15];  end;
try;ab2as(20+1);catch;   ab2as(20+1)=[0.4766d-16];  end;
try;ab2as(21+1);catch;   ab2as(21+1)=[-0.1246d-16];  end;
try;ab2as(22+1);catch;   ab2as(22+1)=[0.316d-17];  end;
try;ab2as(23+1);catch;   ab2as(23+1)=[-0.77d-18];  end;
try;ab2as(24+1);catch;   ab2as(24+1)=[0.18d-18];  end;
try;ab2as(25+1);catch;   ab2as(25+1)=[-0.4d-19];  end;
try;ab2as(26+1);catch;   ab2as(26+1)=[0.1d-19];  end;
if isempty(rt3bpi),   rt3bpi=[0.97720502380583984317d0];  end;
if isempty(rtpib4),   rtpib4=[0.44311346272637900682d0];  end;
if isempty(onerpi),   onerpi=[0.56418958354775628695d0];  end;
%
%   Machine-dependent constants (suitable for IEEE machines)
%
if isempty(xlow),   xlow=[2.22045d-16];  end;
if isempty(xlow1), xlow1=[1.490116d-8];  end;
if isempty(lnxmin), lnxmin=[-708.3964d0];  end;
x = xvalue;
if( x < zero ) ;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'ABRAM2 - Fatal error!');
writef(1,['%s','\n'], '  Argument X < 0.');
abram2Result = zero;
elseif( x == zero ) ;
abram2Result = rtpib4;
elseif( x < xlow ) ;
abram2Result = rtpib4;
elseif( x < xlow1 ) ;
abram2Result = rtpib4 - half .* x + x .* x .* x .* log( x ) ./ six;
elseif( x <= 2.0d+00 ) ;
t =( x .* x ./ two - half ) - half;
[fval , ntermf, ab2f, t ]=cheval( ntermf, ab2f, t );
[gval , ntermg, ab2g, t ]=cheval( ntermg, ab2g, t );
[hval , ntermh, ab2h, t ]=cheval( ntermh, ab2h, t );
abram2Result = fval ./ onerpi + x .*( x .* x .* log( x ) .* hval - gval );
else;
v = three .*(( x ./ two ) .^( two ./ three ) );
t =( six ./ v - half ) - half;
[asval , nterma, ab2as, t ]=cheval( nterma, ab2as, t );
asln = log( asval ./ rt3bpi ) + log( v ./ three ) - v;
if( asln < lnxmin ) ;
abram2Result = zero;
else;
abram2Result = exp( asln );
end;
end;
if ~isempty(inputname(1)), assignin('caller','FUntemp',xvalue); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end
function [n_data, x, fx]=abram2_values( n_data, x, fx ,varargin);
%*******************************************************************************
%
%! ABRAM2_VALUES returns some values of the Abramowitz2 function.
%
%  Discussion:
%
%    The function is defined by:
%
%      ABRAM2(x) = Integral ( 0 <= t < infinity ) t^2 * exp ( -t^2 - x / t ) dt
%
%    The data was reported by McLeod.
%
%  Modified:
%
%    22 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real ( kind = 8 ) X, the argument of the function.
%
%    Output, real ( kind = 8 ) FX, the value of the function.
%
persistent fx_vec n_max x_vec 

if isempty(n_max), n_max = 20; end;
if isempty(fx_vec), fx_vec([1:n_max]) =[real( 0.44213858162107913430d+00),real( 0.43923379545684026308d+00),real( 0.42789857297092602234d+00),real( 0.38652825661854504406d+00),real( 0.26538204413231368110d+00),real( 0.16848734838334595000d+00),real( 0.13609200032513227112d+00),real( 0.11070330027727917352d+00),real( 0.82126019995530382267d-01),real( 0.74538781999594581763d-01),real( 0.67732034377612811390d-01),real( 0.35641808698811851022d-01),real( 0.17956589956618269083d-01),real( 0.94058737143575370625d-02),real( 0.50809356204299213556d-02),real( 0.28149565414209719359d-02),real( 0.53808696422559303431d-03),real( 0.44821756380146327259d-04),real( 0.46890678427324100410d-05),real( 0.20161544850996420504d-08) ]; end;
if isempty(x_vec), x_vec([1:n_max]) =[real( 0.0019531250d+00),real( 0.0078125000d+00),real( 0.0312500000d+00),real( 0.1250000000d+00),real( 0.5000000000d+00),real( 1.0000000000d+00),real( 1.2500000000d+00),real( 1.5000000000d+00),real( 1.8750000000d+00),real( 2.0000000000d+00),real( 2.1250000000d+00),real( 3.0000000000d+00),real( 4.0000000000d+00),real( 5.0000000000d+00),real( 6.0000000000d+00),real( 7.0000000000d+00),real( 10.0000000000d+00),real( 15.0000000000d+00),real( 20.0000000000d+00),real( 40.0000000000d+00) ]; end;
if( n_data < 0 ) ;
n_data = 0;
end;
n_data = n_data + 1;
if( n_max < n_data ) ;
n_data = 0;
x = 0.0d+00;
fx = 0.0d+00;
else;
x  = x_vec(n_data);
fx = fx_vec(n_data);
end;
return;
end
function [airy_ai_intResult, xvalue ]=airy_ai_int( xvalue ,varargin);
%*******************************************************************************
%
%! AIRY_AI_INT calculates the integral of the Airy function Ai.
%
%  Discussion:
%
%    The function is defined by:
%
%      AIRY_AI_INT(x) = Integral ( 0 <= t <= x ) Ai(t) dt
%
%    The program uses Chebyshev expansions, the coefficients of which
%    are given to 20 decimal places.
%
%    This subroutine is set up to work on IEEE machines.
%
%  Modified:
%
%    07 August 2004
%
%  Author:
%
%    Allan McLeod,
%    Department of Mathematics and Statistics,
%    Paisley University, High Street, Paisley, Scotland, PA12BE
%    macl_ms0@paisley.ac.uk
%
%  Reference:
%
%    Allan Mcleod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input, real ( kind = 8 ) XVALUE, the argument of the function.
%
%    Output, real ( kind = 8 ) AIRY_AI_INT, the value of the function.
%
airy_ai_intResult=[];
persistent aaint1 aaint2 aaint3 aaint4 aaint5 airzer arg eight forty1 four fr996 gval hval nine ninhun nterm1 nterm2 nterm3 nterm4 nterm5 one piby4 pitim6 rt2b3p t temp three two x xhigh1 xlow1 xneg1 z zero 

if isempty(arg), arg=0; end;
if isempty(eight), eight = 8.0d+00; end;
if isempty(four), four = 4.0d+00; end;
if isempty(gval), gval=0; end;
if isempty(hval), hval=0; end;
if isempty(nterm1), nterm1 = 22; end;
if isempty(nterm2), nterm2 = 17; end;
if isempty(nterm3), nterm3 = 37; end;
if isempty(one), one = 1.0d+00; end;
if isempty(t), t=0; end;
if isempty(temp), temp=0; end;
if isempty(three), three = 3.0d+00; end;
if isempty(two), two = 2.0d+00; end;
if isempty(x), x=0; end;
if isempty(zero), zero = 0.0d+00; end;
if isempty(z), z=0; end;
try;aaint1(0+1);catch;   aaint1(0+1)=[0.37713517694683695526d0];  end;
try;aaint1(1+1);catch;   aaint1(1+1)=[-0.13318868432407947431d0];  end;
try;aaint1(2+1);catch;   aaint1(2+1)=[0.3152497374782884809d-1];  end;
try;aaint1(3+1);catch;   aaint1(3+1)=[-0.318543076436574077d-2];  end;
try;aaint1(4+1);catch;   aaint1(4+1)=[-0.87398764698621915d-3];  end;
try;aaint1(5+1);catch;   aaint1(5+1)=[0.46699497655396971d-3];  end;
try;aaint1(6+1);catch;   aaint1(6+1)=[-0.9544936738983692d-4];  end;
try;aaint1(7+1);catch;   aaint1(7+1)=[0.542705687156716d-5];  end;
try;aaint1(8+1);catch;   aaint1(8+1)=[0.239496406252188d-5];  end;
try;aaint1(9+1);catch;   aaint1(9+1)=[-0.75690270205649d-6];  end;
try;aaint1(10+1);catch;   aaint1(10+1)=[0.9050138584518d-7];  end;
try;aaint1(11+1);catch;   aaint1(11+1)=[0.320529456043d-8];  end;
try;aaint1(12+1);catch;   aaint1(12+1)=[-0.303825536444d-8];  end;
try;aaint1(13+1);catch;   aaint1(13+1)=[0.48900118596d-9];  end;
try;aaint1(14+1);catch;   aaint1(14+1)=[-0.1839820572d-10];  end;
try;aaint1(15+1);catch;   aaint1(15+1)=[-0.711247519d-11];  end;
try;aaint1(16+1);catch;   aaint1(16+1)=[0.151774419d-11];  end;
try;aaint1(17+1);catch;   aaint1(17+1)=[-0.10801922d-12];  end;
try;aaint1(18+1);catch;   aaint1(18+1)=[-0.963542d-14];  end;
try;aaint1(19+1);catch;   aaint1(19+1)=[0.313425d-14];  end;
try;aaint1(20+1);catch;   aaint1(20+1)=[-0.29446d-15];  end;
try;aaint1(21+1);catch;   aaint1(21+1)=[-0.477d-17];  end;
try;aaint1(22+1);catch;   aaint1(22+1)=[0.461d-17];  end;
try;aaint1(23+1);catch;   aaint1(23+1)=[-0.53d-18];  end;
try;aaint1(24+1);catch;   aaint1(24+1)=[0.1d-19];  end;
try;aaint1(25+1);catch;   aaint1(25+1)=[0.1d-19];  end;
try;aaint2(0+1);catch;   aaint2(0+1)=[1.92002524081984009769d0];  end;
try;aaint2(1+1);catch;   aaint2(1+1)=[-0.4220049417256287021d-1];  end;
try;aaint2(2+1);catch;   aaint2(2+1)=[-0.239457722965939223d-2];  end;
try;aaint2(3+1);catch;   aaint2(3+1)=[-0.19564070483352971d-3];  end;
try;aaint2(4+1);catch;   aaint2(4+1)=[-0.1547252891056112d-4];  end;
try;aaint2(5+1);catch;   aaint2(5+1)=[-0.140490186137889d-5];  end;
try;aaint2(6+1);catch;   aaint2(6+1)=[-0.12128014271367d-6];  end;
try;aaint2(7+1);catch;   aaint2(7+1)=[-0.1179186050192d-7];  end;
try;aaint2(8+1);catch;   aaint2(8+1)=[-0.104315578788d-8];  end;
try;aaint2(9+1);catch;   aaint2(9+1)=[-0.10908209293d-9];  end;
try;aaint2(10+1);catch;   aaint2(10+1)=[-0.929633045d-11];  end;
try;aaint2(11+1);catch;   aaint2(11+1)=[-0.110946520d-11];  end;
try;aaint2(12+1);catch;   aaint2(12+1)=[-0.7816483d-13];  end;
try;aaint2(13+1);catch;   aaint2(13+1)=[-0.1319661d-13];  end;
try;aaint2(14+1);catch;   aaint2(14+1)=[-0.36823d-15];  end;
try;aaint2(15+1);catch;   aaint2(15+1)=[-0.21505d-15];  end;
try;aaint2(16+1);catch;   aaint2(16+1)=[0.1238d-16];  end;
try;aaint2(17+1);catch;   aaint2(17+1)=[-0.557d-17];  end;
try;aaint2(18+1);catch;   aaint2(18+1)=[0.84d-18];  end;
try;aaint2(19+1);catch;   aaint2(19+1)=[-0.21d-18];  end;
try;aaint2(20+1);catch;   aaint2(20+1)=[0.4d-19];  end;
try;aaint2(21+1);catch;   aaint2(21+1)=[-0.1d-19];  end;
try;aaint3(0+1);catch;   aaint3(0+1)=[0.47985893264791052053d0];  end;
try;aaint3(1+1);catch;   aaint3(1+1)=[-0.19272375126169608863d0];  end;
try;aaint3(2+1);catch;   aaint3(2+1)=[0.2051154129525428189d-1];  end;
try;aaint3(3+1);catch;   aaint3(3+1)=[0.6332000070732488786d-1];  end;
try;aaint3(4+1);catch;   aaint3(4+1)=[-0.5093322261845754082d-1];  end;
try;aaint3(5+1);catch;   aaint3(5+1)=[0.1284424078661663016d-1];  end;
try;aaint3(6+1);catch;   aaint3(6+1)=[0.2760137088989479413d-1];  end;
try;aaint3(7+1);catch;   aaint3(7+1)=[-0.1547066673866649507d-1];  end;
try;aaint3(8+1);catch;   aaint3(8+1)=[-0.1496864655389316026d-1];  end;
try;aaint3(9+1);catch;   aaint3(9+1)=[0.336617614173574541d-2];  end;
try;aaint3(10+1);catch;   aaint3(10+1)=[0.530851163518892985d-2];  end;
try;aaint3(11+1);catch;   aaint3(11+1)=[0.41371226458555081d-3];  end;
try;aaint3(12+1);catch;   aaint3(12+1)=[-0.102490579926726266d-2];  end;
try;aaint3(13+1);catch;   aaint3(13+1)=[-0.32508221672025853d-3];  end;
try;aaint3(14+1);catch;   aaint3(14+1)=[0.8608660957169213d-4];  end;
try;aaint3(15+1);catch;   aaint3(15+1)=[0.6671367298120775d-4];  end;
try;aaint3(16+1);catch;   aaint3(16+1)=[0.449205999318095d-5];  end;
try;aaint3(17+1);catch;   aaint3(17+1)=[-0.670427230958249d-5];  end;
try;aaint3(18+1);catch;   aaint3(18+1)=[-0.196636570085009d-5];  end;
try;aaint3(19+1);catch;   aaint3(19+1)=[0.22229677407226d-6];  end;
try;aaint3(20+1);catch;   aaint3(20+1)=[0.22332222949137d-6];  end;
try;aaint3(21+1);catch;   aaint3(21+1)=[0.2803313766457d-7];  end;
try;aaint3(22+1);catch;   aaint3(22+1)=[-0.1155651663619d-7];  end;
try;aaint3(23+1);catch;   aaint3(23+1)=[-0.433069821736d-8];  end;
try;aaint3(24+1);catch;   aaint3(24+1)=[-0.6227777938d-10];  end;
try;aaint3(25+1);catch;   aaint3(25+1)=[0.26432664903d-9];  end;
try;aaint3(26+1);catch;   aaint3(26+1)=[0.5333881114d-10];  end;
try;aaint3(27+1);catch;   aaint3(27+1)=[-0.522957269d-11];  end;
try;aaint3(28+1);catch;   aaint3(28+1)=[-0.382229283d-11];  end;
try;aaint3(29+1);catch;   aaint3(29+1)=[-0.40958233d-12];  end;
try;aaint3(30+1);catch;   aaint3(30+1)=[0.11515622d-12];  end;
try;aaint3(31+1);catch;   aaint3(31+1)=[0.3875766d-13];  end;
try;aaint3(32+1);catch;   aaint3(32+1)=[0.140283d-14];  end;
try;aaint3(33+1);catch;   aaint3(33+1)=[-0.141526d-14];  end;
try;aaint3(34+1);catch;   aaint3(34+1)=[-0.28746d-15];  end;
try;aaint3(35+1);catch;   aaint3(35+1)=[0.923d-17];  end;
try;aaint3(36+1);catch;   aaint3(36+1)=[0.1224d-16];  end;
try;aaint3(37+1);catch;   aaint3(37+1)=[0.157d-17];  end;
try;aaint3(38+1);catch;   aaint3(38+1)=[-0.19d-18];  end;
try;aaint3(39+1);catch;   aaint3(39+1)=[-0.8d-19];  end;
try;aaint3(40+1);catch;   aaint3(40+1)=[-0.1d-19];  end;
if isempty(aaint4),   aaint4=[1.99653305828522730048d0,-0.187541177605417759d-2,-0.15377536280305750d-3,-0.1283112967682349d-4,-0.108128481964162d-5,-0.9182131174057d-7,-0.784160590960d-8,-0.67292453878d-9,-0.5796325198d-10,-0.501040991d-11,-0.43420222d-12,-0.3774305d-13,-0.328473d-14,-0.28700d-15,-0.2502d-16,-0.220d-17,-0.19d-18,-0.2d-19];  end;
if isempty(aaint5),   aaint5=[1.13024602034465716133d0,-0.464718064639872334d-2,-0.35137413382693203d-3,-0.2768117872545185d-4,-0.222057452558107d-5,-0.18089142365974d-6,-0.1487613383373d-7,-0.123515388168d-8,-0.10310104257d-9,-0.867493013d-11,-0.73080054d-12,-0.6223561d-13,-0.525128d-14,-0.45677d-15,-0.3748d-16,-0.356d-17,-0.23d-18,-0.4d-19];  end;
if isempty(nine),   nine=[9.0d0];  end;
if isempty(forty1), forty1=[41.0d0];  end;
if isempty(ninhun),   ninhun=[900.0d0];  end;
if isempty(fr996), fr996=[4996.0d0];  end;
if isempty(piby4),   piby4=[0.78539816339744830962d0];  end;
if isempty(pitim6),   pitim6=[18.84955592153875943078d0];  end;
if isempty(rt2b3p),   rt2b3p=[0.46065886596178063902d0];  end;
if isempty(airzer),   airzer=[0.35502805388781723926d0];  end;
%
%   Machine-dependant constants (suitable for IEEE machines)
%
if isempty(nterm4),   nterm4=[15];  end;
if isempty(nterm5), nterm5=[15];  end;
if isempty(xlow1),   xlow1=[2.22045d-16];  end;
if isempty(xhigh1), xhigh1=[14.480884d0];  end;
if isempty(xneg1), xneg1=[-2.727134d10];  end;
x = xvalue;
if( x < xneg1 ) ;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'AIRY_AI_INT - Fatal error!');
writef(1,['%s','\n'], '  X too negative for accurate computation.');
airy_ai_intResult = -two ./ three;
if ~isempty(inputname(1)), assignin('caller','FUntemp',xvalue); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
elseif( x < -eight ) ;
z = -( x + x ) .* sqrt( -x ) ./ three;
arg = z + piby4;
temp = nine .* z .* z;
t =( fr996 - temp ) ./( ninhun + temp );
[gval , nterm4, aaint4, t ]=cheval( nterm4, aaint4, t );
[hval , nterm5, aaint5, t ]=cheval( nterm5, aaint5, t );
temp = gval .* cos( arg ) + hval .* sin( arg ) ./ z;
airy_ai_intResult = rt2b3p .* temp ./ sqrt( z ) - two ./ three;
elseif( x <= -xlow1 );
t = -x ./ four - one;
airy_ai_intResult = x .* cheval( nterm3, aaint3, t );
elseif( x < xlow1 ) ;
airy_ai_intResult = airzer .* x;
elseif( x <= four ) ;
t = x ./ two - one;
airy_ai_intResult = cheval( nterm1, aaint1, t ) .* x;
elseif( x <= xhigh1 ) ;
z =( x + x ) .* sqrt( x ) ./ three;
temp = three .* z;
t =( forty1 - temp ) ./( nine + temp );
temp = exp( -z ) .* cheval( nterm2, aaint2, t ) ./ sqrt( pitim6 .* z );
airy_ai_intResult = one ./ three - temp;
else;
airy_ai_intResult = one ./ three;
end;
if ~isempty(inputname(1)), assignin('caller','FUntemp',xvalue); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end
function [n_data, x, fx]=airy_ai_int_values( n_data, x, fx ,varargin);
%*******************************************************************************
%
%! AIRY_AI_INT_VALUES returns some values of the integral of the Airy function.
%
%  Discussion:
%
%    The function is defined by:
%
%      AIRY_AI_INT(x) = Integral ( 0 <= t <= x ) Ai(t) dt
%
%    The data was reported by McLeod.
%
%  Modified:
%
%    22 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real ( kind = 8 ) X, the argument of the function.
%
%    Output, real ( kind = 8 ) FX, the value of the function.
%
persistent fx_vec n_max x_vec 

if isempty(n_max), n_max = 20; end;
if isempty(fx_vec), fx_vec([1:n_max]) =[real( -0.75228838916610124300d+00),real( -0.57348350185854889466d+00),real( -0.76569840313421291743d+00),real( -0.65181015505382467421d+00),real( -0.55881974894471876922d+00),real( -0.56902352870716815309d+00),real( -0.47800749642926168100d+00),real( -0.46567398346706861416d+00),real( -0.96783140945618013679d-01),real( -0.34683049857035607494d-03),real(  0.34658366917927930790d-03),real(  0.27657581846051227124d-02),real(  0.14595330491185717833d+00),real(  0.23631734191710977960d+00),real(  0.33289264538612212697d+00),real(  0.33318759129779422976d+00),real(  0.33332945170523851439d+00),real(  0.33333331724248357420d+00),real(  0.33333333329916901594d+00),real(  0.33333333333329380187d+00) ]; end;
if isempty(x_vec), x_vec([1:n_max]) =[real( -12.0000000000d+00),real( -11.0000000000d+00),real( -10.0000000000d+00),real(  -9.5000000000d+00),real(  -9.0000000000d+00),real(  -6.5000000000d+00),real(  -4.0000000000d+00),real(  -1.0000000000d+00),real(  -0.2500000000d+00),real(  -0.0009765625d+00),real(   0.0009765625d+00),real(   0.0078125000d+00),real(   0.5000000000d+00),real(   1.0000000000d+00),real(   4.0000000000d+00),real(   4.5000000000d+00),real(   6.0000000000d+00),real(   8.0000000000d+00),real(  10.0000000000d+00),real(  12.0000000000d+00) ]; end;
if( n_data < 0 ) ;
n_data = 0;
end;
n_data = n_data + 1;
if( n_max < n_data ) ;
n_data = 0;
x = 0.0d+00;
fx = 0.0d+00;
else;
x  = x_vec(n_data);
fx = fx_vec(n_data);
end;
return;
end
function [airy_bi_intResult, xvalue ]=airy_bi_int( xvalue ,varargin);
%*******************************************************************************
%
%! AIRY_BI_INT calculates the integral of the Airy function Bi.
%
%  Discussion:
%
%    The function is defined by:
%
%      AIRY_BI_INT(x) = Integral ( 0 <= t <= x ) Bi(t) dt
%
%    The program uses Chebyshev expansions, the coefficients of which
%    are given to 20 decimal places.
%
%    This subroutine is set up to work on IEEE machines.
%
%  Modified:
%
%    07 August 2004
%
%  Author:
%
%    Allan McLeod,
%    Department of Mathematics and Statistics,
%    Paisley University, High Street, Paisley, Scotland, PA12BE
%    macl_ms0@paisley.ac.uk
%
%  Reference:
%
%    Allan Mcleod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input, real ( kind = 8 ) XVALUE, the argument of the function.
%
%    Output, real ( kind = 8 ) AIRY_BI_INT, the value of the function.
%
global realmlv
airy_bi_intResult=[];
persistent abint1 abint2 abint3 abint4 abint5 arg birzer eight f1 f2 four nine ninhun nterm1 nterm2 nterm3 nterm4 nterm5 one onept5 piby4 rt2b3p seven sixten t temp thr644 three x xhigh1 xlow1 xmax xneg1 z zero 

if isempty(eight), eight = real( 8.0); end;
if isempty(four), four = real( 4.0d+00); end;
if isempty(nterm1), nterm1 = 33; end;
if isempty(nterm2), nterm2 = 30; end;
if isempty(nterm3), nterm3 = 34; end;
if isempty(one), one = 1.0d+00; end;
if isempty(three), three = 3.0d+00; end;
if isempty(x), x=0; end;
if isempty(zero), zero = 0.0d+00; end;
if isempty(arg), arg=0; end;
if isempty(f1), f1=0; end;
if isempty(f2), f2=0; end;
if isempty(t), t=0; end;
if isempty(temp), temp=0; end;
if isempty(z), z=0; end;
try;abint1(0+1);catch;   abint1(0+1)=[0.38683352445038543350d0];  end;
try;abint1(1+1);catch;   abint1(1+1)=[-0.8823213550888908821d-1];  end;
try;abint1(2+1);catch;   abint1(2+1)=[0.21463937440355429239d0];  end;
try;abint1(3+1);catch;   abint1(3+1)=[-0.4205347375891315126d-1];  end;
try;abint1(4+1);catch;   abint1(4+1)=[0.5932422547496086771d-1];  end;
try;abint1(5+1);catch;   abint1(5+1)=[-0.840787081124270210d-2];  end;
try;abint1(6+1);catch;   abint1(6+1)=[0.871824772778487955d-2];  end;
try;abint1(7+1);catch;   abint1(7+1)=[-0.12191600199613455d-3];  end;
try;abint1(8+1);catch;   abint1(8+1)=[0.44024821786023234d-3];  end;
try;abint1(9+1);catch;   abint1(9+1)=[0.27894686666386678d-3];  end;
try;abint1(10+1);catch;   abint1(10+1)=[-0.7052804689785537d-4];  end;
try;abint1(11+1);catch;   abint1(11+1)=[0.5901080066770100d-4];  end;
try;abint1(12+1);catch;   abint1(12+1)=[-0.1370862587982142d-4];  end;
try;abint1(13+1);catch;   abint1(13+1)=[0.505962573749073d-5];  end;
try;abint1(14+1);catch;   abint1(14+1)=[-0.51598837766735d-6];  end;
try;abint1(15+1);catch;   abint1(15+1)=[0.397511312349d-8];  end;
try;abint1(16+1);catch;   abint1(16+1)=[0.9524985978055d-7];  end;
try;abint1(17+1);catch;   abint1(17+1)=[-0.3681435887321d-7];  end;
try;abint1(18+1);catch;   abint1(18+1)=[0.1248391688136d-7];  end;
try;abint1(19+1);catch;   abint1(19+1)=[-0.249097619137d-8];  end;
try;abint1(20+1);catch;   abint1(20+1)=[0.31775245551d-9];  end;
try;abint1(21+1);catch;   abint1(21+1)=[0.5434365270d-10];  end;
try;abint1(22+1);catch;   abint1(22+1)=[-0.4024566915d-10];  end;
try;abint1(23+1);catch;   abint1(23+1)=[0.1393855527d-10];  end;
try;abint1(24+1);catch;   abint1(24+1)=[-0.303817509d-11];  end;
try;abint1(25+1);catch;   abint1(25+1)=[0.40809511d-12];  end;
try;abint1(26+1);catch;   abint1(26+1)=[0.1634116d-13];  end;
try;abint1(27+1);catch;   abint1(27+1)=[-0.2683809d-13];  end;
try;abint1(28+1);catch;   abint1(28+1)=[0.896641d-14];  end;
try;abint1(29+1);catch;   abint1(29+1)=[-0.183089d-14];  end;
try;abint1(30+1);catch;   abint1(30+1)=[0.21333d-15];  end;
try;abint1(31+1);catch;   abint1(31+1)=[0.1108d-16];  end;
try;abint1(32+1);catch;   abint1(32+1)=[-0.1276d-16];  end;
try;abint1(33+1);catch;   abint1(33+1)=[0.363d-17];  end;
try;abint1(34+1);catch;   abint1(34+1)=[-0.62d-18];  end;
try;abint1(35+1);catch;   abint1(35+1)=[0.5d-19];  end;
try;abint1(36+1);catch;   abint1(36+1)=[0.1d-19];  end;
try;abint2(0+1);catch;   abint2(0+1)=[2.04122078602516135181d0];  end;
try;abint2(1+1);catch;   abint2(1+1)=[0.2124133918621221230d-1];  end;
try;abint2(2+1);catch;   abint2(2+1)=[0.66617599766706276d-3];  end;
try;abint2(3+1);catch;   abint2(3+1)=[0.3842047982808254d-4];  end;
try;abint2(4+1);catch;   abint2(4+1)=[0.362310366020439d-5];  end;
try;abint2(5+1);catch;   abint2(5+1)=[0.50351990115074d-6];  end;
try;abint2(6+1);catch;   abint2(6+1)=[0.7961648702253d-7];  end;
try;abint2(7+1);catch;   abint2(7+1)=[0.717808442336d-8];  end;
try;abint2(8+1);catch;   abint2(8+1)=[-0.267770159104d-8];  end;
try;abint2(9+1);catch;   abint2(9+1)=[-0.168489514699d-8];  end;
try;abint2(10+1);catch;   abint2(10+1)=[-0.36811757255d-9];  end;
try;abint2(11+1);catch;   abint2(11+1)=[0.4757128727d-10];  end;
try;abint2(12+1);catch;   abint2(12+1)=[0.5263621945d-10];  end;
try;abint2(13+1);catch;   abint2(13+1)=[0.778973500d-11];  end;
try;abint2(14+1);catch;   abint2(14+1)=[-0.460546143d-11];  end;
try;abint2(15+1);catch;   abint2(15+1)=[-0.183433736d-11];  end;
try;abint2(16+1);catch;   abint2(16+1)=[0.32191249d-12];  end;
try;abint2(17+1);catch;   abint2(17+1)=[0.29352060d-12];  end;
try;abint2(18+1);catch;   abint2(18+1)=[-0.1657935d-13];  end;
try;abint2(19+1);catch;   abint2(19+1)=[-0.4483808d-13];  end;
try;abint2(20+1);catch;   abint2(20+1)=[0.27907d-15];  end;
try;abint2(21+1);catch;   abint2(21+1)=[0.711921d-14];  end;
try;abint2(22+1);catch;   abint2(22+1)=[-0.1042d-16];  end;
try;abint2(23+1);catch;   abint2(23+1)=[-0.119591d-14];  end;
try;abint2(24+1);catch;   abint2(24+1)=[0.4606d-16];  end;
try;abint2(25+1);catch;   abint2(25+1)=[0.20884d-15];  end;
try;abint2(26+1);catch;   abint2(26+1)=[-0.2416d-16];  end;
try;abint2(27+1);catch;   abint2(27+1)=[-0.3638d-16];  end;
try;abint2(28+1);catch;   abint2(28+1)=[0.863d-17];  end;
try;abint2(29+1);catch;   abint2(29+1)=[0.591d-17];  end;
try;abint2(30+1);catch;   abint2(30+1)=[-0.256d-17];  end;
try;abint2(31+1);catch;   abint2(31+1)=[-0.77d-18];  end;
try;abint2(32+1);catch;   abint2(32+1)=[0.66d-18];  end;
try;abint2(33+1);catch;   abint2(33+1)=[0.3d-19];  end;
try;abint2(34+1);catch;   abint2(34+1)=[-0.15d-18];  end;
try;abint2(35+1);catch;   abint2(35+1)=[0.2d-19];  end;
try;abint2(36+1);catch;   abint2(36+1)=[0.3d-19];  end;
try;abint2(37+1);catch;   abint2(37+1)=[-0.1d-19];  end;
try;abint3(0+1);catch;   abint3(0+1)=[0.31076961598640349251d0];  end;
try;abint3(1+1);catch;   abint3(1+1)=[-0.27528845887452542718d0];  end;
try;abint3(2+1);catch;   abint3(2+1)=[0.17355965706136543928d0];  end;
try;abint3(3+1);catch;   abint3(3+1)=[-0.5544017909492843130d-1];  end;
try;abint3(4+1);catch;   abint3(4+1)=[-0.2251265478295950941d-1];  end;
try;abint3(5+1);catch;   abint3(5+1)=[0.4107347447812521894d-1];  end;
try;abint3(6+1);catch;   abint3(6+1)=[0.984761275464262480d-2];  end;
try;abint3(7+1);catch;   abint3(7+1)=[-0.1555618141666041932d-1];  end;
try;abint3(8+1);catch;   abint3(8+1)=[-0.560871870730279234d-2];  end;
try;abint3(9+1);catch;   abint3(9+1)=[0.246017783322230475d-2];  end;
try;abint3(10+1);catch;   abint3(10+1)=[0.165740392292336978d-2];  end;
try;abint3(11+1);catch;   abint3(11+1)=[-0.3277587501435402d-4];  end;
try;abint3(12+1);catch;   abint3(12+1)=[-0.24434680860514925d-3];  end;
try;abint3(13+1);catch;   abint3(13+1)=[-0.5035305196152321d-4];  end;
try;abint3(14+1);catch;   abint3(14+1)=[0.1630264722247854d-4];  end;
try;abint3(15+1);catch;   abint3(15+1)=[0.851914057780934d-5];  end;
try;abint3(16+1);catch;   abint3(16+1)=[0.29790363004664d-6];  end;
try;abint3(17+1);catch;   abint3(17+1)=[-0.64389707896401d-6];  end;
try;abint3(18+1);catch;   abint3(18+1)=[-0.15046988145803d-6];  end;
try;abint3(19+1);catch;   abint3(19+1)=[0.1587013535823d-7];  end;
try;abint3(20+1);catch;   abint3(20+1)=[0.1276766299622d-7];  end;
try;abint3(21+1);catch;   abint3(21+1)=[0.140578534199d-8];  end;
try;abint3(22+1);catch;   abint3(22+1)=[-0.46564739741d-9];  end;
try;abint3(23+1);catch;   abint3(23+1)=[-0.15682748791d-9];  end;
try;abint3(24+1);catch;   abint3(24+1)=[-0.403893560d-11];  end;
try;abint3(25+1);catch;   abint3(25+1)=[0.666708192d-11];  end;
try;abint3(26+1);catch;   abint3(26+1)=[0.128869380d-11];  end;
try;abint3(27+1);catch;   abint3(27+1)=[-0.6968663d-13];  end;
try;abint3(28+1);catch;   abint3(28+1)=[-0.6254319d-13];  end;
try;abint3(29+1);catch;   abint3(29+1)=[-0.718392d-14];  end;
try;abint3(30+1);catch;   abint3(30+1)=[0.115296d-14];  end;
try;abint3(31+1);catch;   abint3(31+1)=[0.42276d-15];  end;
try;abint3(32+1);catch;   abint3(32+1)=[0.2493d-16];  end;
try;abint3(33+1);catch;   abint3(33+1)=[-0.971d-17];  end;
try;abint3(34+1);catch;   abint3(34+1)=[-0.216d-17];  end;
try;abint3(35+1);catch;   abint3(35+1)=[-0.2d-19];  end;
try;abint3(36+1);catch;   abint3(36+1)=[0.6d-19];  end;
try;abint3(37+1);catch;   abint3(37+1)=[0.1d-19];  end;
try;abint4(0+1);catch;   abint4(0+1)=[1.99507959313352047614d0];  end;
try;abint4(1+1);catch;   abint4(1+1)=[-0.273736375970692738d-2];  end;
try;abint4(2+1);catch;   abint4(2+1)=[-0.30897113081285850d-3];  end;
try;abint4(3+1);catch;   abint4(3+1)=[-0.3550101982798577d-4];  end;
try;abint4(4+1);catch;   abint4(4+1)=[-0.412179271520133d-5];  end;
try;abint4(5+1);catch;   abint4(5+1)=[-0.48235892316833d-6];  end;
try;abint4(6+1);catch;   abint4(6+1)=[-0.5678730727927d-7];  end;
try;abint4(7+1);catch;   abint4(7+1)=[-0.671874810365d-8];  end;
try;abint4(8+1);catch;   abint4(8+1)=[-0.79811649857d-9];  end;
try;abint4(9+1);catch;   abint4(9+1)=[-0.9514271478d-10];  end;
try;abint4(10+1);catch;   abint4(10+1)=[-0.1137468966d-10];  end;
try;abint4(11+1);catch;   abint4(11+1)=[-0.136359969d-11];  end;
try;abint4(12+1);catch;   abint4(12+1)=[-0.16381418d-12];  end;
try;abint4(13+1);catch;   abint4(13+1)=[-0.1972575d-13];  end;
try;abint4(14+1);catch;   abint4(14+1)=[-0.237844d-14];  end;
try;abint4(15+1);catch;   abint4(15+1)=[-0.28752d-15];  end;
try;abint4(16+1);catch;   abint4(16+1)=[-0.3475d-16];  end;
try;abint4(17+1);catch;   abint4(17+1)=[-0.422d-17];  end;
try;abint4(18+1);catch;   abint4(18+1)=[-0.51d-18];  end;
try;abint4(19+1);catch;   abint4(19+1)=[-0.6d-19];  end;
try;abint4(20+1);catch;   abint4(20+1)=[-0.1d-19];  end;
try;abint5(0+1);catch;   abint5(0+1)=[1.12672081961782566017d0];  end;
try;abint5(1+1);catch;   abint5(1+1)=[-0.671405567525561198d-2];  end;
try;abint5(2+1);catch;   abint5(2+1)=[-0.69812918017832969d-3];  end;
try;abint5(3+1);catch;   abint5(3+1)=[-0.7561689886425276d-4];  end;
try;abint5(4+1);catch;   abint5(4+1)=[-0.834985574510207d-5];  end;
try;abint5(5+1);catch;   abint5(5+1)=[-0.93630298232480d-6];  end;
try;abint5(6+1);catch;   abint5(6+1)=[-0.10608556296250d-6];  end;
try;abint5(7+1);catch;   abint5(7+1)=[-0.1213128916741d-7];  end;
try;abint5(8+1);catch;   abint5(8+1)=[-0.139631129765d-8];  end;
try;abint5(9+1);catch;   abint5(9+1)=[-0.16178918054d-9];  end;
try;abint5(10+1);catch;   abint5(10+1)=[-0.1882307907d-10];  end;
try;abint5(11+1);catch;   abint5(11+1)=[-0.220272985d-11];  end;
try;abint5(12+1);catch;   abint5(12+1)=[-0.25816189d-12];  end;
try;abint5(13+1);catch;   abint5(13+1)=[-0.3047964d-13];  end;
try;abint5(14+1);catch;   abint5(14+1)=[-0.358370d-14];  end;
try;abint5(15+1);catch;   abint5(15+1)=[-0.42831d-15];  end;
try;abint5(16+1);catch;   abint5(16+1)=[-0.4993d-16];  end;
try;abint5(17+1);catch;   abint5(17+1)=[-0.617d-17];  end;
try;abint5(18+1);catch;   abint5(18+1)=[-0.68d-18];  end;
try;abint5(19+1);catch;   abint5(19+1)=[-0.10d-18];  end;
try;abint5(20+1);catch;   abint5(20+1)=[-0.1d-19];  end;
if isempty(onept5),   onept5=[1.5d0];  end;
if isempty(seven),   seven=[7.0d0];  end;
if isempty(nine),   nine=[9.0d0];  end;
if isempty(sixten), sixten=[16.0d0];  end;
if isempty(ninhun),   ninhun=[900.0d0];  end;
if isempty(thr644), thr644=[3644.0d0];  end;
if isempty(piby4),   piby4=[0.78539816339744830962d0];  end;
if isempty(rt2b3p),   rt2b3p=[0.46065886596178063902d0];  end;
if isempty(birzer),   birzer=[0.61492662744600073515d0];  end;
%
%   Machine-dependent parameters (suitable for IEEE machines)
%
if isempty(nterm4),   nterm4=[17];  end;
if isempty(nterm5), nterm5=[17];  end;
if isempty(xlow1),   xlow1=[2.22044604d-16];  end;
if isempty(xhigh1), xhigh1=[104.587632d0];  end;
if isempty(xneg1),   xneg1=[-2.727134d10];  end;
if isempty(xmax), xmax=[1.79d308];  end;
x = xvalue;
if( x < xneg1 ) ;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'AIRY_BI_INT - Warning!');
writef(1,['%s','\n'], '  Argument is too negative for accurate computation.');
airy_bi_intResult = zero;
elseif( x < -seven ) ;
z = -( x + x ) .* sqrt( -x ) ./ three;
arg = z + piby4;
temp = nine .* z .* z;
t =( thr644 - temp ) ./( ninhun + temp );
f1 = cheval( nterm4, abint4, t ) .* sin( arg );
f2 = cheval( nterm5, abint5, t ) .* cos( arg ) ./ z;
airy_bi_intResult =( f2 - f1 ) .* rt2b3p ./ sqrt( z );
elseif( x <= -xlow1 ) ;
t = -( x + x ) ./ seven - one;
airy_bi_intResult = x .* cheval( nterm3, abint3, t );
elseif( x < xlow1 ) ;
airy_bi_intResult = birzer .* x;
elseif( x <= eight ) ;
t = x ./ four - one;
airy_bi_intResult = x .* exp( onept5 .* x ) .* cheval( nterm1, abint1, t );
elseif( x <= xhigh1 ) ;
t = sixten .* sqrt( eight ./ x ) ./ x - one;
z =( x + x ) .* sqrt( x ) ./ three;
temp = rt2b3p .* cheval( nterm2, abint2, t ) ./ sqrt( z );
temp = z + log( temp );
airy_bi_intResult = exp( temp );
else;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'AIRY_BI_INT - Warning!');
writef(1,['%s','\n'], '  Argument is too large for accurate computation.');
airy_bi_intResult = xmax;
end;
if ~isempty(inputname(1)), assignin('caller','FUntemp',xvalue); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end
function [n_data, x, fx]=airy_bi_int_values( n_data, x, fx ,varargin);
%*******************************************************************************
%
%! AIRY_BI_INT_VALUES returns some values of the integral of the Airy function.
%
%  Discussion:
%
%    The function is defined by:
%
%      AIRY_BI_INT(x) = Integral ( 0 <= t <= x ) Bi(t) dt
%
%    The data was reported by McLeod.
%
%  Modified:
%
%    23 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real ( kind = 8 ) X, the argument of the function.
%
%    Output, real ( kind = 8 ) FX, the value of the function.
%
persistent fx_vec n_max x_vec 

if isempty(n_max), n_max = 20; end;
if isempty(fx_vec), fx_vec([1:n_max]) =[real(  0.17660819031554631869d-01),real( -0.15040424806140020451d-01),real(  0.14756446293227661920d-01),real( -0.11847304264848446271d+00),real( -0.64916741266165856037d-01),real(  0.97260832464381044540d-01),real(  0.50760058495287539119d-01),real( -0.37300500963429492179d+00),real( -0.13962988442666578531d+00),real( -0.12001735266723296160d-02),real(  0.12018836117890354598d-02),real(  0.36533846550952011043d+00),real(  0.87276911673800812196d+00),real(  0.48219475263803429675d+02),real(  0.44006525804904178439d+06),real(  0.17608153976228301458d+07),real(  0.73779211705220007228d+07),real(  0.14780980310740671617d+09),real(  0.97037614223613433849d+11),real(  0.11632737638809878460d+15) ]; end;
if isempty(x_vec), x_vec([1:n_max]) =[real( -12.0000000000d+00),real( -10.0000000000d+00),real(  -8.0000000000d+00),real(  -7.5000000000d+00),real(  -7.0000000000d+00),real(  -6.5000000000d+00),real(  -4.0000000000d+00),real(  -1.0000000000d+00),real(  -0.2500000000d+00),real(  -0.0019531250d+00),real(   0.0019531250d+00),real(   0.5000000000d+00),real(   1.0000000000d+00),real(   4.0000000000d+00),real(   8.0000000000d+00),real(   8.5000000000d+00),real(   9.0000000000d+00),real(  10.0000000000d+00),real(  12.0000000000d+00),real(  14.0000000000d+00) ]; end;
if( n_data < 0 ) ;
n_data = 0;
end;
n_data = n_data + 1;
if( n_max < n_data ) ;
n_data = 0;
x = 0.0d+00;
fx = 0.0d+00;
else;
x  = x_vec(n_data);
fx = fx_vec(n_data);
end;
return;
end
function [airy_giResult, xvalue ]=airy_gi( xvalue ,varargin);
%*******************************************************************************
%
%! AIRY_GI computes the modified Airy function Gi(x).
%
%  Discussion:
%
%    The function is defined by:
%
%      AIRY_GI(x) = Integral ( 0 <= t < infinity ) sin ( x*t+t^3/3) dt / pi
%
%    The approximation uses Chebyshev expansions with the coefficients
%    given to 20 decimal places.
%
%    This subroutine is set up to work on IEEE machines.
%
%  Modified:
%
%    07 August 2004
%
%  Author:
%
%    Allan McLeod,
%    Department of Mathematics and Statistics,
%    Paisley University, High Street, Paisley, Scotland, PA12BE
%    macl_ms0@paisley.ac.uk
%
%  Reference:
%
%    Allan Mcleod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input, real ( kind = 8 ) XVALUE, the argument of the function.
%
%    Output, real ( kind = 8 ) AIRY_GI, the value of the function.
%
global realmlv
airy_giResult=[];
persistent arbin1 arbin2 argin1 argip1 argip2 arhin1 bi cheb1 cheb2 cosz five five14 four gizero minate nine nterm1 nterm2 nterm3 nterm4 nterm5 nterm6 one one024 one76 onebpi piby4 rtpiin seven seven2 sinz t temp three twelhu twent8 x xcube xhigh1 xhigh2 xhigh3 xlow1 xminus zero zeta 

if isempty(four), four = real( 4.0d+00); end;
if isempty(nterm1), nterm1 = 28; end;
if isempty(nterm2), nterm2 = 23; end;
if isempty(nterm3), nterm3 = 39; end;
if isempty(one), one = real( 1.0); end;
if isempty(three), three = 3.0d+00; end;
if isempty(x), x=0; end;
if isempty(zero), zero = 0.0d+00; end;
if isempty(bi), bi=0; end;
if isempty(cheb1), cheb1=0; end;
if isempty(cheb2), cheb2=0; end;
if isempty(cosz), cosz=0; end;
if isempty(sinz), sinz=0; end;
if isempty(t), t=0; end;
if isempty(temp), temp=0; end;
if isempty(xcube), xcube=0; end;
if isempty(xminus), xminus=0; end;
if isempty(zeta), zeta=0; end;
try;argip1(0+1);catch;   argip1(0+1)=[0.26585770795022745082d0];  end;
try;argip1(1+1);catch;   argip1(1+1)=[-0.10500333097501922907d0];  end;
try;argip1(2+1);catch;   argip1(2+1)=[0.841347475328454492d-2];  end;
try;argip1(3+1);catch;   argip1(3+1)=[0.2021067387813439541d-1];  end;
try;argip1(4+1);catch;   argip1(4+1)=[-0.1559576113863552234d-1];  end;
try;argip1(5+1);catch;   argip1(5+1)=[0.564342939043256481d-2];  end;
try;argip1(6+1);catch;   argip1(6+1)=[-0.59776844826655809d-3];  end;
try;argip1(7+1);catch;   argip1(7+1)=[-0.42833850264867728d-3];  end;
try;argip1(8+1);catch;   argip1(8+1)=[0.22605662380909027d-3];  end;
try;argip1(9+1);catch;   argip1(9+1)=[-0.3608332945592260d-4];  end;
try;argip1(10+1);catch;   argip1(10+1)=[-0.785518988788901d-5];  end;
try;argip1(11+1);catch;   argip1(11+1)=[0.473252480746370d-5];  end;
try;argip1(12+1);catch;   argip1(12+1)=[-0.59743513977694d-6];  end;
try;argip1(13+1);catch;   argip1(13+1)=[-0.15917609165602d-6];  end;
try;argip1(14+1);catch;   argip1(14+1)=[0.6336129065570d-7];  end;
try;argip1(15+1);catch;   argip1(15+1)=[-0.276090232648d-8];  end;
try;argip1(16+1);catch;   argip1(16+1)=[-0.256064154085d-8];  end;
try;argip1(17+1);catch;   argip1(17+1)=[0.47798676856d-9];  end;
try;argip1(18+1);catch;   argip1(18+1)=[0.4488131863d-10];  end;
try;argip1(19+1);catch;   argip1(19+1)=[-0.2346508882d-10];  end;
try;argip1(20+1);catch;   argip1(20+1)=[0.76839085d-12];  end;
try;argip1(21+1);catch;   argip1(21+1)=[0.73227985d-12];  end;
try;argip1(22+1);catch;   argip1(22+1)=[-0.8513687d-13];  end;
try;argip1(23+1);catch;   argip1(23+1)=[-0.1630201d-13];  end;
try;argip1(24+1);catch;   argip1(24+1)=[0.356769d-14];  end;
try;argip1(25+1);catch;   argip1(25+1)=[0.25001d-15];  end;
try;argip1(26+1);catch;   argip1(26+1)=[-0.10859d-15];  end;
try;argip1(27+1);catch;   argip1(27+1)=[-0.158d-17];  end;
try;argip1(28+1);catch;   argip1(28+1)=[0.275d-17];  end;
try;argip1(29+1);catch;   argip1(29+1)=[-0.5d-19];  end;
try;argip1(30+1);catch;   argip1(30+1)=[-0.6d-19];  end;
try;argip2(0+1);catch;   argip2(0+1)=[2.00473712275801486391d0];  end;
try;argip2(1+1);catch;   argip2(1+1)=[0.294184139364406724d-2];  end;
try;argip2(2+1);catch;   argip2(2+1)=[0.71369249006340167d-3];  end;
try;argip2(3+1);catch;   argip2(3+1)=[0.17526563430502267d-3];  end;
try;argip2(4+1);catch;   argip2(4+1)=[0.4359182094029882d-4];  end;
try;argip2(5+1);catch;   argip2(5+1)=[0.1092626947604307d-4];  end;
try;argip2(6+1);catch;   argip2(6+1)=[0.272382418399029d-5];  end;
try;argip2(7+1);catch;   argip2(7+1)=[0.66230900947687d-6];  end;
try;argip2(8+1);catch;   argip2(8+1)=[0.15425323370315d-6];  end;
try;argip2(9+1);catch;   argip2(9+1)=[0.3418465242306d-7];  end;
try;argip2(10+1);catch;   argip2(10+1)=[0.728157724894d-8];  end;
try;argip2(11+1);catch;   argip2(11+1)=[0.151588525452d-8];  end;
try;argip2(12+1);catch;   argip2(12+1)=[0.30940048039d-9];  end;
try;argip2(13+1);catch;   argip2(13+1)=[0.6149672614d-10];  end;
try;argip2(14+1);catch;   argip2(14+1)=[0.1202877045d-10];  end;
try;argip2(15+1);catch;   argip2(15+1)=[0.233690586d-11];  end;
try;argip2(16+1);catch;   argip2(16+1)=[0.43778068d-12];  end;
try;argip2(17+1);catch;   argip2(17+1)=[0.7996447d-13];  end;
try;argip2(18+1);catch;   argip2(18+1)=[0.1494075d-13];  end;
try;argip2(19+1);catch;   argip2(19+1)=[0.246790d-14];  end;
try;argip2(20+1);catch;   argip2(20+1)=[0.37672d-15];  end;
try;argip2(21+1);catch;   argip2(21+1)=[0.7701d-16];  end;
try;argip2(22+1);catch;   argip2(22+1)=[0.354d-17];  end;
try;argip2(23+1);catch;   argip2(23+1)=[-0.49d-18];  end;
try;argip2(24+1);catch;   argip2(24+1)=[0.62d-18];  end;
try;argip2(25+1);catch;   argip2(25+1)=[-0.40d-18];  end;
try;argip2(26+1);catch;   argip2(26+1)=[-0.1d-19];  end;
try;argip2(27+1);catch;   argip2(27+1)=[0.2d-19];  end;
try;argip2(28+1);catch;   argip2(28+1)=[-0.3d-19];  end;
try;argip2(29+1);catch;   argip2(29+1)=[0.1d-19];  end;
try;argin1(0+1);catch;   argin1(0+1)=[-0.20118965056732089130d0];  end;
try;argin1(1+1);catch;   argin1(1+1)=[-0.7244175303324530499d-1];  end;
try;argin1(2+1);catch;   argin1(2+1)=[0.4505018923894780120d-1];  end;
try;argin1(3+1);catch;   argin1(3+1)=[-0.24221371122078791099d0];  end;
try;argin1(4+1);catch;   argin1(4+1)=[0.2717884964361678294d-1];  end;
try;argin1(5+1);catch;   argin1(5+1)=[-0.5729321004818179697d-1];  end;
try;argin1(6+1);catch;   argin1(6+1)=[-0.18382107860337763587d0];  end;
try;argin1(7+1);catch;   argin1(7+1)=[0.7751546082149475511d-1];  end;
try;argin1(8+1);catch;   argin1(8+1)=[0.18386564733927560387d0];  end;
try;argin1(9+1);catch;   argin1(9+1)=[0.2921504250185567173d-1];  end;
try;argin1(10+1);catch;   argin1(10+1)=[-0.6142294846788018811d-1];  end;
try;argin1(11+1);catch;   argin1(11+1)=[-0.2999312505794616238d-1];  end;
try;argin1(12+1);catch;   argin1(12+1)=[0.585937118327706636d-2];  end;
try;argin1(13+1);catch;   argin1(13+1)=[0.822221658497402529d-2];  end;
try;argin1(14+1);catch;   argin1(14+1)=[0.132579817166846893d-2];  end;
try;argin1(15+1);catch;   argin1(15+1)=[-0.96248310766565126d-3];  end;
try;argin1(16+1);catch;   argin1(16+1)=[-0.45065515998211807d-3];  end;
try;argin1(17+1);catch;   argin1(17+1)=[0.772423474325474d-5];  end;
try;argin1(18+1);catch;   argin1(18+1)=[0.5481874134758052d-4];  end;
try;argin1(19+1);catch;   argin1(19+1)=[0.1245898039742876d-4];  end;
try;argin1(20+1);catch;   argin1(20+1)=[-0.246196891092083d-5];  end;
try;argin1(21+1);catch;   argin1(21+1)=[-0.169154183545285d-5];  end;
try;argin1(22+1);catch;   argin1(22+1)=[-0.16769153169442d-6];  end;
try;argin1(23+1);catch;   argin1(23+1)=[0.9636509337672d-7];  end;
try;argin1(24+1);catch;   argin1(24+1)=[0.3253314928030d-7];  end;
try;argin1(25+1);catch;   argin1(25+1)=[0.5091804231d-10];  end;
try;argin1(26+1);catch;   argin1(26+1)=[-0.209180453553d-8];  end;
try;argin1(27+1);catch;   argin1(27+1)=[-0.41237387870d-9];  end;
try;argin1(28+1);catch;   argin1(28+1)=[0.4163338253d-10];  end;
try;argin1(29+1);catch;   argin1(29+1)=[0.3032532117d-10];  end;
try;argin1(30+1);catch;   argin1(30+1)=[0.340580529d-11];  end;
try;argin1(31+1);catch;   argin1(31+1)=[-0.88444592d-12];  end;
try;argin1(32+1);catch;   argin1(32+1)=[-0.31639612d-12];  end;
try;argin1(33+1);catch;   argin1(33+1)=[-0.1505076d-13];  end;
try;argin1(34+1);catch;   argin1(34+1)=[0.1104148d-13];  end;
try;argin1(35+1);catch;   argin1(35+1)=[0.246508d-14];  end;
try;argin1(36+1);catch;   argin1(36+1)=[-0.3107d-16];  end;
try;argin1(37+1);catch;   argin1(37+1)=[-0.9851d-16];  end;
try;argin1(38+1);catch;   argin1(38+1)=[-0.1453d-16];  end;
try;argin1(39+1);catch;   argin1(39+1)=[0.118d-17];  end;
try;argin1(40+1);catch;   argin1(40+1)=[0.67d-18];  end;
try;argin1(41+1);catch;   argin1(41+1)=[0.6d-19];  end;
try;argin1(42+1);catch;   argin1(42+1)=[-0.1d-19];  end;
if isempty(arbin1),   arbin1=[1.99983763583586155980d0,-0.8104660923669418d-4,0.13475665984689d-6,-0.70855847143d-9,0.748184187d-11,-0.12902774d-12,0.322504d-14,-0.10809d-15,0.460d-17,-0.24d-18,0.1d-19];  end;
if isempty(arbin2),   arbin2=[0.13872356453879120276d0,-0.8239286225558228d-4,0.26720919509866d-6,-0.207423685368d-8,0.2873392593d-10,-0.60873521d-12,0.1792489d-13,-0.68760d-15,0.3280d-16,-0.188d-17,0.13d-18,-0.1d-19];  end;
if isempty(arhin1),   arhin1=[1.99647720399779650525d0,-0.187563779407173213d-2,-0.12186470897787339d-3,-0.814021609659287d-5,-0.55050925953537d-6,-0.3763008043303d-7,-0.258858362365d-8,-0.17931829265d-9,-0.1245916873d-10,-0.87171247d-12,-0.6084943d-13,-0.431178d-14,-0.29787d-15,-0.2210d-16,-0.136d-17,-0.14d-18];  end;
if isempty(five),   five=[5.0d0];  end;
if isempty(seven), seven=[7.0d0];  end;
if isempty(minate), minate=[-8.0d0];  end;
if isempty(nine),   nine=[9.0d0];  end;
if isempty(twent8), twent8=[28.0d0];  end;
if isempty(seven2), seven2=[72.0d0];  end;
if isempty(one76),   one76=[176.0d0];  end;
if isempty(five14), five14=[514.0d0];  end;
if isempty(one024),   one024=[1024.0d0];  end;
if isempty(twelhu), twelhu=[1200.0d0];  end;
if isempty(gizero),   gizero=[0.20497554248200024505d0];  end;
if isempty(onebpi),   onebpi=[0.31830988618379067154d0];  end;
if isempty(piby4),   piby4=[0.78539816339744830962d0];  end;
if isempty(rtpiin),   rtpiin=[0.56418958354775628695d0];  end;
%
%   Machine-dependent constants (suitable for IEEE machines)
%
if isempty(nterm4),   nterm4=[9];  end;
if isempty(nterm5), nterm5=[10];  end;
if isempty(nterm6), nterm6=[14];  end;
if isempty(xlow1),   xlow1=[2.22045d-16];  end;
if isempty(xhigh1), xhigh1=[208063.8307d0];  end;
if isempty(xhigh2),   xhigh2=[0.14274d308];  end;
if isempty(xhigh3), xhigh3=[-2097152.0d0];  end;
x = xvalue;
if( x < -xhigh1 .* xhigh1 ) ;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'AIRY_GI - Fatal error!');
writef(1,['%s','\n'], '  Argument too negative for accurate computation.');
airy_giResult = zero;
elseif( x <= xhigh3 ) ;
xminus = -x;
t = xminus .* sqrt( xminus );
zeta =( t + t ) ./ three;
temp = rtpiin ./ sqrt( sqrt( xminus ) );
cosz = cos( zeta + piby4 );
sinz = sin( zeta + piby4 ) ./ zeta;
xcube = x .* x .* x;
bi =( cosz + sinz .* five ./ seven2 ) .* temp;
t =( xcube + twelhu ) ./( one76 - xcube );
airy_giResult = bi + cheval( nterm6, arhin1, t ) .* onebpi ./ x;
elseif( x < minate ) ;
xminus = -x;
t = xminus .* sqrt( xminus );
zeta =( t + t ) ./ three;
temp = rtpiin ./ sqrt( sqrt( xminus ) );
cosz = cos( zeta + piby4 );
sinz = sin( zeta + piby4 ) ./ zeta;
xcube = x .* x .* x;
t = -( one024 ./( xcube ) + one );
[cheb1 , nterm4, arbin1, t ]=cheval( nterm4, arbin1, t );
[cheb2 , nterm5, arbin2, t ]=cheval( nterm5, arbin2, t );
bi =( cosz .* cheb1 + sinz .* cheb2 ) .* temp;
t =( xcube + twelhu ) ./( one76 - xcube );
airy_giResult = bi + cheval( nterm6, arhin1, t ) .* onebpi ./ x;
elseif( x <= -xlow1 ) ;
t = -( x + four ) ./ four;
[airy_giResult , nterm3, argin1, t ]=cheval( nterm3, argin1, t );
elseif( x < xlow1 ) ;
airy_giResult = gizero;
elseif( x <= seven ) ;
t =( nine .* x - twent8 ) ./( x + twent8 );
[airy_giResult , nterm1, argip1, t ]=cheval( nterm1, argip1, t );
elseif( x <= xhigh1 ) ;
xcube = x .* x .* x;
t =( twelhu - xcube ) ./( five14 + xcube );
airy_giResult = onebpi .* cheval( nterm2, argip2, t ) ./ x;
elseif( x <= xhigh2 ) ;
airy_giResult = onebpi ./ x;
else;
airy_giResult = zero;
end;
if ~isempty(inputname(1)), assignin('caller','FUntemp',xvalue); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end
function [n_data, x, fx]=airy_gi_values( n_data, x, fx ,varargin);
%*******************************************************************************
%
%! AIRY_GI_VALUES returns some values of the Airy Gi function.
%
%  Discussion:
%
%    The function is defined by:
%
%      AIRY_GI(x) = Integral ( 0 <= t < infinity ) sin ( x*t+t^3/3) dt / pi
%
%    The data was reported by McLeod.
%
%  Modified:
%
%    24 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real ( kind = 8 ) X, the argument of the function.
%
%    Output, real ( kind = 8 ) FX, the value of the function.
%
persistent fx_vec n_max x_vec 

if isempty(n_max), n_max = 20; end;
if isempty(fx_vec), fx_vec([1:n_max]) =[real(  0.20468308070040542435d+00),real(  0.18374662832557904078d+00),real( -0.11667221729601528265d+00),real(  0.31466934902729557596d+00),real( -0.37089040722426257729d+00),real( -0.25293059772424019694d+00),real(  0.28967410658692701936d+00),real( -0.34644836492634090590d+00),real(  0.28076035913873049496d+00),real(  0.21814994508094865815d+00),real(  0.20526679000810503329d+00),real(  0.22123695363784773258d+00),real(  0.23521843981043793760d+00),real(  0.82834303363768729338d-01),real(  0.45757385490989281893d-01),real(  0.44150012014605159922d-01),real(  0.39951133719508907541d-01),real(  0.35467706833949671483d-01),real(  0.31896005100679587981d-01),real(  0.26556892713512410405d-01) ]; end;
if isempty(x_vec), x_vec([1:n_max]) =[real(  -0.0019531250d+00),real(  -0.1250000000d+00),real(  -1.0000000000d+00),real(  -4.0000000000d+00),real(  -8.0000000000d+00),real(  -8.2500000000d+00),real(  -9.0000000000d+00),real( -10.0000000000d+00),real( -11.0000000000d+00),real( -13.0000000000d+00),real(   0.0019531250d+00),real(   0.1250000000d+00),real(   1.0000000000d+00),real(   4.0000000000d+00),real(   7.0000000000d+00),real(   7.2500000000d+00),real(   8.0000000000d+00),real(   9.0000000000d+00),real(  10.0000000000d+00),real(  12.0000000000d+00) ]; end;
if( n_data < 0 ) ;
n_data = 0;
end;
n_data = n_data + 1;
if( n_max < n_data ) ;
n_data = 0;
x = 0.0d+00;
fx = 0.0d+00;
else;
x  = x_vec(n_data);
fx = fx_vec(n_data);
end;
return;
end
function [airy_hiResult, xvalue ]=airy_hi( xvalue ,varargin);
%*******************************************************************************
%
%! AIRY_HI computes the modified Airy function Hi(x).
%
%  Discussion:
%
%    The function is defined by:
%
%      AIRY_HI(x) = Integral ( 0 <= t < infinity ) exp(x*t-t^3/3) dt / pi
%
%    The approximation uses Chebyshev expansions with the coefficients
%    given to 20 decimal places.
%
%    This subroutine is set up to work on IEEE machines.
%
%  Modified:
%
%    07 August 2004
%
%  Author:
%
%    Allan McLeod,
%    Department of Mathematics and Statistics,
%    Paisley University, High Street, Paisley, Scotland, PA12BE
%    macl_ms0@paisley.ac.uk
%
%  Reference:
%
%    Allan Mcleod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input, real ( kind = 8 ) XVALUE, the argument of the function.
%
%    Output, real ( kind = 8 ) AIRY_HI, the value of the function.
%
global realmlv
airy_hiResult=[];
persistent arbip argip1 arhin1 arhin2 arhip bi five14 four gi hizero lnrtpi minate nterm1 nterm2 nterm3 nterm4 nterm5 one one76 onebpi seven t temp thre43 three twelhu twelve two x xcube xhigh1 xlow1 xmax xneg1 xneg2 zero zeta 

if isempty(four), four = real( 4.0d+00); end;
if isempty(nterm1), nterm1 = 29; end;
if isempty(nterm2), nterm2 = 17; end;
if isempty(nterm3), nterm3 = 22; end;
if isempty(one), one = real( 1.0); end;
if isempty(three), three = 3.0d+00; end;
if isempty(two), two = 2.0d+00; end;
if isempty(x), x=0; end;
if isempty(zero), zero = 0.0d+00; end;
if isempty(bi), bi=0; end;
if isempty(gi), gi=0; end;
if isempty(t), t=0; end;
if isempty(temp), temp=0; end;
if isempty(xcube), xcube=0; end;
if isempty(zeta), zeta=0; end;
try;arhip(0+1);catch;   arhip(0+1)=[1.24013562561762831114d0];  end;
try;arhip(1+1);catch;   arhip(1+1)=[0.64856341973926535804d0];  end;
try;arhip(2+1);catch;   arhip(2+1)=[0.55236252592114903246d0];  end;
try;arhip(3+1);catch;   arhip(3+1)=[0.20975122073857566794d0];  end;
try;arhip(4+1);catch;   arhip(4+1)=[0.12025669118052373568d0];  end;
try;arhip(5+1);catch;   arhip(5+1)=[0.3768224931095393785d-1];  end;
try;arhip(6+1);catch;   arhip(6+1)=[0.1651088671548071651d-1];  end;
try;arhip(7+1);catch;   arhip(7+1)=[0.455922755211570993d-2];  end;
try;arhip(8+1);catch;   arhip(8+1)=[0.161828480477635013d-2];  end;
try;arhip(9+1);catch;   arhip(9+1)=[0.40841282508126663d-3];  end;
try;arhip(10+1);catch;   arhip(10+1)=[0.12196479721394051d-3];  end;
try;arhip(11+1);catch;   arhip(11+1)=[0.2865064098657610d-4];  end;
try;arhip(12+1);catch;   arhip(12+1)=[0.742221556424344d-5];  end;
try;arhip(13+1);catch;   arhip(13+1)=[0.163536231932831d-5];  end;
try;arhip(14+1);catch;   arhip(14+1)=[0.37713908188749d-6];  end;
try;arhip(15+1);catch;   arhip(15+1)=[0.7815800336008d-7];  end;
try;arhip(16+1);catch;   arhip(16+1)=[0.1638447121370d-7];  end;
try;arhip(17+1);catch;   arhip(17+1)=[0.319857665992d-8];  end;
try;arhip(18+1);catch;   arhip(18+1)=[0.61933905307d-9];  end;
try;arhip(19+1);catch;   arhip(19+1)=[0.11411161191d-9];  end;
try;arhip(20+1);catch;   arhip(20+1)=[0.2064923454d-10];  end;
try;arhip(21+1);catch;   arhip(21+1)=[0.360018664d-11];  end;
try;arhip(22+1);catch;   arhip(22+1)=[0.61401849d-12];  end;
try;arhip(23+1);catch;   arhip(23+1)=[0.10162125d-12];  end;
try;arhip(24+1);catch;   arhip(24+1)=[0.1643701d-13];  end;
try;arhip(25+1);catch;   arhip(25+1)=[0.259084d-14];  end;
try;arhip(26+1);catch;   arhip(26+1)=[0.39931d-15];  end;
try;arhip(27+1);catch;   arhip(27+1)=[0.6014d-16];  end;
try;arhip(28+1);catch;   arhip(28+1)=[0.886d-17];  end;
try;arhip(29+1);catch;   arhip(29+1)=[0.128d-17];  end;
try;arhip(30+1);catch;   arhip(30+1)=[0.18d-18];  end;
try;arhip(31+1);catch;   arhip(31+1)=[0.3d-19];  end;
try;arbip(0+1);catch;   arbip(0+1)=[2.00582138209759064905d0];  end;
try;arbip(1+1);catch;   arbip(1+1)=[0.294478449170441549d-2];  end;
try;arbip(2+1);catch;   arbip(2+1)=[0.3489754514775355d-4];  end;
try;arbip(3+1);catch;   arbip(3+1)=[0.83389733374343d-6];  end;
try;arbip(4+1);catch;   arbip(4+1)=[0.3136215471813d-7];  end;
try;arbip(5+1);catch;   arbip(5+1)=[0.167865306015d-8];  end;
try;arbip(6+1);catch;   arbip(6+1)=[0.12217934059d-9];  end;
try;arbip(7+1);catch;   arbip(7+1)=[0.1191584139d-10];  end;
try;arbip(8+1);catch;   arbip(8+1)=[0.154142553d-11];  end;
try;arbip(9+1);catch;   arbip(9+1)=[0.24844455d-12];  end;
try;arbip(10+1);catch;   arbip(10+1)=[0.4213012d-13];  end;
try;arbip(11+1);catch;   arbip(11+1)=[0.505293d-14];  end;
try;arbip(12+1);catch;   arbip(12+1)=[-0.60032d-15];  end;
try;arbip(13+1);catch;   arbip(13+1)=[-0.65474d-15];  end;
try;arbip(14+1);catch;   arbip(14+1)=[-0.22364d-15];  end;
try;arbip(15+1);catch;   arbip(15+1)=[-0.3015d-16];  end;
try;arbip(16+1);catch;   arbip(16+1)=[0.959d-17];  end;
try;arbip(17+1);catch;   arbip(17+1)=[0.616d-17];  end;
try;arbip(18+1);catch;   arbip(18+1)=[0.97d-18];  end;
try;arbip(19+1);catch;   arbip(19+1)=[-0.37d-18];  end;
try;arbip(20+1);catch;   arbip(20+1)=[-0.21d-18];  end;
try;arbip(21+1);catch;   arbip(21+1)=[-0.1d-19];  end;
try;arbip(22+1);catch;   arbip(22+1)=[0.2d-19];  end;
try;arbip(23+1);catch;   arbip(23+1)=[0.1d-19];  end;
try;argip1(0+1);catch;   argip1(0+1)=[2.00473712275801486391d0];  end;
try;argip1(1+1);catch;   argip1(1+1)=[0.294184139364406724d-2];  end;
try;argip1(2+1);catch;   argip1(2+1)=[0.71369249006340167d-3];  end;
try;argip1(3+1);catch;   argip1(3+1)=[0.17526563430502267d-3];  end;
try;argip1(4+1);catch;   argip1(4+1)=[0.4359182094029882d-4];  end;
try;argip1(5+1);catch;   argip1(5+1)=[0.1092626947604307d-4];  end;
try;argip1(6+1);catch;   argip1(6+1)=[0.272382418399029d-5];  end;
try;argip1(7+1);catch;   argip1(7+1)=[0.66230900947687d-6];  end;
try;argip1(8+1);catch;   argip1(8+1)=[0.15425323370315d-6];  end;
try;argip1(9+1);catch;   argip1(9+1)=[0.3418465242306d-7];  end;
try;argip1(10+1);catch;   argip1(10+1)=[0.728157724894d-8];  end;
try;argip1(11+1);catch;   argip1(11+1)=[0.151588525452d-8];  end;
try;argip1(12+1);catch;   argip1(12+1)=[0.30940048039d-9];  end;
try;argip1(13+1);catch;   argip1(13+1)=[0.6149672614d-10];  end;
try;argip1(14+1);catch;   argip1(14+1)=[0.1202877045d-10];  end;
try;argip1(15+1);catch;   argip1(15+1)=[0.233690586d-11];  end;
try;argip1(16+1);catch;   argip1(16+1)=[0.43778068d-12];  end;
try;argip1(17+1);catch;   argip1(17+1)=[0.7996447d-13];  end;
try;argip1(18+1);catch;   argip1(18+1)=[0.1494075d-13];  end;
try;argip1(19+1);catch;   argip1(19+1)=[0.246790d-14];  end;
try;argip1(20+1);catch;   argip1(20+1)=[0.37672d-15];  end;
try;argip1(21+1);catch;   argip1(21+1)=[0.7701d-16];  end;
try;argip1(22+1);catch;   argip1(22+1)=[0.354d-17];  end;
try;argip1(23+1);catch;   argip1(23+1)=[-0.49d-18];  end;
try;argip1(24+1);catch;   argip1(24+1)=[0.62d-18];  end;
try;argip1(25+1);catch;   argip1(25+1)=[-0.40d-18];  end;
try;argip1(26+1);catch;   argip1(26+1)=[-0.1d-19];  end;
try;argip1(27+1);catch;   argip1(27+1)=[0.2d-19];  end;
try;argip1(28+1);catch;   argip1(28+1)=[-0.3d-19];  end;
try;argip1(29+1);catch;   argip1(29+1)=[0.1d-19];  end;
try;arhin1(0+1);catch;   arhin1(0+1)=[0.31481017206423404116d0];  end;
try;arhin1(1+1);catch;   arhin1(1+1)=[-0.16414499216588964341d0];  end;
try;arhin1(2+1);catch;   arhin1(2+1)=[0.6176651597730913071d-1];  end;
try;arhin1(3+1);catch;   arhin1(3+1)=[-0.1971881185935933028d-1];  end;
try;arhin1(4+1);catch;   arhin1(4+1)=[0.536902830023331343d-2];  end;
try;arhin1(5+1);catch;   arhin1(5+1)=[-0.124977068439663038d-2];  end;
try;arhin1(6+1);catch;   arhin1(6+1)=[0.24835515596994933d-3];  end;
try;arhin1(7+1);catch;   arhin1(7+1)=[-0.4187024096746630d-4];  end;
try;arhin1(8+1);catch;   arhin1(8+1)=[0.590945437979124d-5];  end;
try;arhin1(9+1);catch;   arhin1(9+1)=[-0.68063541184345d-6];  end;
try;arhin1(10+1);catch;   arhin1(10+1)=[0.6072897629164d-7];  end;
try;arhin1(11+1);catch;   arhin1(11+1)=[-0.367130349242d-8];  end;
try;arhin1(12+1);catch;   arhin1(12+1)=[0.7078017552d-10];  end;
try;arhin1(13+1);catch;   arhin1(13+1)=[0.1187894334d-10];  end;
try;arhin1(14+1);catch;   arhin1(14+1)=[-0.120898723d-11];  end;
try;arhin1(15+1);catch;   arhin1(15+1)=[0.1189656d-13];  end;
try;arhin1(16+1);catch;   arhin1(16+1)=[0.594128d-14];  end;
try;arhin1(17+1);catch;   arhin1(17+1)=[-0.32257d-15];  end;
try;arhin1(18+1);catch;   arhin1(18+1)=[-0.2290d-16];  end;
try;arhin1(19+1);catch;   arhin1(19+1)=[0.253d-17];  end;
try;arhin1(20+1);catch;   arhin1(20+1)=[0.9d-19];  end;
try;arhin1(21+1);catch;   arhin1(21+1)=[-0.2d-19];  end;
if isempty(arhin2),   arhin2=[1.99647720399779650525d0,-0.187563779407173213d-2,-0.12186470897787339d-3,-0.814021609659287d-5,-0.55050925953537d-6,-0.3763008043303d-7,-0.258858362365d-8,-0.17931829265d-9,-0.1245916873d-10,-0.87171247d-12,-0.6084943d-13,-0.431178d-14,-0.29787d-15,-0.2210d-16,-0.136d-17,-0.14d-18];  end;
if isempty(seven),   seven=[7.0d0];  end;
if isempty(minate),   minate=[-8.0d0];  end;
if isempty(twelve), twelve=[12.0d0];  end;
if isempty(one76), one76=[176.0d0];  end;
if isempty(thre43),   thre43=[343.0d0];  end;
if isempty(five14), five14=[514.0d0];  end;
if isempty(twelhu), twelhu=[1200.0d0];  end;
if isempty(hizero),   hizero=[0.40995108496400049010d0];  end;
if isempty(lnrtpi),   lnrtpi=[0.57236494292470008707d0];  end;
if isempty(onebpi),   onebpi=[0.31830988618379067154d0];  end;
%
%   Machine-dependent constants (suitable for IEEE machines)
%
if isempty(nterm4),   nterm4=[19];  end;
if isempty(nterm5), nterm5=[14];  end;
if isempty(xlow1),   xlow1=[2.220446d-16];  end;
if isempty(xhigh1), xhigh1=[104.4175d0];  end;
if isempty(xneg1),   xneg1=[-0.14274d308];  end;
if isempty(xneg2), xneg2=[-208063.831d0];  end;
if isempty(xmax), xmax=[1.79d308];  end;
x = xvalue;
if( xhigh1 < x ) ;
writef(1,['%s','\n'], ' ');
writef(1,['%s','\n'], 'AIRY_HI - Fatal error!');
writef(1,['%s','\n'], '  Argument too large.');
airy_hiResult = xmax;
if ~isempty(inputname(1)), assignin('caller','FUntemp',xvalue); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
%  Code for x < 0.0
%
if( x < zero ) ;
if( x < minate ) ;
if( x < xneg1 ) ;
airy_hiResult = zero;
else;
if( x < xneg2 ) ;
temp = one;
airy_hiResult = - temp .* onebpi ./ x;
else;
xcube = x .* x .* x;
t =( xcube + twelhu ) ./( one76 - xcube );
[temp , nterm5, arhin2, t ]=cheval( nterm5, arhin2, t );
airy_hiResult = - temp .* onebpi ./ x;
end;
end;
else;
if( -xlow1 < x ) ;
airy_hiResult = hizero;
else;
t =( four .* x + twelve ) ./( x - twelve );
[airy_hiResult , nterm4, arhin1, t ]=cheval( nterm4, arhin1, t );
end;
end;
%
%   Code for x >= 0.0
%
else;
if( x <= seven ) ;
if( x < xlow1 ) ;
airy_hiResult = hizero;
else;
t =( x + x ) ./ seven - one;
temp =( x + x + x ) ./ two;
airy_hiResult = exp( temp ) .* cheval( nterm1, arhip, t );
end;
else;
xcube = x .* x .* x;
temp = sqrt( xcube );
zeta =( temp + temp ) ./ three;
t = two .*( sqrt( thre43 ./ xcube ) ) - one;
[temp , nterm2, arbip, t ]=cheval( nterm2, arbip, t );
temp = zeta + log( temp ) - log( x ) ./ four - lnrtpi;
bi = exp( temp );
t =( twelhu - xcube ) ./( xcube + five14 );
gi = cheval( nterm3, argip1, t ) .* onebpi ./ x;
airy_hiResult = bi - gi;
end;
end;
if ~isempty(inputname(1)), assignin('caller','FUntemp',xvalue); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end
function [n_data, x, fx]=airy_hi_values( n_data, x, fx ,varargin);
%*******************************************************************************
%
%! AIRY_HI_VALUES returns some values of the Airy Hi function.
%
%  Discussion:
%
%    The function is defined by:
%
%      AIRY_HI(x) = Integral ( 0 <= t < infinity ) exp ( x * t - t^3 / 3 ) dt / pi
%
%    The data was reported by McLeod.
%
%  Modified:
%
%    24 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real ( kind = 8 ) X, the argument of the function.
%
%    Output, real ( kind = 8 ) FX, the value of the function.
%
persistent fx_vec n_max x_vec 

if isempty(n_max), n_max = 20; end;
if isempty(fx_vec), fx_vec([1:n_max]) =[real( 0.40936798278458884024d+00),real( 0.37495291608048868619d+00),real( 0.22066960679295989454d+00),real( 0.77565356679703713590d-01),real( 0.39638826473124717315d-01),real( 0.38450072575004151871d-01),real( 0.35273216868317898556d-01),real( 0.31768535282502272742d-01),real( 0.28894408288051391369d-01),real( 0.24463284011678541180d-01),real( 0.41053540139998941517d+00),real( 0.44993502381204990817d+00),real( 0.97220515514243332184d+00),real( 0.83764237105104371193d+02),real( 0.80327744952044756016d+05),real( 0.15514138847749108298d+06),real( 0.11995859641733262114d+07),real( 0.21472868855967642259d+08),real( 0.45564115351632913590d+09),real( 0.32980722582904761929d+12) ]; end;
if isempty(x_vec), x_vec([1:n_max]) =[real(  -0.0019531250d+00),real(  -0.1250000000d+00),real(  -1.0000000000d+00),real(  -4.0000000000d+00),real(  -8.0000000000d+00),real(  -8.2500000000d+00),real(  -9.0000000000d+00),real( -10.0000000000d+00),real( -11.0000000000d+00),real( -13.0000000000d+00),real(   0.0019531250d+00),real(   0.1250000000d+00),real(   1.0000000000d+00),real(   4.0000000000d+00),real(   7.0000000000d+00),real(   7.2500000000d+00),real(   8.0000000000d+00),real(   9.0000000000d+00),real(  10.0000000000d+00),real(  12.0000000000d+00) ]; end;
if( n_data < 0 ) ;
n_data = 0;
end;
n_data = n_data + 1;
if( n_max < n_data ) ;
n_data = 0;
x = 0.0d+00;
fx = 0.0d+00;
else;
x  = x_vec(n_data);
fx = fx_vec(n_data);
end;
return;
end
function [arctan_intResult, xvalue ]=arctan_int( xvalue ,varargin);
%*******************************************************************************
%
%! ARCTAN_INT calculates the inverse tangent integral.
%
%  Discussion:
%
%    The function is defined by:
%
%      ARCTAN_INT(x) = Integral ( 0 <= t <= x ) arctan ( t ) / t dt
%
%    The approximation uses Chebyshev series with the coefficients
%    given to an accuracy of 20D.
%
%    This subroutine is set up to work on IEEE machines.
%
%  Modified:
%
%    24 August 2004
%
%  Author:
%
%    Allan McLeod,
%    Department of Mathematics and Statistics,
%    Paisley University, High Street, Paisley, Scotland, PA12BE
%    macl_ms0@paisley.ac.uk
%
%  Reference:
%
%    Allan Mcleod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input, real ( kind = 8 ) XVALUE, the argument of the function.
%
%    Output, real ( kind = 8 ) ARCTAN_INT, the value of the function.
%
arctan_intResult=[];
persistent atnina half ind nterms one t twobpi x xlow xupper zero 

if isempty(atnina), atnina([0:22]+1) =[1.91040361296235937512d0,-0.4176351437656746940d-1,0.275392550786367434d-2,-0.25051809526248881d-3,0.2666981285121171d-4,-0.311890514107001d-5,0.38833853132249d-6,-0.5057274584964d-7,0.681225282949d-8,-0.94212561654d-9,0.13307878816d-9,-0.1912678075d-10,0.278912620d-11,-0.41174820d-12,0.6142987d-13,-0.924929d-14,0.140387d-14,-0.21460d-15,0.3301d-16,-0.511d-17,0.79d-18,-0.12d-18,0.2d-19 ]; end;
if isempty(half), half = 0.5d+00; end;
if isempty(ind), ind=0; end;
if isempty(nterms), nterms = 19; end;
if isempty(one), one = real( 1.0); end;
if isempty(t), t=0; end;
if isempty(twobpi), twobpi = real( 0.63661977236758134308d0); end;
if isempty(x), x=0; end;
if isempty(xlow), xlow = real( 7.4505806d-9); end;
if isempty(xupper), xupper = real( 4.5036d15); end;
if isempty(zero), zero = 0.0d+00; end;
ind = 1;
x = xvalue;
if( x < zero ) ;
x = -x;
ind = -1;
end;
if( x < xlow ) ;
arctan_intResult = x;
elseif( x <= one ) ;
t = x .* x;
t =( t - half ) +( t - half );
arctan_intResult = x .* cheval( nterms, atnina, t );
elseif( x <= xupper ) ;
t = one ./( x .* x );
t =( t - half ) +( t - half );
arctan_intResult = log( x ) ./ twobpi + cheval( nterms, atnina, t ) ./ x;
else;
arctan_intResult = log( x ) ./ twobpi;
end;
if( ind < 0 ) ;
arctan_intResult = -arctan_intResult;
end;
if ~isempty(inputname(1)), assignin('caller','FUntemp',xvalue); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end
function [n_data, x, fx]=arctan_int_values( n_data, x, fx ,varargin);
%*******************************************************************************
%
%! ARCTAN_INT_VALUES returns some values of the inverse tangent integral.
%
%  Discussion:
%
%    The function is defined by:
%
%      ARCTAN_INT(x) = Integral ( 0 <= t <= x ) arctan ( t ) / t dt
%
%    The data was reported by McLeod.
%
%  Modified:
%
%    25 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real ( kind = 8 ) X, the argument of the function.
%
%    Output, real ( kind = 8 ) FX, the value of the function.
%
persistent fx_vec n_max x_vec 

if isempty(n_max), n_max = 20; end;
if isempty(fx_vec), fx_vec([1:n_max]) =[real(  0.19531241721588483191d-02),real( -0.39062433772980711281d-02),real(  0.78124470192576499535d-02),real(  0.15624576181996527280d-01),real( -0.31246610349485401551d-01),real(  0.62472911335014397321d-01),real(  0.12478419717389654039d+00),real( -0.24830175098230686908d+00),real(  0.48722235829452235711d+00),real(  0.91596559417721901505d+00),real(  0.12749694484943800618d+01),real( -0.15760154034463234224d+01),real(  0.24258878412859089996d+01),real(  0.33911633326292997361d+01),real(  0.44176450919422186583d+01),real( -0.47556713749547247774d+01),real(  0.50961912150934111303d+01),real(  0.53759175735714876256d+01),real( -0.61649904785027487422d+01),real(  0.72437843013083534973d+01) ]; end;
if isempty(x_vec), x_vec([1:n_max]) =[real(   0.0019531250d+00),real(  -0.0039062500d+00),real(   0.0078125000d+00),real(   0.0156250000d+00),real(  -0.0312500000d+00),real(   0.0625000000d+00),real(   0.1250000000d+00),real(  -0.2500000000d+00),real(   0.5000000000d+00),real(   1.0000000000d+00),real(   1.5000000000d+00),real(  -2.0000000000d+00),real(   4.0000000000d+00),real(   8.0000000000d+00),real(  16.0000000000d+00),real( -20.0000000000d+00),real(  25.0000000000d+00),real(  30.0000000000d+00),real( -50.0000000000d+00),real( 100.0000000000d+00) ]; end;
if( n_data < 0 ) ;
n_data = 0;
end;
n_data = n_data + 1;
if( n_max < n_data ) ;
n_data = 0;
x = 0.0d+00;
fx = 0.0d+00;
else;
x  = x_vec(n_data);
fx = fx_vec(n_data);
end;
return;
end
function [bessel_i0_intResult, xvalue ]=bessel_i0_int( xvalue ,varargin);
%*******************************************************************************
%
%! BESSEL_I0_INT computes the integral of the modified Bessel function I0(X).
%
%  Discussion:
%
%    The function is defined by:
%
%      I0_INT(x) = Integral ( 0 <= t <= x ) I0(t) dt
%
%    The program uses Chebyshev expansions, the coefficients of
%    which are given to 20 decimal places.
%
%    This subroutine is set up to work on IEEE machines.
%
%  Modified:
%
%    29 August 2004
%
%  Author:
%
%    Allan McLeod,
%    Department of Mathematics and Statistics,
%    Paisley University, High Street, Paisley, Scotland, PA12BE
%    macl_ms0@paisley.ac.uk
%
%  Reference:
%
%    Allan Mcleod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input, real ( kind = 8 ) XVALUE, the argument of the function.
%
%    Output, real ( kind = 8 ) BESSEL_I0_INT, the value of the function.
%
global realmlv
bessel_i0_intResult=[];
persistent ari01 ari0a ateen half ind lnr2pi nterm1 nterm2 t temp thirt6 three x xhigh xlow zero 

if isempty(ateen), ateen = real( 18.0d+00); end;
if isempty(half), half = 0.5d+00; end;
if isempty(ind), ind=0; end;
if isempty(lnr2pi), lnr2pi =real( 0.91893853320467274178d0); end;
if isempty(nterm1), nterm1 = 25; end;
if isempty(nterm2), nterm2 = 27; end;
if isempty(t), t=0; end;
if isempty(temp), temp=0; end;
if isempty(thirt6), thirt6 = real( 36.0d+00); end;
if isempty(three), three = 3.0d+00; end;
if isempty(x), x=0; end;
if isempty(xhigh), xhigh = real( 713.758339d0); end;
if isempty(xlow), xlow = real( 0.5161914d-7); end;
if isempty(zero), zero = 0.0d+00; end;
try;ari01(0+1);catch;   ari01(0+1)=[0.41227906926781516801d0];  end;
try;ari01(1+1);catch;   ari01(1+1)=[-0.34336345150081519562d0];  end;
try;ari01(2+1);catch;   ari01(2+1)=[0.22667588715751242585d0];  end;
try;ari01(3+1);catch;   ari01(3+1)=[-0.12608164718742260032d0];  end;
try;ari01(4+1);catch;   ari01(4+1)=[0.6012484628777990271d-1];  end;
try;ari01(5+1);catch;   ari01(5+1)=[-0.2480120462913358248d-1];  end;
try;ari01(6+1);catch;   ari01(6+1)=[0.892773389565563897d-2];  end;
try;ari01(7+1);catch;   ari01(7+1)=[-0.283253729936696605d-2];  end;
try;ari01(8+1);catch;   ari01(8+1)=[0.79891339041712994d-3];  end;
try;ari01(9+1);catch;   ari01(9+1)=[-0.20053933660964890d-3];  end;
try;ari01(10+1);catch;   ari01(10+1)=[0.4416816783014313d-4];  end;
try;ari01(11+1);catch;   ari01(11+1)=[-0.822377042246068d-5];  end;
try;ari01(12+1);catch;   ari01(12+1)=[0.120059794219015d-5];  end;
try;ari01(13+1);catch;   ari01(13+1)=[-0.11350865004889d-6];  end;
try;ari01(14+1);catch;   ari01(14+1)=[0.69606014466d-9];  end;
try;ari01(15+1);catch;   ari01(15+1)=[0.180622772836d-8];  end;
try;ari01(16+1);catch;   ari01(16+1)=[-0.26039481370d-9];  end;
try;ari01(17+1);catch;   ari01(17+1)=[-0.166188103d-11];  end;
try;ari01(18+1);catch;   ari01(18+1)=[0.510500232d-11];  end;
try;ari01(19+1);catch;   ari01(19+1)=[-0.41515879d-12];  end;
try;ari01(20+1);catch;   ari01(20+1)=[-0.7368138d-13];  end;
try;ari01(21+1);catch;   ari01(21+1)=[0.1279323d-13];  end;
try;ari01(22+1);catch;   ari01(22+1)=[0.103247d-14];  end;
try;ari01(23+1);catch;   ari01(23+1)=[-0.30379d-15];  end;
try;ari01(24+1);catch;   ari01(24+1)=[-0.1789d-16];  end;
try;ari01(25+1);catch;   ari01(25+1)=[0.673d-17];  end;
try;ari01(26+1);catch;   ari01(26+1)=[0.44d-18];  end;
try;ari01(27+1);catch;   ari01(27+1)=[-0.14d-18];  end;
try;ari01(28+1);catch;   ari01(28+1)=[-0.1d-19];  end;
try;ari0a(0+1);catch;   ari0a(0+1)=[2.03739654571143287070d0];  end;
try;ari0a(1+1);catch;   ari0a(1+1)=[0.1917631647503310248d-1];  end;
try;ari0a(2+1);catch;   ari0a(2+1)=[0.49923334519288147d-3];  end;
try;ari0a(3+1);catch;   ari0a(3+1)=[0.2263187103659815d-4];  end;
try;ari0a(4+1);catch;   ari0a(4+1)=[0.158682108285561d-5];  end;
try;ari0a(5+1);catch;   ari0a(5+1)=[0.16507855636318d-6];  end;
try;ari0a(6+1);catch;   ari0a(6+1)=[0.2385058373640d-7];  end;
try;ari0a(7+1);catch;   ari0a(7+1)=[0.392985182304d-8];  end;
try;ari0a(8+1);catch;   ari0a(8+1)=[0.46042714199d-9];  end;
try;ari0a(9+1);catch;   ari0a(9+1)=[-0.7072558172d-10];  end;
try;ari0a(10+1);catch;   ari0a(10+1)=[-0.6747183961d-10];  end;
try;ari0a(11+1);catch;   ari0a(11+1)=[-0.2026962001d-10];  end;
try;ari0a(12+1);catch;   ari0a(12+1)=[-0.87320338d-12];  end;
try;ari0a(13+1);catch;   ari0a(13+1)=[0.175520014d-11];  end;
try;ari0a(14+1);catch;   ari0a(14+1)=[0.60383944d-12];  end;
try;ari0a(15+1);catch;   ari0a(15+1)=[-0.3977983d-13];  end;
try;ari0a(16+1);catch;   ari0a(16+1)=[-0.8049048d-13];  end;
try;ari0a(17+1);catch;   ari0a(17+1)=[-0.1158955d-13];  end;
try;ari0a(18+1);catch;   ari0a(18+1)=[0.827318d-14];  end;
try;ari0a(19+1);catch;   ari0a(19+1)=[0.282290d-14];  end;
try;ari0a(20+1);catch;   ari0a(20+1)=[-0.77667d-15];  end;
try;ari0a(21+1);catch;   ari0a(21+1)=[-0.48731d-15];  end;
try;ari0a(22+1);catch;   ari0a(22+1)=[0.7279d-16];  end;
try;ari0a(23+1);catch;   ari0a(23+1)=[0.7873d-16];  end;
try;ari0a(24+1);catch;   ari0a(24+1)=[-0.785d-17];  end;
try;ari0a(25+1);catch;   ari0a(25+1)=[-0.1281d-16];  end;
try;ari0a(26+1);catch;   ari0a(26+1)=[0.121d-17];  end;
try;ari0a(27+1);catch;   ari0a(27+1)=[0.214d-17];  end;
try;ari0a(28+1);catch;   ari0a(28+1)=[-0.27d-18];  end;
try;ari0a(29+1);catch;   ari0a(29+1)=[-0.36d-18];  end;
try;ari0a(30+1);catch;   ari0a(30+1)=[0.7d-19];  end;
try;ari0a(31+1);catch;   ari0a(31+1)=[0.6d-19];  end;
try;ari0a(32+1);catch;   ari0a(32+1)=[-0.2d-19];  end;
try;ari0a(33+1);catch;   ari0a(33+1)=[-0.1d-19];  end;
ind = 1;
x = xvalue;
if( x