subroutine abram0_values ( n_data, x, fx ) c*********************************************************************72 c cc abram0_values() returns some values of the Abramowitz0 function. c c Discussion: c c The function is defined by: c c ABRAM0(X) = integral ( 0 <= T < +oo ) exp ( -T * T - X / T ) dT c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 21 May 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.87377726306985360531D+00, & 0.84721859650456925922D+00, & 0.77288934483988301615D+00, & 0.59684345853450151603D+00, & 0.29871735283675888392D+00, & 0.15004596450516388138D+00, & 0.11114662419157955096D+00, & 0.83909567153151897766D-01, & 0.56552321717943417515D-01, & 0.49876496603033790206D-01, & 0.44100889219762791328D-01, & 0.19738535180254062496D-01, & 0.86193088287161479900D-02, & 0.40224788162540127227D-02, & 0.19718658458164884826D-02, & 0.10045868340133538505D-02, & 0.15726917263304498649D-03, & 0.10352666912350263437D-04, & 0.91229759190956745069D-06, & 0.25628287737952698742D-09 / data x_vec / & 0.0019531250D+00, & 0.0078125000D+00, & 0.0312500000D+00, & 0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 1.2500000000D+00, & 1.5000000000D+00, & 1.8750000000D+00, & 2.0000000000D+00, & 2.1250000000D+00, & 3.0000000000D+00, & 4.0000000000D+00, & 5.0000000000D+00, & 6.0000000000D+00, & 7.0000000000D+00, & 10.0000000000D+00, & 15.0000000000D+00, & 20.0000000000D+00, & 40.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine abram1_values ( n_data, x, fx ) c*********************************************************************72 c cc abram1_values() returns some values of the Abramowitz1 function. c c Discussion: c c The function is defined by: c c ABRAM1(x) = integral ( 0 <= t < +oo ) t * exp ( -t^2 - x / t ) dt c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 21 August 2004 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.49828219848799921792D+00, & 0.49324391773047288556D+00, & 0.47431612784691234649D+00, & 0.41095983258760410149D+00, & 0.25317617388227035867D+00, & 0.14656338138597777543D+00, & 0.11421547056018366587D+00, & 0.90026307383483764795D-01, & 0.64088214170742303375D-01, & 0.57446614314166191085D-01, & 0.51581624564800730959D-01, & 0.25263719555776416016D-01, & 0.11930803330196594536D-01, & 0.59270542280915272465D-02, & 0.30609215358017829567D-02, & 0.16307382136979552833D-02, & 0.28371851916959455295D-03, & 0.21122150121323238154D-04, & 0.20344578892601627337D-05, & 0.71116517236209642290D-09 / data x_vec / & 0.0019531250D+00, & 0.0078125000D+00, & 0.0312500000D+00, & 0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 1.2500000000D+00, & 1.5000000000D+00, & 1.8750000000D+00, & 2.0000000000D+00, & 2.1250000000D+00, & 3.0000000000D+00, & 4.0000000000D+00, & 5.0000000000D+00, & 6.0000000000D+00, & 7.0000000000D+00, & 10.0000000000D+00, & 15.0000000000D+00, & 20.0000000000D+00, & 40.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine abram2_values ( n_data, x, fx ) c*********************************************************************72 c cc abram2_values() returns some values of the Abramowitz2 function. c c Discussion: c c The function is defined by: c c ABRAM2(x) = Integral ( 0 <= t < +oo ) t^2 * exp( -t^2 - x / t ) dt c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 February 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.44213858162107913430D+00, & 0.43923379545684026308D+00, & 0.42789857297092602234D+00, & 0.38652825661854504406D+00, & 0.26538204413231368110D+00, & 0.16848734838334595000D+00, & 0.13609200032513227112D+00, & 0.11070330027727917352D+00, & 0.82126019995530382267D-01, & 0.74538781999594581763D-01, & 0.67732034377612811390D-01, & 0.35641808698811851022D-01, & 0.17956589956618269083D-01, & 0.94058737143575370625D-02, & 0.50809356204299213556D-02, & 0.28149565414209719359D-02, & 0.53808696422559303431D-03, & 0.44821756380146327259D-04, & 0.46890678427324100410D-05, & 0.20161544850996420504D-08 / data x_vec / & 0.0019531250D+00, & 0.0078125000D+00, & 0.0312500000D+00, & 0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 1.2500000000D+00, & 1.5000000000D+00, & 1.8750000000D+00, & 2.0000000000D+00, & 2.1250000000D+00, & 3.0000000000D+00, & 4.0000000000D+00, & 5.0000000000D+00, & 6.0000000000D+00, & 7.0000000000D+00, & 10.0000000000D+00, & 15.0000000000D+00, & 20.0000000000D+00, & 40.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine agm_values ( n_data, a, b, fx ) c*********************************************************************72 c cc agm_values() returns some values of the arithmetic geometric mean. c c Discussion: c c The AGM is defined for nonnegative A and B. c c The AGM of numbers A and B is defined by setting c c A(0) = A, c B(0) = B c c A(N+1) = ( A(N) + B(N) ) / 2 c B(N+1) = sqrt ( A(N) * B(N) ) c c The two sequences both converge to AGM(A,B). c c In Mathematica, the AGM can be evaluated by c c ArithmeticGeometricMean [ a, b ] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 February 2008 c c Author: c c John Burkardt c c Reference: c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision A, B, the numbers whose AGM is desired. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 15 ) double precision a double precision a_vec(n_max) double precision b double precision b_vec(n_max) double precision fx double precision fx_vec(n_max) integer n_data save a_vec save b_vec save fx_vec data a_vec / & 22.0D+00, & 83.0D+00, & 42.0D+00, & 26.0D+00, & 4.0D+00, & 6.0D+00, & 40.0D+00, & 80.0D+00, & 90.0D+00, & 9.0D+00, & 53.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.5D+00 / data b_vec / & 96.0D+00, & 56.0D+00, & 7.0D+00, & 11.0D+00, & 63.0D+00, & 45.0D+00, & 75.0D+00, & 0.0D+00, & 35.0D+00, & 1.0D+00, & 53.0D+00, & 2.0D+00, & 4.0D+00, & 8.0D+00, & 8.0D+00 / data fx_vec / & 52.274641198704240049D+00, & 68.836530059858524345D+00, & 20.659301196734009322D+00, & 17.696854873743648823D+00, & 23.867049721753300163D+00, & 20.717015982805991662D+00, & 56.127842255616681863D+00, & 0.000000000000000000D+00, & 59.269565081229636528D+00, & 3.9362355036495554780D+00, & 53.000000000000000000D+00, & 1.4567910310469068692D+00, & 2.2430285802876025701D+00, & 3.6157561775973627487D+00, & 4.0816924080221632670D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 a = 0.0D+00 b = 0.0D+00 fx = 0.0D+00 else a = a_vec(n_data) b = b_vec(n_data) fx = fx_vec(n_data) end if return end subroutine airy_ai_values ( n_data, x, ai ) c*********************************************************************72 c cc airy_ai_values() returns some values of the Airy Ai(x) function. c c Discussion: c c The Airy functions Ai(X) and Bi(X) are a pair of linearly independent c solutions of the differential equation: c c W'' - X * W = 0 c c In Mathematica, the function can be evaluated by: c c AiryAi[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 February 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision AI, the value of the Airy AI function. c implicit none integer n_max parameter ( n_max = 11 ) double precision ai double precision ai_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save ai_vec save x_vec data ai_vec / & 0.3550280538878172D+00, & 0.3292031299435381D+00, & 0.3037031542863820D+00, & 0.2788064819550049D+00, & 0.2547423542956763D+00, & 0.2316936064808335D+00, & 0.2098000616663795D+00, & 0.1891624003981501D+00, & 0.1698463174443649D+00, & 0.1518868036405444D+00, & 0.1352924163128814D+00 / data x_vec / & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 ai = 0.0D+00 else x = x_vec(n_data) ai = ai_vec(n_data) end if return end subroutine airy_ai_int_values ( n_data, x, fx ) c*********************************************************************72 c cc airy_ai_int_values() returns some values of the integral of the Airy function. c c Discussion: c c The function is defined by: c c AIRY_AI_INT(x) = Integral ( 0 <= t <= x ) Ai(t) dt c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 02 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & -0.75228838916610124300D+00, & -0.57348350185854889466D+00, & -0.76569840313421291743D+00, & -0.65181015505382467421D+00, & -0.55881974894471876922D+00, & -0.56902352870716815309D+00, & -0.47800749642926168100D+00, & -0.46567398346706861416D+00, & -0.96783140945618013679D-01, & -0.34683049857035607494D-03, & 0.34658366917927930790D-03, & 0.27657581846051227124D-02, & 0.14595330491185717833D+00, & 0.23631734191710977960D+00, & 0.33289264538612212697D+00, & 0.33318759129779422976D+00, & 0.33332945170523851439D+00, & 0.33333331724248357420D+00, & 0.33333333329916901594D+00, & 0.33333333333329380187D+00 / data x_vec / & -12.0000000000D+00, & -11.0000000000D+00, & -10.0000000000D+00, & -9.5000000000D+00, & -9.0000000000D+00, & -6.5000000000D+00, & -4.0000000000D+00, & -1.0000000000D+00, & -0.2500000000D+00, & -0.0009765625D+00, & 0.0009765625D+00, & 0.0078125000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 4.0000000000D+00, & 4.5000000000D+00, & 6.0000000000D+00, & 8.0000000000D+00, & 10.0000000000D+00, & 12.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine airy_ai_prime_values ( n_data, x, aip ) c*********************************************************************72 c cc airy_ai_prime_values() returns some values of the Airy function Ai'(x). c c Discussion: c c The Airy functions Ai(X) and Bi(X) are a pair of linearly independent c solutions of the differential equation: c c W'' - X * W = 0 c c In Mathematica, the function can be evaluated by: c c AiryAiPrime[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 06 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision AIP, the derivative of the Airy AI function. c implicit none integer n_max parameter ( n_max = 11 ) double precision aip double precision aip_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save aip_vec save x_vec data aip_vec / & -0.2588194037928068D+00, & -0.2571304219075862D+00, & -0.2524054702856195D+00, & -0.2451463642190548D+00, & -0.2358320344192082D+00, & -0.2249105326646839D+00, & -0.2127932593891585D+00, & -0.1998511915822805D+00, & -0.1864128638072717D+00, & -0.1727638434616347D+00, & -0.1591474412967932D+00 / data x_vec / & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 aip = 0.0D+00 else x = x_vec(n_data) aip = aip_vec(n_data) end if return end subroutine airy_bi_values ( n_data, x, bi ) c*********************************************************************72 c cc airy_bi_values() returns some values of the Airy Bi(x) function. c c Discussion: c c The Airy functions Ai(X) and Bi(X) are a pair of linearly independent c solutions of the differential equation: c c W'' - X * W = 0 c c In Mathematica, the function can be evaluated by: c c AiryBi[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 06 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision BI, the value of the Airy BI function. c implicit none integer n_max parameter ( n_max = 11 ) double precision bi double precision bi_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save bi_vec save x_vec data bi_vec / & 0.6149266274460007D+00, & 0.6598616901941892D+00, & 0.7054642029186612D+00, & 0.7524855850873156D+00, & 0.8017730000135972D+00, & 0.8542770431031555D+00, & 0.9110633416949405D+00, & 0.9733286558781659D+00, & 0.1042422171231561D+01, & 0.1119872813134447D+01, & 0.1207423594952871D+01 / data x_vec / & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 bi = 0.0D+00 else x = x_vec(n_data) bi = bi_vec(n_data) end if return end subroutine airy_bi_int_values ( n_data, x, fx ) c*********************************************************************72 c cc airy_bi_int_values() returns some values of the integral of the Airy function. c c Discussion: c c The function is defined by: c c AIRY_BI_INT(x) = Integral ( 0 <= t <= x ) Bi(t) dt c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.17660819031554631869D-01, & -0.15040424806140020451D-01, & 0.14756446293227661920D-01, & -0.11847304264848446271D+00, & -0.64916741266165856037D-01, & 0.97260832464381044540D-01, & 0.50760058495287539119D-01, & -0.37300500963429492179D+00, & -0.13962988442666578531D+00, & -0.12001735266723296160D-02, & 0.12018836117890354598D-02, & 0.36533846550952011043D+00, & 0.87276911673800812196D+00, & 0.48219475263803429675D+02, & 0.44006525804904178439D+06, & 0.17608153976228301458D+07, & 0.73779211705220007228D+07, & 0.14780980310740671617D+09, & 0.97037614223613433849D+11, & 0.11632737638809878460D+15 / data x_vec / & -12.0000000000D+00, & -10.0000000000D+00, & -8.0000000000D+00, & -7.5000000000D+00, & -7.0000000000D+00, & -6.5000000000D+00, & -4.0000000000D+00, & -1.0000000000D+00, & -0.2500000000D+00, & -0.0019531250D+00, & 0.0019531250D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 4.0000000000D+00, & 8.0000000000D+00, & 8.5000000000D+00, & 9.0000000000D+00, & 10.0000000000D+00, & 12.0000000000D+00, & 14.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine airy_bi_prime_values ( n_data, x, bip ) c*********************************************************************72 c cc airy_bi_prime_values() returns some values of the Airy function Bi'(x). c c Discussion: c c The Airy functions Ai(X) and Bi(X) are a pair of linearly independent c solutions of the differential equation: c c W'' - X * W = 0 c c In Mathematica, the function can be evaluated by: c c AiryBiPrime[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision BIP, the derivative of the Airy BI function. c implicit none integer n_max parameter ( n_max = 11 ) double precision bip double precision bip_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save bip_vec save x_vec data bip_vec / & 0.4482883573538264D+00, & 0.4515126311496465D+00, & 0.4617892843621509D+00, & 0.4800490287524480D+00, & 0.5072816760506224D+00, & 0.5445725641405923D+00, & 0.5931444786342857D+00, & 0.6544059191721400D+00, & 0.7300069016152518D+00, & 0.8219038903072090D+00, & 0.9324359333927756D+00 / data x_vec / & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 bip = 0.0D+00 else x = x_vec(n_data) bip = bip_vec(n_data) end if return end subroutine airy_cai_values ( n_data, x, cai ) c*********************************************************************72 c cc airy_cai_values() returns some values of the Airy Ai(x) with complex argument. c c Discussion: c c The Airy functions Ai(X) and Bi(X) are a pair of linearly independent c solutions of the differential equation: c c W'' - X * W = 0 c c In Mathematica, the function can be evaluated by: c c AiryAi[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 13 April 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 c before the first call. On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double complex X, the argument of the function. c c Output, double complex CAI, the value of the Airy AI function. c implicit none integer n_max parameter ( n_max = 10 ) double complex cai double complex cai_vec(n_max) integer n_data double complex x double complex x_vec(n_max) save cai_vec save x_vec data cai_vec / & (0.1352924163128814D+00, 0.0000000000000000D+00), & (0.1433824486882056D+00, -0.1092193342707378D+00), & (0.2215404472324631D+00, -0.2588711788891803D+00), & (0.4763929771766866D+00, -0.3036484220291284D+00), & (0.5983692170633874D+00, -0.08154602160771214D+00), & (0.5355608832923521D+00, 0.00000000000000000D+00), & (0.5983692170633874D+00, 0.08154602160771214D+00), & (0.4763929771766866D+00, 0.3036484220291284D+00), & (0.2215404472324631D+00, 0.2588711788891803D+00), & (0.1433824486882056D+00, 0.1092193342707378D+00) / data x_vec / & ( 1.000000000000000D+00, 0.0000000000000000D+00), & ( 0.8090169943749474D+00, 0.5877852522924731D+00), & ( 0.3090169943749474D+00, 0.9510565162951536D+00), & (-0.3090169943749474D+00, 0.9510565162951536D+00), & (-0.8090169943749474D+00, 0.5877852522924731D+00), & (-1.0000000000000000D+00, 0.0000000000000000D+00), & (-0.8090169943749474D+00, -0.5877852522924731D+00), & (-0.3090169943749474D+00, -0.9510565162951536D+00), & ( 0.3090169943749474D+00, -0.9510565162951536D+00), & ( 0.8090169943749474D+00, -0.5877852522924731D+00) / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = ( 0.0D+00, 0.0D+00 ) cai = ( 0.0D+00, 0.0D+00 ) else x = x_vec(n_data) cai = cai_vec(n_data) end if return end subroutine airy_cbi_values ( n_data, x, cbi ) c*********************************************************************72 c cc airy_cbi_values() returns some values of the Airy Bi(x) with complex argument. c c Discussion: c c The Airy functions Ai(X) and Bi(X) are a pair of linearly independent c solutions of the differential equation: c c W'' - X * W = 0 c c In Mathematica, the function can be evaluated by: c c AiryBi[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 13 April 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 c before the first call. On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double complex X, the argument of the function. c c Output, double complex CBI, the value of the Airy BI function. c implicit none integer n_max parameter ( n_max = 10 ) double complex cbi double complex cbi_vec(n_max) integer n_data double complex x double complex x_vec(n_max) save cbi_vec save x_vec data cbi_vec / & ( 1.207423594952871D+00, 0.0000000000000000D+00), & ( 0.9127160108293936D+00, 0.3800456133135556D+00), & ( 0.6824453575635721D+00, 0.3343047153635002D+00), & ( 0.5726265660086474D+00, 0.3988641086982559D+00), & ( 0.2511841251049547D+00, 0.3401447690712719D+00), & ( 0.1039973894969446D+00, 0.0000000000000000D+00), & ( 0.2511841251049547D+00, -0.3401447690712719D+00), & ( 0.5726265660086474D+00, -0.3988641086982559D+00), & ( 0.6824453575635721D+00, -0.3343047153635002D+00), & ( 0.9127160108293936D+00, -0.3800456133135556D+00) / data x_vec / & ( 1.000000000000000D+00, 0.0000000000000000D+00), & ( 0.8090169943749474D+00, 0.5877852522924731D+00), & ( 0.3090169943749474D+00, 0.9510565162951536D+00), & (-0.3090169943749474D+00, 0.9510565162951536D+00), & (-0.8090169943749474D+00, 0.5877852522924731D+00), & (-1.0000000000000000D+00, 0.0000000000000000D+00), & (-0.8090169943749474D+00, -0.5877852522924731D+00), & (-0.3090169943749474D+00, -0.9510565162951536D+00), & ( 0.3090169943749474D+00, -0.9510565162951536D+00), & ( 0.8090169943749474D+00, -0.5877852522924731D+00) / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = ( 0.0D+00, 0.0D+00 ) cbi = ( 0.0D+00, 0.0D+00 ) else x = x_vec(n_data) cbi = cbi_vec(n_data) end if return end subroutine airy_gi_values ( n_data, x, fx ) c*****************************************************************************80 c cc airy_gi_values() returns some values of the Airy Gi function. c c Discussion: c c The function is defined by: c c AIRY_GI(x) = Integral ( 0 <= t < +oo ) sin ( x*t+t^3/3) dt / pi c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 11 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.20468308070040542435D+00, & 0.18374662832557904078D+00, & -0.11667221729601528265D+00, & 0.31466934902729557596D+00, & -0.37089040722426257729D+00, & -0.25293059772424019694D+00, & 0.28967410658692701936D+00, & -0.34644836492634090590D+00, & 0.28076035913873049496D+00, & 0.21814994508094865815D+00, & 0.20526679000810503329D+00, & 0.22123695363784773258D+00, & 0.23521843981043793760D+00, & 0.82834303363768729338D-01, & 0.45757385490989281893D-01, & 0.44150012014605159922D-01, & 0.39951133719508907541D-01, & 0.35467706833949671483D-01, & 0.31896005100679587981D-01, & 0.26556892713512410405D-01 / data x_vec / & -0.0019531250D+00, & -0.1250000000D+00, & -1.0000000000D+00, & -4.0000000000D+00, & -8.0000000000D+00, & -8.2500000000D+00, & -9.0000000000D+00, & -10.0000000000D+00, & -11.0000000000D+00, & -13.0000000000D+00, & 0.0019531250D+00, & 0.1250000000D+00, & 1.0000000000D+00, & 4.0000000000D+00, & 7.0000000000D+00, & 7.2500000000D+00, & 8.0000000000D+00, & 9.0000000000D+00, & 10.0000000000D+00, & 12.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine airy_hi_values ( n_data, x, fx ) c*****************************************************************************80 c cc airy_hi_values() returns some values of the Airy Hi function. c c Discussion: c c The function is defined by: c c AIRY_HI(x) = Integral ( 0 <= t < +oo ) exp(x*t-t^3/3) dt / pi c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 11 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.40936798278458884024D+00, & 0.37495291608048868619D+00, & 0.22066960679295989454D+00, & 0.77565356679703713590D-01, & 0.39638826473124717315D-01, & 0.38450072575004151871D-01, & 0.35273216868317898556D-01, & 0.31768535282502272742D-01, & 0.28894408288051391369D-01, & 0.24463284011678541180D-01, & 0.41053540139998941517D+00, & 0.44993502381204990817D+00, & 0.97220515514243332184D+00, & 0.83764237105104371193D+02, & 0.80327744952044756016D+05, & 0.15514138847749108298D+06, & 0.11995859641733262114D+07, & 0.21472868855967642259D+08, & 0.45564115351632913590D+09, & 0.32980722582904761929D+12 / data x_vec / & -0.0019531250D+00, & -0.1250000000D+00, & -1.0000000000D+00, & -4.0000000000D+00, & -8.0000000000D+00, & -8.2500000000D+00, & -9.0000000000D+00, & -10.0000000000D+00, & -11.0000000000D+00, & -13.0000000000D+00, & 0.0019531250D+00, & 0.1250000000D+00, & 1.0000000000D+00, & 4.0000000000D+00, & 7.0000000000D+00, & 7.2500000000D+00, & 8.0000000000D+00, & 9.0000000000D+00, & 10.0000000000D+00, & 12.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine arccos_values ( n_data, x, fx ) c*********************************************************************72 c cc arccos_values() returns some values of the arc cosine function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c ArcCos[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 10 January 2008 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 c before the first call. On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 12 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 1.6709637479564564156D+00, & 1.5707963267948966192D+00, & 1.4706289056333368229D+00, & 1.3694384060045658278D+00, & 1.2661036727794991113D+00, & 1.1592794807274085998D+00, & 1.0471975511965977462D+00, & 0.92729521800161223243D+00, & 0.79539883018414355549D+00, & 0.64350110879328438680D+00, & 0.45102681179626243254D+00, & 0.00000000000000000000D+00 / data x_vec / & -0.1D+00, & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine arccosh_values ( n_data, x, fx ) c*********************************************************************72 c cc arccosh_values() returns some values of the hyperbolic arc cosine function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c ArcCosh[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 10 January 2008 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 c before the first call. On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 15 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.0000000000000000000D+00, & 0.14130376948564857735D+00, & 0.44356825438511518913D+00, & 0.62236250371477866781D+00, & 0.75643291085695958624D+00, & 0.86701472649056510395D+00, & 0.96242365011920689500D+00, & 1.3169578969248167086D+00, & 1.7627471740390860505D+00, & 1.8115262724608531070D+00, & 2.0634370688955605467D+00, & 2.2924316695611776878D+00, & 2.9932228461263808979D+00, & 5.2982923656104845907D+00, & 7.6009022095419886114D+00 / data x_vec / & 1.0D+00, & 1.01D+00, & 1.1D+00, & 1.2D+00, & 1.3D+00, & 1.4D+00, & 1.5D+00, & 2.0D+00, & 3.0D+00, & 3.1415926535897932385D+00, & 4.0D+00, & 5.0D+00, & 10.0D+00, & 100.0D+00, & 1000.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine arcsin_values ( n_data, x, fx ) c*********************************************************************72 c cc arcsin_values() returns some values of the arc sine function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c ArcSin[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 10 January 2008 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 c before the first call. On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 12 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & -0.10016742116155979635D+00, & 0.00000000000000000000D+00, & 0.10016742116155979635D+00, & 0.20135792079033079146D+00, & 0.30469265401539750797D+00, & 0.41151684606748801938D+00, & 0.52359877559829887308D+00, & 0.64350110879328438680D+00, & 0.77539749661075306374D+00, & 0.92729521800161223243D+00, & 1.1197695149986341867D+00, & 1.5707963267948966192D+00 / data x_vec / & -0.1D+00, & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine arcsinh_values ( n_data, x, fx ) c*********************************************************************72 c cc arcsinh_values() returns some values of the hyperbolic arc sine function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c ArcSinh[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 10 January 2008 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 c before the first call. On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & -2.3124383412727526203D+00, & -0.88137358701954302523D+00, & 0.00000000000000000000D+00, & 0.099834078899207563327D+00, & 0.19869011034924140647D+00, & 0.29567304756342243910D+00, & 0.39003531977071527608D+00, & 0.48121182505960344750D+00, & 0.56882489873224753010D+00, & 0.65266656608235578681D+00, & 0.73266825604541086415D+00, & 0.80886693565278246251D+00, & 0.88137358701954302523D+00, & 1.4436354751788103425D+00, & 1.8184464592320668235D+00, & 2.0947125472611012942D+00, & 2.3124383412727526203D+00, & 2.9982229502979697388D+00, & 5.2983423656105887574D+00, & 7.6009027095419886115D+00 / data x_vec / & -5.0D+00, & -1.0D+00, & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 10.0D+00, & 100.0D+00, & 1000.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine arctan_values ( n_data, x, fx ) c*********************************************************************72 c cc arctan_values() returns some values of the arc tangent function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c ArcTan[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 10 January 2008 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 c before the first call. On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 11 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.00000000000000000000D+00, & 0.24497866312686415417D+00, & 0.32175055439664219340D+00, & 0.46364760900080611621D+00, & 0.78539816339744830962D+00, & 1.1071487177940905030D+00, & 1.2490457723982544258D+00, & 1.3258176636680324651D+00, & 1.3734007669450158609D+00, & 1.4711276743037345919D+00, & 1.5208379310729538578D+00 / data x_vec / & 0.00000000000000000000D+00, & 0.25000000000000000000D+00, & 0.33333333333333333333D+00, & 0.50000000000000000000D+00, & 1.0000000000000000000D+00, & 2.0000000000000000000D+00, & 3.0000000000000000000D+00, & 4.0000000000000000000D+00, & 5.0000000000000000000D+00, & 10.000000000000000000D+00, & 20.000000000000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine arctan2_values ( n_data, x, y, f ) c*********************************************************************72 c cc arctan2_values(): arc tangent function of two arguments. c c Discussion: c c In Mathematica, the function can be evaluated by: c c ArcTan[x,y] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 28 March 2010 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 c before the first call. On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, Y, the arguments of the function. c c Output, double precision F, the value of the function. c implicit none integer n_max parameter ( n_max = 19 ) double precision f double precision f_vec(n_max) integer n_data double precision x double precision x_vec(n_max) double precision y double precision y_vec(n_max) save f_vec save x_vec save y_vec data f_vec / & -1.5707963267948966192D+00, & -1.0471975511965977462D+00, & -0.52359877559829887308D+00, & 0.00000000000000000000D+00, & 0.52359877559829887308D+00, & 1.0471975511965977462D+00, & 1.5707963267948966192D+00, & 2.0943951023931954923D+00, & 2.6179938779914943654D+00, & 3.1415926535897932385D+00, & -2.6179938779914943654D+00, & -2.0943951023931954923D+00, & -1.5707963267948966192D+00, & -1.0471975511965977462D+00, & -0.52359877559829887308D+00, & 0.00000000000000000000D+00, & 0.52359877559829887308D+00, & 1.0471975511965977462D+00, & 1.5707963267948966192D+00 / data x_vec / & 0.00000000000000000000D+00, & 0.50000000000000000000D+00, & 0.86602540378443864676D+00, & 1.00000000000000000000D+00, & 0.86602540378443864676D+00, & 0.50000000000000000000D+00, & 0.00000000000000000000D+00, & -0.50000000000000000000D+00, & -0.86602540378443864676D+00, & -1.00000000000000000000D+00, & -0.86602540378443864676D+00, & -0.50000000000000000000D+00, & 0.00000000000000000000D+00, & 0.50000000000000000000D+00, & 0.86602540378443864676D+00, & 1.00000000000000000000D+00, & 0.86602540378443864676D+00, & 0.50000000000000000000D+00, & 0.00000000000000000000D+00 / data y_vec / & -1.00000000000000000000D+00, & -0.86602540378443864676D+00, & -0.50000000000000000000D+00, & 0.00000000000000000000D+00, & 0.50000000000000000000D+00, & 0.86602540378443864676D+00, & 1.00000000000000000000D+00, & 0.86602540378443864676D+00, & 0.50000000000000000000D+00, & 0.00000000000000000000D+00, & -0.50000000000000000000D+00, & -0.86602540378443864676D+00, & -1.00000000000000000000D+00, & -0.86602540378443864676D+00, & -0.50000000000000000000D+00, & 0.00000000000000000000D+00, & 0.50000000000000000000D+00, & 0.86602540378443864676D+00, & 1.00000000000000000000D+00/ if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 y = 0.0D+00 f = 0.0D+00 else x = x_vec(n_data) y = y_vec(n_data) f = f_vec(n_data) end if return end subroutine arctan_int_values ( n_data, x, fx ) c*********************************************************************72 c cc arctan_int_values() returns some values of the inverse tangent integral. c c Discussion: c c The function is defined by: c c ARCTAN_INT(x) = Integral ( 0 <= t <= x ) arctan ( t ) / t dt c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 18 March 2007 c c Author: c c John Burkardt c c Reference: c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.19531241721588483191D-02, & -0.39062433772980711281D-02, & 0.78124470192576499535D-02, & 0.15624576181996527280D-01, & -0.31246610349485401551D-01, & 0.62472911335014397321D-01, & 0.12478419717389654039D+00, & -0.24830175098230686908D+00, & 0.48722235829452235711D+00, & 0.91596559417721901505D+00, & 0.12749694484943800618D+01, & -0.15760154034463234224D+01, & 0.24258878412859089996D+01, & 0.33911633326292997361D+01, & 0.44176450919422186583D+01, & -0.47556713749547247774D+01, & 0.50961912150934111303D+01, & 0.53759175735714876256D+01, & -0.61649904785027487422D+01, & 0.72437843013083534973D+01 / data x_vec / & 0.0019531250D+00, & -0.0039062500D+00, & 0.0078125000D+00, & 0.0156250000D+00, & -0.0312500000D+00, & 0.0625000000D+00, & 0.1250000000D+00, & -0.2500000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 1.5000000000D+00, & -2.0000000000D+00, & 4.0000000000D+00, & 8.0000000000D+00, & 16.0000000000D+00, & -20.0000000000D+00, & 25.0000000000D+00, & 30.0000000000D+00, & -50.0000000000D+00, & 100.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine arctanh_values ( n_data, x, fx ) c*********************************************************************72 c cc arctanh_values() returns some values of the hyperbolic arc tangent function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c ArcTanh[x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 23 June 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 c before the first call. On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 15 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & -0.54930614433405484570D+00, & 0.00000000000000000000D+00, & 0.0010000003333335333335D+00, & 0.10033534773107558064D+00, & 0.20273255405408219099D+00, & 0.30951960420311171547D+00, & 0.42364893019360180686D+00, & 0.54930614433405484570D+00, & 0.69314718055994530942D+00, & 0.86730052769405319443D+00, & 1.0986122886681096914D+00, & 1.4722194895832202300D+00, & 2.6466524123622461977D+00, & 3.8002011672502000318D+00, & 7.2543286192620472067D+00 / data x_vec / & -0.5D+00, & 0.0D+00, & 0.001D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 0.99D+00, & 0.999D+00, & 0.999999D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bei0_values ( n_data, x, fx ) c*********************************************************************72 c cc bei0_values() returns some values of the Kelvin BEI function of order NU = 0. c c Discussion: c c The function is defined by: c c BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4)) c c where J(NU,X) is the J Bessel function. c c In Mathematica, BEI(NU,X) can be defined by: c c Im [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 28 June 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 11 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.0000000000000000D+00, & 0.06249321838219946D+00, & 0.2495660400366597D+00, & 0.5575600623030867D+00, & 0.9722916273066612D+00, & 1.457182044159804D+00, & 1.937586785266043D+00, & 2.283249966853915D+00, & 2.292690322699300D+00, & 1.686017203632139D+00, & 0.1160343815502004D+00 / data x_vec / & 0.0D+00, & 0.5D+00, & 1.0D+00, & 1.5D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 3.5D+00, & 4.0D+00, & 4.5D+00, & 5.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bei1_values ( n_data, x, fx ) c*********************************************************************72 c cc bei1_values() returns some values of the Kelvin BEI function of order NU = 1. c c Discussion: c c The function is defined by: c c BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4)) c c where J(NU,X) is the J Bessel function. c c In Mathematica, BEI(NU,X) can be defined by: c c Im [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 28 June 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 11 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.0000000000000000D+00, & 0.1711951797170153D+00, & 0.3075566313755366D+00, & 0.3678649890020899D+00, & 0.2997754370020335D+00, & 0.03866844396595048D+00, & -0.4874541770160708D+00, & -1.344042373111174D+00, & -2.563821688561078D+00, & -4.105685408400878D+00, & -5.797907901792625D+00 / data x_vec / & 0.0D+00, & 0.5D+00, & 1.0D+00, & 1.5D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 3.5D+00, & 4.0D+00, & 4.5D+00, & 5.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bell_values ( n_data, n, c ) c*********************************************************************72 c cc bell_values() returns some values of the Bell numbers for testing. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 03 January 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Parameters: c c Input/output, integer N_DATA. c On input, if N_DATA is 0, the first test data is returned, and N_DATA c is set to 1. On each subsequent call, the input value of N_DATA is c incremented and that test data item is returned, if available. When c there is no more test data, N_DATA is set to 0. c c Output, integer N, the order of the Bell number. c c Output, integer C, the value of the Bell number. c implicit none integer n_max parameter ( n_max = 11 ) integer c integer c_vec(n_max) integer n integer n_data integer n_vec(n_max) save c_vec save n_vec data c_vec / & 1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975 / data n_vec / & 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 n = 0 c = 0 else n = n_vec(n_data) c = c_vec(n_data) end if return end subroutine ber0_values ( n_data, x, fx ) c*********************************************************************72 c cc ber0_values() returns some values of the Kelvin BER function of order NU = 0. c c Discussion: c c The function is defined by: c c BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4)) c c where J(NU,X) is the J Bessel function. c c In Mathematica, BER(NU,X) can be defined by: c c Re [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 27 June 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 11 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 1.0000000000000000D+00, & 0.9990234639908383D+00, & 0.9843817812130869D+00, & 0.9210721835462558D+00, & 0.7517341827138082D+00, & 0.3999684171295313D+00, & -0.2213802495986939D+00, & -1.193598179589928D+00, & -2.563416557258580D+00, & -4.299086551599756D+00, & -6.230082478666358D+00 / data x_vec / & 0.0D+00, & 0.5D+00, & 1.0D+00, & 1.5D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 3.5D+00, & 4.0D+00, & 4.5D+00, & 5.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine ber1_values ( n_data, x, fx ) c*********************************************************************72 c cc ber1_values() returns some values of the Kelvin BER function of order NU = 1. c c Discussion: c c The function is defined by: c c BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4)) c c where J(NU,X) is the J Bessel function. c c In Mathematica, BER(NU,X) can be defined by: c c Re [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 28 June 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 11 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.0000000000000000D+00, & -0.1822431237551121D+00, & -0.3958682610197114D+00, & -0.6648654179597691D+00, & -0.9970776519264285D+00, & -1.373096897645111D+00, & -1.732644221128481D+00, & -1.959644131289749D+00, & -1.869248459031899D+00, & -1.202821631480086D+00, & 0.3597766667766728D+00 / data x_vec / & 0.0D+00, & 0.5D+00, & 1.0D+00, & 1.5D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 3.5D+00, & 4.0D+00, & 4.5D+00, & 5.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bernoulli_number_values ( n_data, n, c ) c*********************************************************************72 c cc bernoulli_number_values() returns some values of the Bernoulli numbers. c c Discussion: c c The Bernoulli numbers are rational. c c If we define the sum of the M-th powers of the first N integers as: c c SIGMA(M,N) = sum ( 0 <= I <= N ) I**M c c and let C(I,J) be the combinatorial coefficient: c c C(I,J) = Ic / ( ( I - J )c * Jc ) c c then the Bernoulli numbers B(J) satisfy: c c SIGMA(M,N) = 1/(M+1) * sum ( 0 <= J <= M ) C(M+1,J) B(J) * (N+1)**(M+1-J) c c In Mathematica, the function can be evaluated by: c c BernoulliB[n] c c First values: c c B0 1 = 1.00000000000 c B1 -1/2 = -0.50000000000 c B2 1/6 = 1.66666666666 c B3 0 = 0 c B4 -1/30 = -0.03333333333 c B5 0 = 0 c B6 1/42 = 0.02380952380 c B7 0 = 0 c B8 -1/30 = -0.03333333333 c B9 0 = 0 c B10 5/66 = 0.07575757575 c B11 0 = 0 c B12 -691/2730 = -0.25311355311 c B13 0 = 0 c B14 7/6 = 1.16666666666 c B15 0 = 0 c B16 -3617/510 = -7.09215686274 c B17 0 = 0 c B18 43867/798 = 54.97117794486 c B19 0 = 0 c B20 -174611/330 = -529.12424242424 c B21 0 = 0 c B22 854,513/138 = 6192.123 c B23 0 = 0 c B24 -236364091/2730 = -86580.257 c B25 0 = 0 c B26 8553103/6 = 1425517.16666 c B27 0 = 0 c B28 -23749461029/870 = -27298231.0678 c B29 0 = 0 c B30 8615841276005/14322 = 601580873.901 c c Recursion: c c With C(N+1,K) denoting the standard binomial coefficient, c c B(0) = 1.0 c B(N) = - ( sum ( 0 <= K .lt. N ) C(N+1,K) * B(K) ) / C(N+1,N) c c Special Values: c c Except for B(1), all Bernoulli numbers of odd index are 0. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 19 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer N, the order of the Bernoulli number. c c Output, double precision C, the value of the Bernoulli number. c implicit none integer n_max parameter ( n_max = 10 ) double precision c double precision c_vec(n_max) integer n integer n_data integer n_vec(n_max) save c_vec save n_vec data c_vec / & 0.1000000000000000D+01, & -0.5000000000000000D+00, & 0.1666666666666667D+00, & 0.0000000000000000D+00, & -0.3333333333333333D-01, & -0.2380952380952380D-01, & -0.3333333333333333D-01, & 0.7575757575757575D-01, & -0.5291242424242424D+03, & 0.6015808739006424D+09 / data n_vec / & 0, 1, 2, 3, 4, 6, 8, 10, 20, 30 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 n = 0 c = 0.0D+00 else n = n_vec(n_data) c = c_vec(n_data) end if return end subroutine bernoulli_poly_values ( n_data, n, x, fx ) c*********************************************************************72 c cc bernoulli_poly_values() returns some values of the Bernoulli polynomials. c c Discussion: c c In Mathematica, the function can be evaluated by: c c BernoulliB[n,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 19 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer N, the order of the Bernoulli polynomial. c c Output, double precision X, the argument of the Bernoulli polynomial. c c Output, double precision FX, the value of the Bernoulli polynomial. c implicit none integer n_max parameter ( n_max = 27 ) double precision fx double precision fx_vec(n_max) integer n integer n_data integer n_vec(n_max) double precision x double precision x_vec(n_max) save fx_vec save n_vec save x_vec data fx_vec / & 0.1000000000000000D+01, & -0.3000000000000000D+00, & 0.6666666666666667D-02, & 0.4800000000000000D-01, & -0.7733333333333333D-02, & -0.2368000000000000D-01, & 0.6913523809523810D-02, & 0.2490880000000000D-01, & -0.1014997333333333D-01, & -0.4527820800000000D-01, & 0.2332631815757576D-01, & -0.3125000000000000D+00, & -0.1142400000000000D+00, & -0.0176800000000000D+00, & 0.0156800000000000D+00, & 0.0147400000000000D+00, & 0.0000000000000000D+00, & -0.1524000000000000D-01, & -0.2368000000000000D-01, & -0.2282000000000000D-01, & -0.1376000000000000D-01, & 0.0000000000000000D+01, & 0.1376000000000000D-01, & 0.2282000000000000D-01, & 0.2368000000000000D-01, & 0.1524000000000000D-01, & 0.0000000000000000D+01 / data n_vec / & 0, & 1, & 2, & 3, & 4, & 5, & 6, & 7, & 8, & 9, & 10, & 5, & 5, & 5, & 5, & 5, & 5, & 5, & 5, & 5, & 5, & 5, & 5, & 5, & 5, & 5, & 5 / data x_vec / & 0.2D+00, & 0.2D+00, & 0.2D+00, & 0.2D+00, & 0.2D+00, & 0.2D+00, & 0.2D+00, & 0.2D+00, & 0.2D+00, & 0.2D+00, & 0.2D+00, & -0.5D+00, & -0.4D+00, & -0.3D+00, & -0.2D+00, & -0.1D+00, & 0.0D+00, & 0.1D+00, & 0.2D+00, & 0.3D+00, & 0.4D+00, & 0.5D+00, & 0.6D+00, & 0.7D+00, & 0.8D+00, & 0.9D+00, & 1.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 n = 0 x = 0.0D+00 fx = 0.0D+00 else n = n_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bernstein_poly_01_values ( n_data, n, k, x, b ) c*********************************************************************72 c cc bernstein_poly_01_values() returns some values of the Bernstein polynomials. c c Discussion: c c The Bernstein polynomials are assumed to be based on [0,1]. c c The formula for the Bernstein polynomials is c c B(N,I)(X) = [N!/(I!*(N-I)!)] * (1-X)^(N-I) * X^I c c In Mathematica, the function can be evaluated by: c c Binomial[n,i] * (1-x)^(n-i) * x^i c c First values: c c B(0,0)(X) = 1 c c B(1,0)(X) = 1-X c B(1,1)(X) = X c c B(2,0)(X) = (1-X)^2 c B(2,1)(X) = 2 * (1-X) * X c B(2,2)(X) = X^2 c c B(3,0)(X) = (1-X)^3 c B(3,1)(X) = 3 * (1-X)^2 * X c B(3,2)(X) = 3 * (1-X) * X^2 c B(3,3)(X) = X^3 c c B(4,0)(X) = (1-X)^4 c B(4,1)(X) = 4 * (1-X)^3 * X c B(4,2)(X) = 6 * (1-X)^2 * X^2 c B(4,3)(X) = 4 * (1-X) * X^3 c B(4,4)(X) = X^4 c c Special values: c c B(N,I)(X) has a unique maximum value at X = I/N. c c B(N,I)(X) has an I-fold zero at 0 and and N-I fold zero at 1. c c B(N,I)(1/2) = C(N,K) / 2^N c c For a fixed X and N, the polynomials add up to 1: c c Sum ( 0 <= I <= N ) B(N,I)(X) = 1 c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 20 March 2007 c c Author: c c John Burkardt c c Reference: c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer N, the degree of the polynomial. c c Output, integer K, the index of the polynomial. c c Output, double precision X, the argument of the polynomial. c c Output, double precision B, the value of the polynomial B(N,K)(X). c implicit none integer n_max parameter ( n_max = 15 ) double precision b double precision b_vec(n_max) integer k integer k_vec(n_max) integer n integer n_data integer n_vec(n_max) double precision x double precision x_vec(n_max) save b_vec save k_vec save n_vec save x_vec data b_vec / & 0.1000000000000000D+01, & 0.7500000000000000D+00, & 0.2500000000000000D+00, & 0.5625000000000000D+00, & 0.3750000000000000D+00, & 0.6250000000000000D-01, & 0.4218750000000000D+00, & 0.4218750000000000D+00, & 0.1406250000000000D+00, & 0.1562500000000000D-01, & 0.3164062500000000D+00, & 0.4218750000000000D+00, & 0.2109375000000000D+00, & 0.4687500000000000D-01, & 0.3906250000000000D-02 / data k_vec / & 0, & 0, 1, & 0, 1, 2, & 0, 1, 2, 3, & 0, 1, 2, 3, 4 / data n_vec / & 0, & 1, 1, & 2, 2, 2, & 3, 3, 3, 3, & 4, 4, 4, 4, 4 / data x_vec / & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 n = 0 k = 0 x = 0.0D+00 b = 0.0D+00 else n = n_vec(n_data) k = k_vec(n_data) x = x_vec(n_data) b = b_vec(n_data) end if return end subroutine bessel_i0_int_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_i0_int_values() returns some values of the Bessel I0 integral. c c Discussion: c c The function is defined by: c c I0_INT(x) = Integral ( 0 <= t <= x ) I0(t) dt c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.19531256208818052282D-02, & -0.39062549670565734544D-02, & 0.62520348032546565850D-01, & 0.12516285581366971819D+00, & -0.51051480879740303760D+00, & 0.10865210970235898158D+01, & 0.27750019054282535299D+01, & -0.13775208868039716639D+02, & 0.46424372058106108576D+03, & 0.64111867658021584522D+07, & -0.10414860803175857953D+08, & 0.44758598913855743089D+08, & -0.11852985311558287888D+09, & 0.31430078220715992752D+09, & -0.83440212900794309620D+09, & 0.22175367579074298261D+10, & 0.58991731842803636487D+10, & -0.41857073244691522147D+11, & 0.79553885818472357663D+12, & 0.15089715082719201025D+17 / data x_vec / & 0.0019531250D+00, & -0.0039062500D+00, & 0.0625000000D+00, & 0.1250000000D+00, & -0.5000000000D+00, & 1.0000000000D+00, & 2.0000000000D+00, & -4.0000000000D+00, & 8.0000000000D+00, & 18.0000000000D+00, & -18.5000000000D+00, & 20.0000000000D+00, & -21.0000000000D+00, & 22.0000000000D+00, & -23.0000000000D+00, & 24.0000000000D+00, & 25.0000000000D+00, & -27.0000000000D+00, & 30.0000000000D+00, & 40.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_i0_spherical_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_i0_spherical_values() returns some values of the Spherical Bessel function i0. c c Discussion: c c In Mathematica, the function can be evaluated by: c c Sqrt[Pi/(2*x)] * BesselI[1/2,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 06 January 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 21 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 1.001667500198440D+00, & 1.006680012705470D+00, & 1.026880814507039D+00, & 1.061089303580402D+00, & 1.110132477734529D+00, & 1.175201193643801D+00, & 1.257884462843477D+00, & 1.360215358179667D+00, & 1.484729970750144D+00, & 1.634541271164267D+00, & 1.813430203923509D+00, & 2.025956895698133D+00, & 2.277595505698373D+00, & 2.574897010920645D+00, & 2.925685126512827D+00, & 3.339291642469967D+00, & 3.826838748926716D+00, & 4.401577467270101D+00, & 5.079293155726485D+00, & 5.878791279137455D+00, & 6.822479299281938D+00 / data x_vec / & 0.1D+00, & 0.2D+00, & 0.4D+00, & 0.6D+00, & 0.8D+00, & 1.0D+00, & 1.2D+00, & 1.4D+00, & 1.6D+00, & 1.8D+00, & 2.0D+00, & 2.2D+00, & 2.4D+00, & 2.6D+00, & 2.8D+00, & 3.0D+00, & 3.2D+00, & 3.4D+00, & 3.6D+00, & 3.8D+00, & 4.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_i0_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_i0_values() returns some values of the I0 Bessel function. c c Discussion: c c The modified Bessel functions In(Z) and Kn(Z) are solutions of c the differential equation c c Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. c c The modified Bessel function I0(Z) corresponds to N = 0. c c In Mathematica, the function can be evaluated by: c c BesselI[0,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.1000000000000000D+01, & 0.1010025027795146D+01, & 0.1040401782229341D+01, & 0.1092045364317340D+01, & 0.1166514922869803D+01, & 0.1266065877752008D+01, & 0.1393725584134064D+01, & 0.1553395099731217D+01, & 0.1749980639738909D+01, & 0.1989559356618051D+01, & 0.2279585302336067D+01, & 0.3289839144050123D+01, & 0.4880792585865024D+01, & 0.7378203432225480D+01, & 0.1130192195213633D+02, & 0.1748117185560928D+02, & 0.2723987182360445D+02, & 0.6723440697647798D+02, & 0.4275641157218048D+03, & 0.2815716628466254D+04 / data x_vec / & 0.00D+00, & 0.20D+00, & 0.40D+00, & 0.60D+00, & 0.80D+00, & 0.10D+01, & 0.12D+01, & 0.14D+01, & 0.16D+01, & 0.18D+01, & 0.20D+01, & 0.25D+01, & 0.30D+01, & 0.35D+01, & 0.40D+01, & 0.45D+01, & 0.50D+01, & 0.60D+01, & 0.80D+01, & 0.10D+02 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_i1_spherical_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_i1_spherical_values() returns some values of the Spherical Bessel function i1. c c Discussion: c c In Mathematica, the function can be evaluated by: c c Sqrt[Pi/(2*x)] * BesselJ[3/2,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 06 January 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 21 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.03336667857363341D+00, & 0.06693371456802954D+00, & 0.1354788933285401D+00, & 0.2072931911031093D+00, & 0.2841280857128948D+00, & 0.3678794411714423D+00, & 0.4606425870674146D+00, & 0.5647736480096238D+00, & 0.6829590627779635D+00, & 0.8182955028627777D+00, & 0.9743827435800610D+00, & 1.155432469636406D+00, & 1.366396525527973D+00, & 1.613118767572064D+00, & 1.902515460838681D+00, & 2.242790117769266D+00, & 2.643689828630357D+00, & 3.116811526884873D+00, & 3.675968313148932D+00, & 4.337627987747642D+00, & 5.121438384183637D+00 / data x_vec / & 0.1D+00, & 0.2D+00, & 0.4D+00, & 0.6D+00, & 0.8D+00, & 1.0D+00, & 1.2D+00, & 1.4D+00, & 1.6D+00, & 1.8D+00, & 2.0D+00, & 2.2D+00, & 2.4D+00, & 2.6D+00, & 2.8D+00, & 3.0D+00, & 3.2D+00, & 3.4D+00, & 3.6D+00, & 3.8D+00, & 4.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_i1_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_i1_values() returns some values of the I1 Bessel function. c c Discussion: c c The modified Bessel functions In(Z) and Kn(Z) are solutions of c the differential equation c c Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. c c In Mathematica, the function can be evaluated by: c c BesselI[1,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.0000000000000000D+00, & 0.1005008340281251D+00, & 0.2040267557335706D+00, & 0.3137040256049221D+00, & 0.4328648026206398D+00, & 0.5651591039924850D+00, & 0.7146779415526431D+00, & 0.8860919814143274D+00, & 0.1084810635129880D+01, & 0.1317167230391899D+01, & 0.1590636854637329D+01, & 0.2516716245288698D+01, & 0.3953370217402609D+01, & 0.6205834922258365D+01, & 0.9759465153704450D+01, & 0.1538922275373592D+02, & 0.2433564214245053D+02, & 0.6134193677764024D+02, & 0.3998731367825601D+03, & 0.2670988303701255D+04 / data x_vec / & 0.00D+00, & 0.20D+00, & 0.40D+00, & 0.60D+00, & 0.80D+00, & 0.10D+01, & 0.12D+01, & 0.14D+01, & 0.16D+01, & 0.18D+01, & 0.20D+01, & 0.25D+01, & 0.30D+01, & 0.35D+01, & 0.40D+01, & 0.45D+01, & 0.50D+01, & 0.60D+01, & 0.80D+01, & 0.10D+02 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_in_values ( n_data, nu, x, fx ) c*********************************************************************72 c cc bessel_in_values() returns some values of the In Bessel function. c c Discussion: c c The modified Bessel functions In(Z) and Kn(Z) are solutions of c the differential equation c c Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. c c In Mathematica, the function can be evaluated by: c c BesselI[n,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer NU, the order of the function. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 28 ) double precision fx double precision fx_vec(n_max) integer n_data integer nu integer nu_vec(n_max) double precision x double precision x_vec(n_max) save fx_vec save nu_vec save x_vec data fx_vec / & 0.5016687513894678D-02, & 0.1357476697670383D+00, & 0.6889484476987382D+00, & 0.1276466147819164D+01, & 0.2245212440929951D+01, & 0.1750561496662424D+02, & 0.2281518967726004D+04, & 0.3931278522104076D+08, & 0.2216842492433190D-01, & 0.2127399592398527D+00, & 0.1033115016915114D+02, & 0.1758380716610853D+04, & 0.2677764138883941D+21, & 0.2714631559569719D-03, & 0.9825679323131702D-02, & 0.2157974547322546D+01, & 0.7771882864032600D+03, & 0.2278548307911282D+21, & 0.2752948039836874D-09, & 0.3016963879350684D-06, & 0.4580044419176051D-02, & 0.2189170616372337D+02, & 0.1071597159477637D+21, & 0.3966835985819020D-24, & 0.4310560576109548D-18, & 0.5024239357971806D-10, & 0.1250799735644948D-03, & 0.5442008402752998D+19 / data nu_vec / & 2, 2, 2, 2, & 2, 2, 2, 2, & 3, 3, 3, 3, & 3, 5, 5, 5, & 5, 5, 10, 10, & 10, 10, 10, 20, & 20, 20, 20, 20 / data x_vec / & 0.2D+00, & 1.0D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 5.0D+00, & 10.0D+00, & 20.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 nu = 0 x = 0.0D+00 fx = 0.0D+00 else nu = nu_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_ix_values ( n_data, nu, x, fx ) c*********************************************************************72 c cc bessel_ix_values() returns some values of the Ix Bessel function. c c Discussion: c c This set of data considers the less common case in which the c index of the Bessel function In is actually not an integer. c We may suggest this case by occasionally replacing the symbol c "In" by "Ix". c c The modified Bessel functions In(Z) and Kn(Z) are solutions of c the differential equation c c Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. c c In Mathematica, the function can be evaluated by: c c BesselI[n,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 02 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision NU, the order of the function. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 28 ) double precision fx double precision fx_vec(n_max) integer n_data double precision nu double precision nu_vec(n_max) double precision x double precision x_vec(n_max) save fx_vec save nu_vec save x_vec data fx_vec / & 0.3592084175833614D+00, & 0.9376748882454876D+00, & 2.046236863089055D+00, & 3.053093538196718D+00, & 4.614822903407601D+00, & 26.47754749755907D+00, & 2778.784603874571D+00, & 4.327974627242893D+07, & 0.2935253263474798D+00, & 1.099473188633110D+00, & 21.18444226479414D+00, & 2500.906154942118D+00, & 2.866653715931464D+20, & 0.05709890920304825D+00, & 0.3970270801393905D+00, & 13.76688213868258D+00, & 2028.512757391936D+00, & 2.753157630035402D+20, & 0.4139416015642352D+00, & 1.340196758982897D+00, & 22.85715510364670D+00, & 2593.006763432002D+00, & 2.886630075077766D+20, & 0.03590910483251082D+00, & 0.2931108636266483D+00, & 11.99397010023068D+00, & 1894.575731562383D+00, & 2.716911375760483D+20 / data nu_vec / & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00 / data x_vec / & 0.2D+00, & 1.0D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 5.0D+00, & 10.0D+00, & 20.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 nu = 0.0D+00 x = 0.0D+00 fx = 0.0D+00 else nu = nu_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_j0_int_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_j0_int_values() returns some values of the Bessel J0 integral. c c Discussion: c c The function is defined by: c c J0_INT(x) = Integral ( 0 <= t <= x ) J0(t) dt c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.97656242238978822427D-03, & 0.39062450329491108875D-02, & -0.62479657927917933620D-01, & 0.12483733492120479139D+00, & -0.48968050664604505505D+00, & 0.91973041008976023931D+00, & -0.14257702931970265690D+01, & 0.10247341594606064818D+01, & -0.12107468348304501655D+01, & 0.11008652032736190799D+01, & -0.10060334829904124192D+01, & 0.81330572662485953519D+00, & -0.10583788214211277585D+01, & 0.87101492116545875169D+00, & -0.88424908882547488420D+00, & 0.11257761503599914603D+01, & -0.90141212258183461184D+00, & 0.91441344369647797803D+00, & -0.94482281938334394886D+00, & 0.92266255696016607257D+00 / data x_vec / & 0.0009765625D+00, & 0.0039062500D+00, & -0.0625000000D+00, & 0.1250000000D+00, & -0.5000000000D+00, & 1.0000000000D+00, & -2.0000000000D+00, & 4.0000000000D+00, & -8.0000000000D+00, & 16.0000000000D+00, & -16.5000000000D+00, & 18.0000000000D+00, & -20.0000000000D+00, & 25.0000000000D+00, & -30.0000000000D+00, & 40.0000000000D+00, & -50.0000000000D+00, & 75.0000000000D+00, & -80.0000000000D+00, & 100.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_j0_spherical_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_j0_spherical_values() returns some values of the Spherical Bessel function j0. c c Discussion: c c In Mathematica, the function can be evaluated by: c c Sqrt[Pi/(2*x)] * BesselJ[1/2,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 21 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.9983341664682815D+00, & 0.9933466539753061D+00, & 0.9735458557716262D+00, & 0.9410707889917256D+00, & 0.8966951136244035D+00, & 0.8414709848078965D+00, & 0.7766992383060220D+00, & 0.7038926642774716D+00, & 0.6247335019009407D+00, & 0.5410264615989973D+00, & 0.4546487134128408D+00, & 0.3674983653725410D+00, & 0.2814429918963129D+00, & 0.1982697583928709D+00, & 0.1196386250556803D+00, & 0.4704000268662241D-01, & -0.1824191982111872D-01, & -0.7515914765495039D-01, & -0.1229223453596812D+00, & -0.1610152344586103D+00, & -0.1892006238269821D+00 / data x_vec / & 0.1D+00, & 0.2D+00, & 0.4D+00, & 0.6D+00, & 0.8D+00, & 1.0D+00, & 1.2D+00, & 1.4D+00, & 1.6D+00, & 1.8D+00, & 2.0D+00, & 2.2D+00, & 2.4D+00, & 2.6D+00, & 2.8D+00, & 3.0D+00, & 3.2D+00, & 3.4D+00, & 3.6D+00, & 3.8D+00, & 4.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_j0_zero_values ( n_data, k, x ) c*****************************************************************************80 c cc bessel_j0_zero_values() returns some values of J0 Bessel zeroes. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 18 March 2024 c c Author: c c John Burkardt c c Input: c c integer n_data: The user sets N_DATA to 0 before the first call. c c Output: c c integer n_data: On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c integer k: the index of the zero. c c double precision x: the location of the zero.. c implicit none integer n_max parameter ( n_max = 30 ) integer k integer k_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save x_vec save k_vec data x_vec / & 2.40482555769577276862163187933D+00, & 5.52007811028631064959660411281D+00, & 8.65372791291101221695419871266D+00, & 11.7915344390142816137430449119D+00, & 14.9309177084877859477625939974D+00, & 18.0710639679109225431478829756D+00, & 21.2116366298792589590783933505D+00, & 24.3524715307493027370579447632D+00, & 27.4934791320402547958772882346D+00, & 30.6346064684319751175495789269D+00, & 33.7758202135735686842385463467D+00, & 36.9170983536640439797694930633D+00, & 40.0584257646282392947993073740D+00, & 43.1997917131767303575240727287D+00, & 46.3411883716618140186857888791D+00, & 49.4826098973978171736027615332D+00, & 52.6240518411149960292512853804D+00, & 55.7655107550199793116834927735D+00, & 58.9069839260809421328344066346D+00, & 62.0484691902271698828525002646D+00, & 65.1899648002068604406360337425D+00, & 68.3314693298567982709923038400D+00, & 71.4729816035937328250630738561D+00, & 74.6145006437018378838205404693D+00, & 77.7560256303880550377393718912D+00, & 80.8975558711376278637721434909D+00, & 84.0390907769381901578796383480D+00, & 87.1806298436411536512618050690D+00, & 90.3221726372104800557177667775D+00, & 93.4637187819447741711905915439D+00 / data k_vec / & 1, & 2, & 3, & 4, & 5, & 6, & 7, & 8, & 9, & 10, & 11, & 12, & 13, & 14, & 15, & 16, & 17, & 18, & 19, & 20, & 21, & 22, & 23, & 24, & 25, & 26, & 27, & 28, & 29, & 30 / if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 k = 0 x = 0.0D+00 else k = k_vec(n_data) x = x_vec(n_data) end if return end subroutine bessel_j0_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_j0_values() returns some values of the J0 Bessel function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c BesselJ[0,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 21 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & -0.1775967713143383D+00, & -0.3971498098638474D+00, & -0.2600519549019334D+00, & 0.2238907791412357D+00, & 0.7651976865579666D+00, & 0.1000000000000000D+01, & 0.7651976865579666D+00, & 0.2238907791412357D+00, & -0.2600519549019334D+00, & -0.3971498098638474D+00, & -0.1775967713143383D+00, & 0.1506452572509969D+00, & 0.3000792705195556D+00, & 0.1716508071375539D+00, & -0.9033361118287613D-01, & -0.2459357644513483D+00, & -0.1711903004071961D+00, & 0.4768931079683354D-01, & 0.2069261023770678D+00, & 0.1710734761104587D+00, & -0.1422447282678077D-01 / data x_vec / & -5.0D+00, & -4.0D+00, & -3.0D+00, & -2.0D+00, & -1.0D+00, & 0.0D+00, & 1.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 6.0D+00, & 7.0D+00, & 8.0D+00, & 9.0D+00, & 10.0D+00, & 11.0D+00, & 12.0D+00, & 13.0D+00, & 14.0D+00, & 15.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_j1_spherical_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_j1_spherical_values() returns some values of the Spherical Bessel function j1. c c Discussion: c c In Mathematica, the function can be evaluated by: c c Sqrt[Pi/(2*x)] * BesselJ[3/2,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 21 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.3330001190255757D-01, & 0.6640038067032223D-01, & 0.1312121544218529D+00, & 0.1928919568034122D+00, & 0.2499855053465475D+00, & 0.3011686789397568D+00, & 0.3452845698577903D+00, & 0.3813753724123076D+00, & 0.4087081401263934D+00, & 0.4267936423844913D+00, & 0.4353977749799916D+00, & 0.4345452193763121D+00, & 0.4245152947656493D+00, & 0.4058301968314685D+00, & 0.3792360591872637D+00, & 0.3456774997623560D+00, & 0.3062665174917607D+00, & 0.2622467779189737D+00, & 0.2149544641595738D+00, & 0.1657769677515280D+00, & 0.1161107492591575D+00 / data x_vec / & 0.1D+00, & 0.2D+00, & 0.4D+00, & 0.6D+00, & 0.8D+00, & 1.0D+00, & 1.2D+00, & 1.4D+00, & 1.6D+00, & 1.8D+00, & 2.0D+00, & 2.2D+00, & 2.4D+00, & 2.6D+00, & 2.8D+00, & 3.0D+00, & 3.2D+00, & 3.4D+00, & 3.6D+00, & 3.8D+00, & 4.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_j1_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_j1_values() returns some values of the J1 Bessel function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c BesselJ[1,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 21 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.3275791375914652D+00, & 0.6604332802354914D-01, & -0.3390589585259365D+00, & -0.5767248077568734D+00, & -0.4400505857449335D+00, & 0.0000000000000000D+00, & 0.4400505857449335D+00, & 0.5767248077568734D+00, & 0.3390589585259365D+00, & -0.6604332802354914D-01, & -0.3275791375914652D+00, & -0.2766838581275656D+00, & -0.4682823482345833D-02, & 0.2346363468539146D+00, & 0.2453117865733253D+00, & 0.4347274616886144D-01, & -0.1767852989567215D+00, & -0.2234471044906276D+00, & -0.7031805212177837D-01, & 0.1333751546987933D+00, & 0.2051040386135228D+00 / data x_vec / & -5.0D+00, & -4.0D+00, & -3.0D+00, & -2.0D+00, & -1.0D+00, & 0.0D+00, & 1.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 6.0D+00, & 7.0D+00, & 8.0D+00, & 9.0D+00, & 10.0D+00, & 11.0D+00, & 12.0D+00, & 13.0D+00, & 14.0D+00, & 15.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_jn_values ( n_data, nu, x, fx ) c*********************************************************************72 c cc bessel_jn_values() returns some values of the Jn Bessel function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c BesselJ[n,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer NU, the order of the function. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data integer nu integer nu_vec(n_max) double precision x double precision x_vec(n_max) save fx_vec save nu_vec save x_vec data fx_vec / & 0.1149034849319005D+00, & 0.3528340286156377D+00, & 0.4656511627775222D-01, & 0.2546303136851206D+00, & -0.5971280079425882D-01, & 0.2497577302112344D-03, & 0.7039629755871685D-02, & 0.2611405461201701D+00, & -0.2340615281867936D+00, & -0.8140024769656964D-01, & 0.2630615123687453D-09, & 0.2515386282716737D-06, & 0.1467802647310474D-02, & 0.2074861066333589D+00, & -0.1138478491494694D+00, & 0.3873503008524658D-24, & 0.3918972805090754D-18, & 0.2770330052128942D-10, & 0.1151336924781340D-04, & -0.1167043527595797D+00 / data nu_vec / & 2, 2, 2, 2, & 2, 5, 5, 5, & 5, 5, 10, 10, & 10, 10, 10, 20, & 20, 20, 20, 20 / data x_vec / & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 nu = 0 x = 0.0D+00 fx = 0.0D+00 else nu = nu_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_jx_values ( n_data, nu, x, fx ) c*********************************************************************72 c cc bessel_jx_values() returns some values of the Jx Bessel function. c c Discussion: c c This set of data considers the less common case in which the c index of the Bessel function Jn is actually not an integer. c We may suggest this case by occasionally replacing the symbol c "Jn" by "Jx". c c In Mathematica, the function can be evaluated by: c c BesselJ[n,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 31 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision NU, the order of the function. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 28 ) double precision fx double precision fx_vec(n_max) integer n_data double precision nu double precision nu_vec(n_max) double precision x double precision x_vec(n_max) save fx_vec save nu_vec save x_vec data fx_vec / & 0.3544507442114011D+00, & 0.6713967071418031D+00, & 0.5130161365618278D+00, & 0.3020049060623657D+00, & 0.06500818287737578D+00, & -0.3421679847981618D+00, & -0.1372637357550505D+00, & 0.1628807638550299D+00, & 0.2402978391234270D+00, & 0.4912937786871623D+00, & -0.1696513061447408D+00, & 0.1979824927558931D+00, & -0.1094768729883180D+00, & 0.04949681022847794D+00, & 0.2239245314689158D+00, & 0.2403772011113174D+00, & 0.1966584835818184D+00, & 0.02303721950962553D+00, & 0.3314145508558904D+00, & 0.5461734240402840D+00, & -0.2616584152094124D+00, & 0.1296035513791289D+00, & -0.1117432171933552D+00, & 0.03142623570527935D+00, & 0.1717922192746527D+00, & 0.3126634069544786D+00, & 0.1340289119304364D+00, & 0.06235967135106445D+00 / data nu_vec / & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00 / data x_vec / & 0.2D+00, & 1.0D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 5.0D+00, & 10.0D+00, & 20.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 nu = 0.0D+00 x = 0.0D+00 fx = 0.0D+00 else nu = nu_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_k0_int_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_k0_int_values() returns some values of the Bessel K0 integral. c c Discussion: c c The function is defined by: c c K0_INT(x) = Integral ( 0 <= t <= x ) K0(t) dt c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.78587929563466784589D-02, & 0.26019991617330578111D-01, & 0.24311842237541167904D+00, & 0.39999633750480508861D+00, & 0.92710252093114907345D+00, & 0.12425098486237782662D+01, & 0.14736757343168286825D+01, & 0.15606495706051741364D+01, & 0.15673873907283660493D+01, & 0.15696345532693743714D+01, & 0.15701153443250786355D+01, & 0.15706574852894436220D+01, & 0.15707793116159788598D+01, & 0.15707942066465767196D+01, & 0.15707962315469192247D+01, & 0.15707963262340149876D+01, & 0.15707963267948756308D+01, & 0.15707963267948966192D+01, & 0.15707963267948966192D+01, & 0.15707963267948966192D+01 / data x_vec / & 0.0009765625D+00, & 0.0039062500D+00, & 0.0625000000D+00, & 0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 2.0000000000D+00, & 4.0000000000D+00, & 5.0000000000D+00, & 6.0000000000D+00, & 6.5000000000D+00, & 8.0000000000D+00, & 10.0000000000D+00, & 12.0000000000D+00, & 15.0000000000D+00, & 20.0000000000D+00, & 30.0000000000D+00, & 50.0000000000D+00, & 80.0000000000D+00, & 100.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_k0_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_k0_values() returns some values of the K0 Bessel function. c c Discussion: c c The modified Bessel functions In(Z) and Kn(Z) are solutions of c the differential equation c c Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. c c The modified Bessel function K0(Z) corresponds to N = 0. c c In Mathematica, the function can be evaluated by: c c BesselK[0,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.2427069024702017D+01, & 0.1752703855528146D+01, & 0.1114529134524434D+01, & 0.7775220919047293D+00, & 0.5653471052658957D+00, & 0.4210244382407083D+00, & 0.3185082202865936D+00, & 0.2436550611815419D+00, & 0.1879547519693323D+00, & 0.1459314004898280D+00, & 0.1138938727495334D+00, & 0.6234755320036619D-01, & 0.3473950438627925D-01, & 0.1959889717036849D-01, & 0.1115967608585302D-01, & 0.6399857243233975D-02, & 0.3691098334042594D-02, & 0.1243994328013123D-02, & 0.1464707052228154D-03, & 0.1778006231616765D-04 / data x_vec / & 0.1D+00, & 0.2D+00, & 0.4D+00, & 0.6D+00, & 0.8D+00, & 1.0D+00, & 1.2D+00, & 1.4D+00, & 1.6D+00, & 1.8D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 3.5D+00, & 4.0D+00, & 4.5D+00, & 5.0D+00, & 6.0D+00, & 8.0D+00, & 10.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_k1_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_k1_values() returns some values of the K1 Bessel function. c c Discussion: c c The modified Bessel functions In(Z) and Kn(Z) are solutions of c the differential equation c c Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. c c The modified Bessel function K1(Z) corresponds to N = 1. c c In Mathematica, the function can be evaluated by: c c BesselK[1,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & 0.9853844780870606D+01, & 0.4775972543220472D+01, & 0.2184354424732687D+01, & 0.1302834939763502D+01, & 0.8617816344721803D+00, & 0.6019072301972346D+00, & 0.4345923910607150D+00, & 0.3208359022298758D+00, & 0.2406339113576119D+00, & 0.1826230998017470D+00, & 0.1398658818165224D+00, & 0.7389081634774706D-01, & 0.4015643112819418D-01, & 0.2223939292592383D-01, & 0.1248349888726843D-01, & 0.7078094908968090D-02, & 0.4044613445452164D-02, & 0.1343919717735509D-02, & 0.1553692118050011D-03, & 0.1864877345382558D-04 / data x_vec / & 0.1D+00, & 0.2D+00, & 0.4D+00, & 0.6D+00, & 0.8D+00, & 1.0D+00, & 1.2D+00, & 1.4D+00, & 1.6D+00, & 1.8D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 3.5D+00, & 4.0D+00, & 4.5D+00, & 5.0D+00, & 6.0D+00, & 8.0D+00, & 10.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_kn_values ( n_data, nu, x, fx ) c*********************************************************************72 c cc bessel_kn_values() returns some values of the Kn Bessel function. c c Discussion: c c The modified Bessel functions In(Z) and Kn(Z) are solutions of c the differential equation c c Z^2 * W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. c c In Mathematica, the function can be evaluated by: c c BesselK[n,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer NU, the order of the function. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 28 ) double precision fx double precision fx_vec(n_max) integer n_data integer nu integer nu_vec(n_max) double precision x double precision x_vec(n_max) save fx_vec save nu_vec save x_vec data fx_vec / & 0.4951242928773287D+02, & 0.1624838898635177D+01, & 0.2537597545660559D+00, & 0.1214602062785638D+00, & 0.6151045847174204D-01, & 0.5308943712223460D-02, & 0.2150981700693277D-04, & 0.6329543612292228D-09, & 0.7101262824737945D+01, & 0.6473853909486342D+00, & 0.8291768415230932D-02, & 0.2725270025659869D-04, & 0.3727936773826211D-22, & 0.3609605896012407D+03, & 0.9431049100596467D+01, & 0.3270627371203186D-01, & 0.5754184998531228D-04, & 0.4367182254100986D-22, & 0.1807132899010295D+09, & 0.1624824039795591D+06, & 0.9758562829177810D+01, & 0.1614255300390670D-02, & 0.9150988209987996D-22, & 0.6294369360424535D+23, & 0.5770856852700241D+17, & 0.4827000520621485D+09, & 0.1787442782077055D+03, & 0.1706148379722035D-20 / data nu_vec / & 2, 2, 2, 2, & 2, 2, 2, 2, & 3, 3, 3, 3, & 3, 5, 5, 5, & 5, 5, 10, 10, & 10, 10, 10, 20, & 20, 20, 20, 20 / data x_vec / & 0.2D+00, & 1.0D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 5.0D+00, & 10.0D+00, & 20.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 nu = 0 x = 0.0D+00 fx = 0.0D+00 else nu = nu_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_kx_values ( n_data, nu, x, fx ) c*********************************************************************72 c cc bessel_kx_values() returns some values of the Kx Bessel function. c c Discussion: c c This set of data considers the less common case in which the c index of the Bessel function Kn is actually not an integer. c We may suggest this case by occasionally replacing the symbol c "Kn" by "Kx". c c The modified Bessel functions In(Z) and Kn(Z) are solutions of c the differential equation c c Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. c c In Mathematica, the function can be evaluated by: c c BesselK[n,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 31 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision NU, the order of the function. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 28 ) double precision fx double precision fx_vec(n_max) integer n_data double precision nu double precision nu_vec(n_max) double precision x double precision x_vec(n_max) save fx_vec save nu_vec save x_vec data fx_vec / & 2.294489339798475D+00, & 0.4610685044478946D+00, & 0.1199377719680614D+00, & 0.06506594315400999D+00, & 0.03602598513176459D+00, & 0.003776613374642883D+00, & 0.00001799347809370518D+00, & 5.776373974707445D-10, & 0.9221370088957891D+00, & 0.1799066579520922D+00, & 0.004531936049571459D+00, & 0.00001979282590307570D+00, & 3.486992497366216D-23, & 3.227479531135262D+00, & 0.3897977588961997D+00, & 0.006495775004385758D+00, & 0.00002393132586462789D+00, & 3.627839645299048D-23, & 0.7311451879202114D+00, & 0.1567475478393932D+00, & 0.004257389528177461D+00, & 0.00001915541065869563D+00, & 3.463337593569306D-23, & 4.731184839919541D+00, & 0.4976876225514758D+00, & 0.007300864610941163D+00, & 0.00002546421294106458D+00, & 3.675275677913656D-23 / data nu_vec / & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00 / data x_vec / & 0.2D+00, & 1.0D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 5.0D+00, & 10.0D+00, & 20.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 nu = 0.0D+00 x = 0.0D+00 fx = 0.0D+00 else nu = nu_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_y0_int_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_y0_int_values() returns some values of the Bessel Y0 integral. c c Discussion: c c The function is defined by: c c Y0_INT(x) = Integral ( 0 <= t <= x ) Y0(t) dt c c The data was reported by McLeod. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Allan McLeod, c Algorithm 757: c MISCFUN: A software package to compute uncommon special functions, c ACM Transactions on Mathematical Software, c Volume 22, Number 3, September 1996, pages 288-301. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & -0.91442642860172110926D-02, & -0.29682047390397591290D-01, & -0.25391431276585388961D+00, & -0.56179545591464028187D+00, & -0.63706937660742309754D+00, & -0.28219285008510084123D+00, & 0.38366964785312561103D+00, & -0.12595061285798929390D+00, & 0.24129031832266684828D+00, & 0.17138069757627037938D+00, & 0.18958142627134083732D+00, & 0.17203846136449706946D+00, & -0.16821597677215029611D+00, & -0.93607927351428988679D-01, & 0.88229711948036648408D-01, & -0.89324662736274161841D-02, & -0.54814071000063488284D-01, & -0.94958246003466381588D-01, & -0.19598064853404969850D-01, & -0.83084772357154773468D-02 / data x_vec / & 0.0019531250D+00, & 0.0078125000D+00, & 0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 2.0000000000D+00, & 4.0000000000D+00, & 6.0000000000D+00, & 10.0000000000D+00, & 16.0000000000D+00, & 16.2500000000D+00, & 17.0000000000D+00, & 20.0000000000D+00, & 25.0000000000D+00, & 30.0000000000D+00, & 40.0000000000D+00, & 50.0000000000D+00, & 70.0000000000D+00, & 100.0000000000D+00, & 125.0000000000D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_y0_spherical_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_y0_spherical_values() returns some values of the Spherical Bessel function y0. c c Discussion: c c In Mathematica, the function can be evaluated by: c c Sqrt[Pi/(2*x)] * BesselY[1/2,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 21 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & -0.9950041652780258D+01, & -0.4900332889206208D+01, & -0.2302652485007213D+01, & -0.1375559358182797D+01, & -0.8708833866839568D+00, & -0.5403023058681397D+00, & -0.3019647953972280D+00, & -0.1214051020716007D+00, & 0.1824970143830545D-01, & 0.1262233859406039D+00, & 0.2080734182735712D+00, & 0.2675005078433390D+00, & 0.3072473814755190D+00, & 0.3295725974495951D+00, & 0.3365079788102351D+00, & 0.3299974988668152D+00, & 0.3119671174358603D+00, & 0.2843524095821944D+00, & 0.2490995600928186D+00, & 0.2081493978722149D+00, & 0.1634109052159030D+00 / data x_vec / & 0.1D+00, & 0.2D+00, & 0.4D+00, & 0.6D+00, & 0.8D+00, & 1.0D+00, & 1.2D+00, & 1.4D+00, & 1.6D+00, & 1.8D+00, & 2.0D+00, & 2.2D+00, & 2.4D+00, & 2.6D+00, & 2.8D+00, & 3.0D+00, & 3.2D+00, & 3.4D+00, & 3.6D+00, & 3.8D+00, & 4.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_y0_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_y0_values() returns some values of the Y0 Bessel function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c BesselY[0,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 16 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & -0.1534238651350367D+01, & 0.8825696421567696D-01, & 0.5103756726497451D+00, & 0.3768500100127904D+00, & -0.1694073932506499D-01, & -0.3085176252490338D+00, & -0.2881946839815792D+00, & -0.2594974396720926D-01, & 0.2235214893875662D+00, & 0.2499366982850247D+00, & 0.5567116728359939D-01, & -0.1688473238920795D+00, & -0.2252373126343614D+00, & -0.7820786452787591D-01, & 0.1271925685821837D+00, & 0.2054642960389183D+00 / data x_vec / & 0.1D+00, & 1.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 6.0D+00, & 7.0D+00, & 8.0D+00, & 9.0D+00, & 10.0D+00, & 11.0D+00, & 12.0D+00, & 13.0D+00, & 14.0D+00, & 15.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_y0_zero_values ( n_data, k, x ) c*****************************************************************************80 c cc bessel_y0_zero_values() returns some values of Y0 Bessel zeroes. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 19 March 2024 c c Author: c c John Burkardt c c Input: c c integer n_data: The user sets N_DATA to 0 before the first call. c c Output: c c integer n_data: On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c integer k: the index of the zero. c c double precision x: the location of the zero.. c implicit none integer n_max parameter ( n_max = 32 ) integer k integer k_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save x_vec save k_vec data x_vec / & 0.8935769662791675D+00, & 3.957678419314858D+00, & 7.086051060301773D+00, & 10.22234504349642D+00, & 13.36109747387276D+00, & 16.50092244152809D+00, & 19.64130970088794D+00, & 22.78202804729156D+00, & 25.92295765318092D+00, & 29.06403025272840D+00, & 32.20520411649328D+00, & 35.34645230521432D+00, & 38.48775665308154D+00, & 41.62910446621381D+00, & 44.77048660722199D+00, & 47.91189633151648D+00, & 51.05332855236236D+00, & 54.19477936108705D+00, & 57.33624570476628D+00, & 60.47772516422348D+00, & 63.61921579772038D+00, & 66.76071602872964D+00, & 69.90222456393850D+00, & 73.04374033239207D+00, & 76.18526243968061D+00, & 79.32679013300400D+00, & 82.46832277421550D+00, & 85.60985981879671D+00, & 88.75140079929514D+00, & 91.89294531215718D+00, & 95.03449300717121D+00, & 98.17604357893667D+00 / data k_vec / & 1, & 2, & 3, & 4, & 5, & 6, & 7, & 8, & 9, & 10, & 11, & 12, & 13, & 14, & 15, & 16, & 17, & 18, & 19, & 20, & 21, & 22, & 23, & 24, & 25, & 26, & 27, & 28, & 29, & 30, & 31, & 32 / if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 k = 0 x = 0.0D+00 else k = k_vec(n_data) x = x_vec(n_data) end if return end subroutine bessel_y1_spherical_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_y1_spherical_values() returns some values of the Spherical Bessel function y1. c c Discussion: c c In Mathematica, the function can be evaluated by: c c Sqrt[Pi/(2*x)] * BesselY[3/2,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 21 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & -0.1004987506942709D+03, & -0.2549501110000635D+02, & -0.6730177068289658D+01, & -0.3233669719296388D+01, & -0.1985299346979349D+01, & -0.1381773290676036D+01, & -0.1028336567803712D+01, & -0.7906105943286149D+00, & -0.6133274385019998D+00, & -0.4709023582986618D+00, & -0.3506120042760553D+00, & -0.2459072254437506D+00, & -0.1534232496148467D+00, & -0.7151106706610352D-01, & 0.5427959479750482D-03, & 0.6295916360231598D-01, & 0.1157316440198251D+00, & 0.1587922092967723D+00, & 0.1921166676076864D+00, & 0.2157913917934037D+00, & 0.2300533501309578D+00 / data x_vec / & 0.1D+00, & 0.2D+00, & 0.4D+00, & 0.6D+00, & 0.8D+00, & 1.0D+00, & 1.2D+00, & 1.4D+00, & 1.6D+00, & 1.8D+00, & 2.0D+00, & 2.2D+00, & 2.4D+00, & 2.6D+00, & 2.8D+00, & 3.0D+00, & 3.2D+00, & 3.4D+00, & 3.6D+00, & 3.8D+00, & 4.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_y1_values ( n_data, x, fx ) c*********************************************************************72 c cc bessel_y1_values() returns some values of the Y1 Bessel function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c BesselY[1,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 16 ) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save fx_vec save x_vec data fx_vec / & -0.6458951094702027D+01, & -0.7812128213002887D+00, & -0.1070324315409375D+00, & 0.3246744247918000D+00, & 0.3979257105571000D+00, & 0.1478631433912268D+00, & -0.1750103443003983D+00, & -0.3026672370241849D+00, & -0.1580604617312475D+00, & 0.1043145751967159D+00, & 0.2490154242069539D+00, & 0.1637055374149429D+00, & -0.5709921826089652D-01, & -0.2100814084206935D+00, & -0.1666448418561723D+00, & 0.2107362803687351D-01 / data x_vec / & 0.1D+00, & 1.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 6.0D+00, & 7.0D+00, & 8.0D+00, & 9.0D+00, & 10.0D+00, & 11.0D+00, & 12.0D+00, & 13.0D+00, & 14.0D+00, & 15.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_yn_values ( n_data, nu, x, fx ) c*********************************************************************72 c cc bessel_yn_values() returns some values of the Yn Bessel function. c c Discussion: c c In Mathematica, the function can be evaluated by: c c BesselY[n,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 January 2006 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer NU, the order of the function. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) double precision fx double precision fx_vec(n_max) integer n_data integer nu integer nu_vec(n_max) double precision x double precision x_vec(n_max) save fx_vec save nu_vec save x_vec data fx_vec / & -0.1650682606816254D+01, & -0.6174081041906827D+00, & 0.3676628826055245D+00, & -0.5868082442208615D-02, & 0.9579316872759649D-01, & -0.2604058666258122D+03, & -0.9935989128481975D+01, & -0.4536948224911019D+00, & 0.1354030476893623D+00, & -0.7854841391308165D-01, & -0.1216180142786892D+09, & -0.1291845422080393D+06, & -0.2512911009561010D+02, & -0.3598141521834027D+00, & 0.5723897182053514D-02, & -0.4113970314835505D+23, & -0.4081651388998367D+17, & -0.5933965296914321D+09, & -0.1597483848269626D+04, & 0.1644263394811578D-01 / data nu_vec / & 2, 2, 2, 2, & 2, 5, 5, 5, & 5, 5, 10, 10, & 10, 10, 10, 20, & 20, 20, 20, 20 / data x_vec / & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 nu = 0 x = 0.0D+00 fx = 0.0D+00 else nu = nu_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bessel_yx_values ( n_data, nu, x, fx ) c*********************************************************************72 c cc bessel_yx_values() returns some values of the Yx Bessel function. c c Discussion: c c This set of data considers the less common case in which the c index of the Bessel function Yn is actually not an integer. c We may suggest this case by occasionally replacing the symbol c "Yn" by "Yx". c c In Mathematica, the function can be evaluated by: c c BesselY[n,x] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 31 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision NU, the order of the function. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 28 ) double precision fx double precision fx_vec(n_max) integer n_data double precision nu double precision nu_vec(n_max) double precision x double precision x_vec(n_max) save fx_vec save nu_vec save x_vec data fx_vec / & -1.748560416961876D+00, & -0.4310988680183761D+00, & 0.2347857104062485D+00, & 0.4042783022390569D+00, & 0.4560488207946332D+00, & -0.1012177091851084D+00, & 0.2117088663313982D+00, & -0.07280690478506185D+00, & -1.102495575160179D+00, & -0.3956232813587035D+00, & 0.3219244429611401D+00, & 0.1584346223881903D+00, & 0.02742813676191382D+00, & -2.876387857462161D+00, & -0.8282206324443037D+00, & 0.2943723749617925D+00, & -0.1641784796149411D+00, & 0.1105304445562544D+00, & -0.9319659251969881D+00, & -0.2609445010948933D+00, & 0.2492796362185881D+00, & 0.2174410301416733D+00, & -0.01578576650557229D+00, & -4.023453301501028D+00, & -0.9588998694752389D+00, & 0.2264260361047367D+00, & -0.2193617736566760D+00, & 0.09413988344515077D+00 / data nu_vec / & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 0.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 1.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 2.50D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 1.25D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00, & 2.75D+00 / data x_vec / & 0.2D+00, & 1.0D+00, & 2.0D+00, & 2.5D+00, & 3.0D+00, & 5.0D+00, & 10.0D+00, & 20.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00, & 1.0D+00, & 2.0D+00, & 5.0D+00, & 10.0D+00, & 50.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 nu = 0.0D+00 x = 0.0D+00 fx = 0.0D+00 else nu = nu_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine beta_cdf_values ( n_data, a, b, x, fx ) c*********************************************************************72 c cc beta_cdf_values() returns some values of the Beta CDF. c c Discussion: c c The incomplete Beta function may be written c c BETA_INC(A,B,X) = Integral (0 to X) T^(A-1) * (1-T)^(B-1) dT c / Integral (0 to 1) T^(A-1) * (1-T)^(B-1) dT c c Thus, c c BETA_INC(A,B,0.0) = 0.0 c BETA_INC(A,B,1.0) = 1.0 c c The incomplete Beta function is also sometimes called the c "modified" Beta function, or the "normalized" Beta function c or the Beta CDF (cumulative density function. c c In Mathematica, the function can be evaluated by: c c BETA[X,A,B] / BETA[A,B] c c The function can also be evaluated by using the Statistics package: c c Needs["Statistics`ContinuousDistributions`"] c dist = BetaDistribution [ a, b ] c CDF [ dist, x ] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 28 April 2013 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Karl Pearson, c Tables of the Incomplete Beta Function, c Cambridge University Press, 1968. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision A, B, the parameters of the function. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 45 ) double precision a double precision a_vec(n_max) double precision b double precision b_vec(n_max) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save a_vec save b_vec save fx_vec save x_vec data a_vec / & 0.5D+00, & 0.5D+00, & 0.5D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 5.5D+00, & 10.0D+00, & 10.0D+00, & 10.0D+00, & 10.0D+00, & 20.0D+00, & 20.0D+00, & 20.0D+00, & 20.0D+00, & 20.0D+00, & 30.0D+00, & 30.0D+00, & 40.0D+00, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.2D+01, & 0.3D+01, & 0.4D+01, & 0.5D+01, & 1.30625D+00, & 1.30625D+00, & 1.30625D+00 / data b_vec / & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 1.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 5.0D+00, & 0.5D+00, & 5.0D+00, & 5.0D+00, & 10.0D+00, & 5.0D+00, & 10.0D+00, & 10.0D+00, & 20.0D+00, & 20.0D+00, & 10.0D+00, & 10.0D+00, & 20.0D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.2D+01, & 0.3D+01, & 0.4D+01, & 0.5D+01, & 0.2D+01, & 0.2D+01, & 0.2D+01, & 0.2D+01, & 11.7562D+00, & 11.7562D+00, & 11.7562D+00 / data fx_vec / & 0.6376856085851985D-01, & 0.2048327646991335D+00, & 0.1000000000000000D+01, & 0.0000000000000000D+00, & 0.5012562893380045D-02, & 0.5131670194948620D-01, & 0.2928932188134525D+00, & 0.5000000000000000D+00, & 0.2800000000000000D-01, & 0.1040000000000000D+00, & 0.2160000000000000D+00, & 0.3520000000000000D+00, & 0.5000000000000000D+00, & 0.6480000000000000D+00, & 0.7840000000000000D+00, & 0.8960000000000000D+00, & 0.9720000000000000D+00, & 0.4361908850559777D+00, & 0.1516409096347099D+00, & 0.8978271484375000D-01, & 0.1000000000000000D+01, & 0.5000000000000000D+00, & 0.4598773297575791D+00, & 0.2146816102371739D+00, & 0.9507364826957875D+00, & 0.5000000000000000D+00, & 0.8979413687105918D+00, & 0.2241297491808366D+00, & 0.7586405487192086D+00, & 0.7001783247477069D+00, & 0.5131670194948620D-01, & 0.1055728090000841D+00, & 0.1633399734659245D+00, & 0.2254033307585166D+00, & 0.3600000000000000D+00, & 0.4880000000000000D+00, & 0.5904000000000000D+00, & 0.6723200000000000D+00, & 0.2160000000000000D+00, & 0.8370000000000000D-01, & 0.3078000000000000D-01, & 0.1093500000000000D-01, & 0.918884684620518D+00, & 0.21052977489419D+00, & 0.1824130512500673D+00 / data x_vec / & 0.01D+00, & 0.10D+00, & 1.00D+00, & 0.00D+00, & 0.01D+00, & 0.10D+00, & 0.50D+00, & 0.50D+00, & 0.10D+00, & 0.20D+00, & 0.30D+00, & 0.40D+00, & 0.50D+00, & 0.60D+00, & 0.70D+00, & 0.80D+00, & 0.90D+00, & 0.50D+00, & 0.90D+00, & 0.50D+00, & 1.00D+00, & 0.50D+00, & 0.80D+00, & 0.60D+00, & 0.80D+00, & 0.50D+00, & 0.60D+00, & 0.70D+00, & 0.80D+00, & 0.70D+00, & 0.10D+00, & 0.20D+00, & 0.30D+00, & 0.40D+00, & 0.20D+00, & 0.20D+00, & 0.20D+00, & 0.20D+00, & 0.30D+00, & 0.30D+00, & 0.30D+00, & 0.30D+00, & 0.225609D+00, & 0.0335568D+00, & 0.0295222D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 a = 0.0D+00 b = 0.0D+00 x = 0.0D+00 fx = 0.0D+00 else a = a_vec(n_data) b = b_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine beta_inc_values ( n_data, a, b, x, fx ) c*********************************************************************72 c cc beta_inc_values() returns some values of the incomplete Beta function. c c Discussion: c c The incomplete Beta function may be written c c BETA_INC(A,B,X) = Integral (0 to X) T^(A-1) * (1-T)^(B-1) dT c / Integral (0 to 1) T^(A-1) * (1-T)^(B-1) dT c c Thus, c c BETA_INC(A,B,0.0) = 0.0 c BETA_INC(A,B,1.0) = 1.0 c c The incomplete Beta function is also sometimes called the c "modified" Beta function, or the "normalized" Beta function c or the Beta CDF (cumulative density function. c c In Mathematica, the function can be evaluated by: c c BETA[X,A,B] / BETA[A,B] c c The function can also be evaluated by using the Statistics package: c c Needs["Statistics`ContinuousDistributions`"] c dist = BetaDistribution [ a, b ] c CDF [ dist, x ] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 28 April 2013 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Karl Pearson, c Tables of the Incomplete Beta Function, c Cambridge University Press, 1968. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision A, B, the parameters of the function. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 45 ) double precision a double precision a_vec(n_max) double precision b double precision b_vec(n_max) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save a_vec save b_vec save fx_vec save x_vec data a_vec / & 0.5D+00, & 0.5D+00, & 0.5D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 5.5D+00, & 10.0D+00, & 10.0D+00, & 10.0D+00, & 10.0D+00, & 20.0D+00, & 20.0D+00, & 20.0D+00, & 20.0D+00, & 20.0D+00, & 30.0D+00, & 30.0D+00, & 40.0D+00, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.2D+01, & 0.3D+01, & 0.4D+01, & 0.5D+01, & 1.30625D+00, & 1.30625D+00, & 1.30625D+00 / data b_vec / & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 1.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 2.0D+00, & 5.0D+00, & 0.5D+00, & 5.0D+00, & 5.0D+00, & 10.0D+00, & 5.0D+00, & 10.0D+00, & 10.0D+00, & 20.0D+00, & 20.0D+00, & 10.0D+00, & 10.0D+00, & 20.0D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.5D+00, & 0.2D+01, & 0.3D+01, & 0.4D+01, & 0.5D+01, & 0.2D+01, & 0.2D+01, & 0.2D+01, & 0.2D+01, & 11.7562D+00, & 11.7562D+00, & 11.7562D+00 / data fx_vec / & 0.6376856085851985D-01, & 0.2048327646991335D+00, & 0.1000000000000000D+01, & 0.0000000000000000D+00, & 0.5012562893380045D-02, & 0.5131670194948620D-01, & 0.2928932188134525D+00, & 0.5000000000000000D+00, & 0.2800000000000000D-01, & 0.1040000000000000D+00, & 0.2160000000000000D+00, & 0.3520000000000000D+00, & 0.5000000000000000D+00, & 0.6480000000000000D+00, & 0.7840000000000000D+00, & 0.8960000000000000D+00, & 0.9720000000000000D+00, & 0.4361908850559777D+00, & 0.1516409096347099D+00, & 0.8978271484375000D-01, & 0.1000000000000000D+01, & 0.5000000000000000D+00, & 0.4598773297575791D+00, & 0.2146816102371739D+00, & 0.9507364826957875D+00, & 0.5000000000000000D+00, & 0.8979413687105918D+00, & 0.2241297491808366D+00, & 0.7586405487192086D+00, & 0.7001783247477069D+00, & 0.5131670194948620D-01, & 0.1055728090000841D+00, & 0.1633399734659245D+00, & 0.2254033307585166D+00, & 0.3600000000000000D+00, & 0.4880000000000000D+00, & 0.5904000000000000D+00, & 0.6723200000000000D+00, & 0.2160000000000000D+00, & 0.8370000000000000D-01, & 0.3078000000000000D-01, & 0.1093500000000000D-01, & 0.918884684620518D+00, & 0.21052977489419D+00, & 0.1824130512500673D+00 / data x_vec / & 0.01D+00, & 0.10D+00, & 1.00D+00, & 0.00D+00, & 0.01D+00, & 0.10D+00, & 0.50D+00, & 0.50D+00, & 0.10D+00, & 0.20D+00, & 0.30D+00, & 0.40D+00, & 0.50D+00, & 0.60D+00, & 0.70D+00, & 0.80D+00, & 0.90D+00, & 0.50D+00, & 0.90D+00, & 0.50D+00, & 1.00D+00, & 0.50D+00, & 0.80D+00, & 0.60D+00, & 0.80D+00, & 0.50D+00, & 0.60D+00, & 0.70D+00, & 0.80D+00, & 0.70D+00, & 0.10D+00, & 0.20D+00, & 0.30D+00, & 0.40D+00, & 0.20D+00, & 0.20D+00, & 0.20D+00, & 0.20D+00, & 0.30D+00, & 0.30D+00, & 0.30D+00, & 0.30D+00, & 0.225609D+00, & 0.0335568D+00, & 0.0295222D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 a = 0.0D+00 b = 0.0D+00 x = 0.0D+00 fx = 0.0D+00 else a = a_vec(n_data) b = b_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine beta_log_values ( n_data, x, y, fxy ) c*********************************************************************72 c cc beta_log_values() returns some values of the logarithm of the Beta function. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 27 March 2010 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Parameters: c c Input/output, integer N_DATA. c On input, if N_DATA is 0, the first test data is returned, and N_DATA is set c to the index of the test data. On each subsequent call, N_DATA is c incremented and that test data is returned. When there is no more c test data, N_DATA is set to 0. c c Output, double precision X, Y, the arguments of the function. c c Output, double precision FXY, the value of the function. c implicit none integer n_max parameter ( n_max = 17 ) double precision fxvec ( n_max ) double precision fxy integer n_data double precision x double precision xvec ( n_max ) double precision y double precision yvec ( n_max ) data fxvec / & 1.609437912434100D+00, & 0.9162907318741551D+00, & 0.5108256237659907D+00, & 0.2231435513142098D+00, & 1.609437912434100D+00, & 0.9162907318741551D+00, & 0.000000000000000D+00, & -1.791759469228055D+00, & -3.401197381662155D+00, & -4.941642422609304D+00, & -6.445719819385578D+00, & -3.737669618283368D+00, & -5.123963979403259D+00, & -6.222576268071369D+00, & -7.138866999945524D+00, & -7.927324360309794D+00, & -9.393661429103221D+00 / data xvec / & 0.2D+00, & 0.4D+00, & 0.6D+00, & 0.8D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 6.0D+00, & 6.0D+00, & 6.0D+00, & 6.0D+00, & 6.0D+00, & 7.0D+00 / data yvec / & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 0.2D+00, & 0.4D+00, & 1.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 6.0D+00, & 7.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 y = 0.0D+00 fxy = 0.0D+00 else x = xvec(n_data) y = yvec(n_data) fxy = fxvec(n_data) end if return end subroutine beta_noncentral_cdf_values ( n_data, a, b, lambda, & x, fx ) c*********************************************************************72 c cc beta_noncentral_cdf_values() returns some values of the noncentral Beta CDF. c c Discussion: c c The values presented here are taken from the reference, where they c were given to a limited number of decimal places. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 13 January 2008 c c Author: c c John Burkardt c c Reference: c c R Chattamvelli, R Shanmugam, c Algorithm AS 310: c Computing the Non-central Beta Distribution Function, c Applied Statistics, c Volume 46, Number 1, 1997, pages 146-156. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision A, B, the shape parameters. c c Output, double precision LAMBDA, the noncentrality parameter. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 25 ) double precision a double precision a_vec(n_max) double precision b double precision b_vec(n_max) double precision fx double precision fx_vec(n_max) double precision lambda double precision lambda_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save a_vec save b_vec save fx_vec save lambda_vec save x_vec data a_vec / & 5.0D+00, & 5.0D+00, & 5.0D+00, & 10.0D+00, & 10.0D+00, & 10.0D+00, & 20.0D+00, & 20.0D+00, & 20.0D+00, & 10.0D+00, & 10.0D+00, & 15.0D+00, & 20.0D+00, & 20.0D+00, & 20.0D+00, & 30.0D+00, & 30.0D+00, & 10.0D+00, & 10.0D+00, & 10.0D+00, & 15.0D+00, & 10.0D+00, & 12.0D+00, & 30.0D+00, & 35.0D+00 / data b_vec / & 5.0D+00, & 5.0D+00, & 5.0D+00, & 10.0D+00, & 10.0D+00, & 10.0D+00, & 20.0D+00, & 20.0D+00, & 20.0D+00, & 20.0D+00, & 10.0D+00, & 5.0D+00, & 10.0D+00, & 30.0D+00, & 50.0D+00, & 20.0D+00, & 40.0D+00, & 5.0D+00, & 10.0D+00, & 30.0D+00, & 20.0D+00, & 5.0D+00, & 17.0D+00, & 30.0D+00, & 30.0D+00 / data fx_vec / & 0.4563021D+00, & 0.1041337D+00, & 0.6022353D+00, & 0.9187770D+00, & 0.6008106D+00, & 0.0902850D+00, & 0.9998655D+00, & 0.9925997D+00, & 0.9641112D+00, & 0.9376626573D+00, & 0.7306817858D+00, & 0.1604256918D+00, & 0.1867485313D+00, & 0.6559386874D+00, & 0.9796881486D+00, & 0.1162386423D+00, & 0.9930430054D+00, & 0.0506899273D+00, & 0.1030959706D+00, & 0.9978417832D+00, & 0.2555552369D+00, & 0.0668307064D+00, & 0.0113601067D+00, & 0.7813366615D+00, & 0.8867126477D+00 / data lambda_vec / & 54.0D+00, & 140.0D+00, & 170.0D+00, & 54.0D+00, & 140.0D+00, & 250.0D+00, & 54.0D+00, & 140.0D+00, & 250.0D+00, & 150.0D+00, & 120.0D+00, & 80.0D+00, & 110.0D+00, & 65.0D+00, & 130.0D+00, & 80.0D+00, & 130.0D+00, & 20.0D+00, & 54.0D+00, & 80.0D+00, & 120.0D+00, & 55.0D+00, & 64.0D+00, & 140.0D+00, & 20.0D+00 / data x_vec / & 0.8640D+00, & 0.9000D+00, & 0.9560D+00, & 0.8686D+00, & 0.9000D+00, & 0.9000D+00, & 0.8787D+00, & 0.9000D+00, & 0.9220D+00, & 0.868D+00, & 0.900D+00, & 0.880D+00, & 0.850D+00, & 0.660D+00, & 0.720D+00, & 0.720D+00, & 0.800D+00, & 0.644D+00, & 0.700D+00, & 0.780D+00, & 0.760D+00, & 0.795D+00, & 0.560D+00, & 0.800D+00, & 0.670D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 a = 0.0D+00 b = 0.0D+00 lambda = 0.0D+00 x = 0.0D+00 fx = 0.0D+00 else a = a_vec(n_data) b = b_vec(n_data) lambda = lambda_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine beta_values ( n_data, x, y, fxy ) c*********************************************************************72 c cc beta_values() returns some values of the Beta function. c c Discussion: c c Beta(X,Y) = ( Gamma(X) * Gamma(Y) ) / Gamma(X+Y) c c Both X and Y must be greater than 0. c c In Mathematica, the function can be evaluated by: c c Beta[X,Y] c c Properties: c c Beta(X,Y) = Beta(Y,X). c Beta(X,Y) = Integral ( 0 .lt.= T .lt.= 1 ) T**(X-1) (1-T)**(Y-1) dT. c Beta(X,Y) = Gamma(X) * Gamma(Y) / Gamma(X+Y) c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 20 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision X, Y, the arguments of the function. c c Output, double precision FXY, the value of the function. c implicit none integer n_max parameter ( n_max = 17 ) double precision b_vec(n_max) double precision fxy integer n_data double precision x double precision x_vec(n_max) double precision y double precision y_vec(n_max) save b_vec save x_vec save y_vec data b_vec / & 0.5000000000000000D+01, 7 0.2500000000000000D+01, & 0.1666666666666667D+01, & 0.1250000000000000D+01, & 0.5000000000000000D+01, & 0.2500000000000000D+01, & 0.1000000000000000D+01, & 0.1666666666666667D+00, & 0.3333333333333333D-01, & 0.7142857142857143D-02, & 0.1587301587301587D-02, & 0.2380952380952381D-01, & 0.5952380952380952D-02, & 0.1984126984126984D-02, & 0.7936507936507937D-03, & 0.3607503607503608D-03, & 0.8325008325008325D-04 / data x_vec / & 0.2D+00, & 0.4D+00, & 0.6D+00, & 0.8D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 6.0D+00, & 6.0D+00, & 6.0D+00, & 6.0D+00, & 6.0D+00, & 7.0D+00 / data y_vec / & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 0.2D+00, & 0.4D+00, & 1.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 2.0D+00, & 3.0D+00, & 4.0D+00, & 5.0D+00, & 6.0D+00, & 7.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 x = 0.0D+00 y = 0.0D+00 fxy = 0.0D+00 else x = x_vec(n_data) y = y_vec(n_data) fxy = b_vec(n_data) end if return end subroutine binomial_values ( n_data, a, b, fx ) c*********************************************************************72 c cc binomial_values() returns some values of the binomial coefficients. c c Discussion: c c The formula for the binomial coefficient is c c C(N,K) = N! / ( K! * (N-K)! ) c c In Mathematica, the function can be evaluated by: c c Binomial[n,k] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 21 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer A, B, the arguments of the function. c c Output, integer FX, the value of the function. c implicit none integer n_max parameter ( n_max = 20 ) integer a integer a_vec(n_max) integer b integer b_vec(n_max) integer fx integer fx_vec(n_max) integer n_data save a_vec save b_vec save fx_vec data a_vec / & 1, 6, 6, 6, 15, & 15, 15, 15, 15, 15, & 15, 25, 25, 25, 25, & 25, 25, 25, 25, 25 / data b_vec / & 0, 1, 3, 5, 1, & 3, 5, 7, 9, 11, & 13, 1, 3, 5, 7, & 9, 11, 13, 15, 17 / data fx_vec / & 1, & 6, & 20, & 6, & 15, & 455, & 3003, & 6435, & 5005, & 1365, & 105, & 25, & 2300, & 53130, & 480700, & 2042975, & 4457400, & 5200300, & 3268760, & 1081575 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 a = 0 b = 0 fx = 0 else a = a_vec(n_data) b = b_vec(n_data) fx = fx_vec(n_data) end if return end subroutine binomial_cdf_values ( n_data, a, b, x, fx ) c*********************************************************************72 c cc binomial_cdf_values() returns some values of the binomial CDF. c c Discussion: c c CDF(X)(A,B) is the probability of at most X successes in A trials, c given that the probability of success on a single trial is B. c c In Mathematica, the function can be evaluated by: c c Needs["Statistics`DiscreteDistributions`] c dist = BinomialDistribution [ n, p ] c CDF [ dist, x ] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 21 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Daniel Zwillinger, editor, c CRC Standard Mathematical Tables and Formulae, c 30th Edition, c CRC Press, 1996, c ISBN: 0-8493-2479-3, c LC: QA47.M315. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer A, a parameter of the function. c c Output, double precision B, a parameter of the function. c c Output, integer X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 17 ) integer a integer a_vec(n_max) double precision b double precision b_vec(n_max) double precision fx double precision fx_vec(n_max) integer n_data integer x integer x_vec(n_max) save a_vec save b_vec save fx_vec save x_vec data a_vec / & 2, 2, 2, 2, & 2, 4, 4, 4, & 4, 10, 10, 10, & 10, 10, 10, 10, & 10 / data b_vec / & 0.05D+00, & 0.05D+00, & 0.05D+00, & 0.50D+00, & 0.50D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.25D+00, & 0.05D+00, & 0.10D+00, & 0.15D+00, & 0.20D+00, & 0.25D+00, & 0.30D+00, & 0.40D+00, & 0.50D+00 / data fx_vec / & 0.9025000000000000D+00, & 0.9975000000000000D+00, & 0.1000000000000000D+01, & 0.2500000000000000D+00, & 0.7500000000000000D+00, & 0.3164062500000000D+00, & 0.7382812500000000D+00, & 0.9492187500000000D+00, & 0.9960937500000000D+00, & 0.9999363101685547D+00, & 0.9983650626000000D+00, & 0.9901259090013672D+00, & 0.9672065024000000D+00, & 0.9218730926513672D+00, & 0.8497316674000000D+00, & 0.6331032576000000D+00, & 0.3769531250000000D+00 / data x_vec / & 0, 1, 2, 0, & 1, 0, 1, 2, & 3, 4, 4, 4, & 4, 4, 4, 4, & 4 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 a = 0 b = 0.0D+00 x = 0 fx = 0.0D+00 else a = a_vec(n_data) b = b_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine bivariate_normal_cdf_values ( n_data, x, y, r, fxy ) c*********************************************************************72 c cc bivariate_normal_cdf_values() returns some values of the bivariate normal CDF. c c Discussion: c c FXY is the probability that two variables A and B, which are c related by a bivariate normal distribution with correlation R, c respectively satisfy A .lt.= X and B .lt.= Y. c c Mathematica can evaluate the bivariate normal CDF via the commands: c c < .0 = 0.0D+00 c 1 = 1 => .1 = 0.5 c 2 = 10 => .01 = 0.25 c 3 = 11 => .11 = 0.75 c 4 = 100 => .001 = 0.125 c 5 = 101 => .101 = 0.625 c 6 = 110 => .011 = 0.375 c 7 = 111 => .111 = 0.875 c 8 = 1000 => .0001 = 0.0625 c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 25 March 2007 c c Author: c c John Burkardt c c Reference: c c John Halton, c On the efficiency of certain quasi-random sequences of points c in evaluating multi-dimensional integrals, c Numerische Mathematik, c Volume 2, pages 84-90, 1960. c c John Hammersley, c Monte Carlo methods for solving multivariable problems, c Proceedings of the New York Academy of Science, c Volume 86, pages 844-874, 1960. c c Johannes van der Corput, c Verteilungsfunktionen, c Proc Akad Amsterdam, c Volume 38, 1935, c Volume 39, 1936. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer BASE, the base of the sequence. c c Output, integer SEED, the index of the element of the sequence. c c Output, double precision VALUE, the value of the SEED-th element of the c van der Corput sequence in base BASE. c implicit none integer n_max parameter ( n_max = 75 ) integer base integer base_vec(n_max) integer n_data integer seed integer seed_vec(n_max) double precision value double precision value_vec(n_max) save base_vec save seed_vec save value_vec data base_vec / & 2, 2, 2, 2, 2, & 2, 2, 2, 2, 3, & 3, 3, 3, 3, 3, & 3, 3, 3, 4, 4, & 4, 4, 4, 4, 4, & 4, 4, 2, 3, 4, & 5, 7, 11, 13, 2, & 3, 4, 5, 7, 11, & 13, 2, 3, 4, 5, & 7, 11, 13, 2, 3, & 4, 5, 7, 11, 13, & 29, 29, 29, 29, 29, & 71, 71, 71, 71, 71, & 173, 173, 173, 173, 173, & 409, 409, 409, 409, 409 / data seed_vec / & 0, 1, 2, 3, 4, & 5, 6, 7, 8, 0, & 1, 2, 3, 4, 5, & 6, 7, 8, 0, 1, & 2, 3, 4, 5, 6, & 7, 8, 10, 10, 10, & 10, 10, 10, 10, 100, & 100, 100, 100, 100, 100, & 100, 1000, 1000, 1000, 1000, & 1000, 1000, 1000, 10000, 10000, & 10000, 10000, 10000, 10000, 10000, & 1000, 1001, 1002, 1003, 1004, & 1000, 1001, 1002, 1003, 1004, & 1000, 1001, 1002, 1003, 1004, & 1000, 1001, 1002, 1003, 1004 / data value_vec / & 0.0000000000000000D+00, & 0.5000000000000000D+00, & 0.2500000000000000D+00, & 0.7500000000000000D+00, & 0.1250000000000000D+00, & 0.6250000000000000D+00, & 0.3750000000000000D+00, & 0.8750000000000000D+00, & 0.0625000000000000D+00, & 0.0000000000000000D+00, & 0.3333333333333333D+00, & 0.6666666666666666D+00, & 0.1111111111111111D+00, & 0.4444444444444444D+00, & 0.7777777777777777D+00, & 0.2222222222222222D+00, & 0.5555555555555556D+00, & 0.8888888888888888D+00, & 0.0000000000000000D+00, & 0.2500000000000000D+00, & 0.5000000000000000D+00, & 0.7500000000000000D+00, & 0.0625000000000000D+00, & 0.3125000000000000D+00, & 0.5625000000000000D+00, & 0.8125000000000000D+00, & 0.1250000000000000D+00, & 0.3125000000000000D+00, & 0.3703703703703703D+00, & 0.6250000000000000D+00, & 0.0800000000000000D+00, & 0.4489795918367347D+00, & 0.9090909090909092D+00, & 0.7692307692307693D+00, & 0.1484375000000000D+00, & 0.4115226337448559D+00, & 0.0976562500000000D+00, & 0.0320000000000000D+00, & 0.2915451895043731D+00, & 0.1652892561983471D+00, & 0.7337278106508875D+00, & 0.0927734375000000D+00, & 0.3475080018289895D+00, & 0.1708984375000000D+00, & 0.0051200000000000D+00, & 0.9162848812994586D+00, & 0.9316303531179565D+00, & 0.9904415111515704D+00, & 0.0347290039062500D+00, & 0.3861200020322105D+00, & 0.0189208984375000D+00, & 0.0005120000000000D+00, & 0.5749985125245433D+00, & 0.1529950140017758D+00, & 0.2459297643639929D+00, & 0.4887449259912255D+00, & 0.5232276846119153D+00, & 0.5577104432326049D+00, & 0.5921932018532945D+00, & 0.6266759604739842D+00, & 0.0872842689942472D+00, & 0.1013687760365007D+00, & 0.1154532830787542D+00, & 0.1295377901210077D+00, & 0.1436222971632613D+00, & 0.7805138828560928D+00, & 0.7862942296769020D+00, & 0.7920745764977113D+00, & 0.7978549233185205D+00, & 0.8036352701393298D+00, & 0.4449997309915651D+00, & 0.4474447187666262D+00, & 0.4498897065416874D+00, & 0.4523346943167484D+00, & 0.4547796820918096D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 base = 0 seed = 0 value = 0.0D+00 else base = base_vec(n_data) seed = seed_vec(n_data) value = value_vec(n_data) end if return end subroutine viscosity_values ( n_data, tc, p, eta ) c*********************************************************************72 c cc VISCOSITY_VALUES returns some values of the viscosity function. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 March 2007 c c Author: c c John Burkardt c c Reference: c c Lester Haar, John Gallagher, George Kell, c NBS/NRC Steam Tables: c Thermodynamic and Transport Properties and Computer Programs c for Vapor and Liquid States of Water in SI Units, c Hemisphere Publishing Corporation, Washington, 1984, c ISBN: 0-89116-353-0, c LC: TJ270.H3. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision TC, the temperature, in degrees Celsius. c c Output, double precision P, the pressure, in bar. c c Output, double precision ETA, the viscosity, in MegaPascal seconds. c implicit none integer n_max parameter ( n_max = 34 ) double precision eta double precision eta_vec(n_max) integer n_data double precision p double precision p_vec(n_max) double precision tc double precision tc_vec(n_max) save eta_vec save p_vec save tc_vec data eta_vec / & 1792.0D+00, & 1791.0D+00, & 1790.0D+00, & 1786.0D+00, & 1780.0D+00, & 1775.0D+00, & 1769.0D+00, & 1764.0D+00, & 1759.0D+00, & 1754.0D+00, & 1749.0D+00, & 1744.0D+00, & 1739.0D+00, & 1735.0D+00, & 1731.0D+00, & 1722.0D+00, & 1714.0D+00, & 1707.0D+00, & 1700.0D+00, & 1694.0D+00, & 1687.0D+00, & 1682.0D+00, & 1676.0D+00, & 1667.0D+00, & 1659.0D+00, & 1653.0D+00, & 890.8D+00, & 547.1D+00, & 378.4D+00, & 12.28D+00, & 16.18D+00, & 24.45D+00, & 32.61D+00, & 40.38D+00 / data p_vec / & 1.0D+00, & 5.0D+00, & 10.0D+00, & 25.0D+00, & 50.0D+00, & 75.0D+00, & 100.0D+00, & 125.0D+00, & 150.0D+00, & 175.0D+00, & 200.0D+00, & 225.0D+00, & 250.0D+00, & 275.0D+00, & 300.0D+00, & 350.0D+00, & 400.0D+00, & 450.0D+00, & 500.0D+00, & 550.0D+00, & 600.0D+00, & 650.0D+00, & 700.0D+00, & 800.0D+00, & 900.0D+00, & 1000.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00, & 1.0D+00 / data tc_vec / & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 25.0D+00, & 50.0D+00, & 75.0D+00, & 100.0D+00, & 200.0D+00, & 400.0D+00, & 600.0D+00, & 800.0D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 tc = 0.0D+00 p = 0.0D+00 eta = 0.0D+00 else tc = tc_vec(n_data) p = p_vec(n_data) eta = eta_vec(n_data) end if return end subroutine von_mises_cdf_values ( n_data, a, b, x, fx ) c*********************************************************************72 c cc VON_MISES_CDF_VALUES returns some values of the von Mises CDF. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 March 2007 c c Author: c c John Burkardt c c Reference: c c Kanti Mardia, Peter Jupp, c Directional Statistics, c Wiley, 2000, c LC: QA276.M335 c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision A, a parameter of the PDF. c A is the preferred direction, in radians. c -PI <= A <= PI. c c Output, double precision B, a parameter of the PDF. c B measures the "concentration" of the distribution around the c angle A. B = 0 corresponds to a uniform distribution c (no concentration). Higher values of B cause greater concentration c of probability near A. c 0.0D+00 <= B. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 23 ) double precision a double precision a_vec(n_max) double precision b double precision b_vec(n_max) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save a_vec save b_vec save fx_vec save x_vec data a_vec / & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & -0.2D+01, & -0.1D+01, & 0.0D+01, & 0.1D+01, & 0.2D+01, & 0.3D+01, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00, & 0.0D+00 / data b_vec / & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.1D+01, & 0.2D+01, & 0.2D+01, & 0.2D+01, & 0.2D+01, & 0.2D+01, & 0.2D+01, & 0.3D+01, & 0.3D+01, & 0.3D+01, & 0.3D+01, & 0.3D+01, & 0.3D+01, & 0.0D+00, & 0.1D+01, & 0.2D+01, & 0.3D+01, & 0.4D+01, & 0.5D+01 / data fx_vec / & 0.2535089956281180D-01, & 0.1097539041177346D+00, & 0.5000000000000000D+00, & 0.8043381312498558D+00, & 0.9417460124555197D+00, & 0.5000000000000000D+00, & 0.6018204118446155D+00, & 0.6959356933122230D+00, & 0.7765935901304593D+00, & 0.8410725934916615D+00, & 0.8895777369550366D+00, & 0.9960322705517925D+00, & 0.9404336090170247D+00, & 0.5000000000000000D+00, & 0.5956639098297530D-01, & 0.3967729448207649D-02, & 0.2321953958111930D-03, & 0.6250000000000000D+00, & 0.7438406999109122D+00, & 0.8369224904294019D+00, & 0.8941711407897124D+00, & 0.9291058600568743D+00, & 0.9514289900655436D+00 / data x_vec / & -0.2617993977991494D+01, & -0.1570796326794897D+01, & 0.0000000000000000D+00, & 0.1047197551196598D+01, & 0.2094395102393195D+01, & 0.1000000000000000D+01, & 0.1200000000000000D+01, & 0.1400000000000000D+01, & 0.1600000000000000D+01, & 0.1800000000000000D+01, & 0.2000000000000000D+01, & 0.0000000000000000D+00, & 0.0000000000000000D+00, & 0.0000000000000000D+00, & 0.0000000000000000D+00, & 0.0000000000000000D+00, & 0.0000000000000000D+00, & 0.7853981633974483D+00, & 0.7853981633974483D+00, & 0.7853981633974483D+00, & 0.7853981633974483D+00, & 0.7853981633974483D+00, & 0.7853981633974483D+00 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 a = 0.0D+00 b = 0.0D+00 x = 0.0D+00 fx = 0.0D+00 else a = a_vec(n_data) b = b_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine weekday_values ( n_data, y, m, d, w ) c*********************************************************************72 c cc WEEKDAY_VALUES returns the day of the week for various dates. c c Discussion: c c The CE or Common Era calendar is used, under the c hybrid Julian/Gregorian Calendar, with a transition from Julian c to Gregorian. The day after 04 October 1582 was 15 October 1582. c c The year before 1 AD or CE is 1 BC or BCE. In this data set, c years BC/BCE are indicated by a negative year value. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 26 May 2012 c c Author: c c John Burkardt c c Reference: c c Edward Reingold, Nachum Dershowitz, c Calendrical Calculations: The Millennium Edition, c Cambridge University Press, 2001, c ISBN: 0 521 77752 6 c LC: CE12.R45. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 c before the first call. On each call, the routine increments N_DATA by 1, c and returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer Y, M, D, the Common Era date. c c Output, integer W, the day of the week. Sunday = 1. c implicit none integern_max parameter ( n_max = 34 ) integer d integer d_vec(n_max) integer m integer m_vec(n_max) integer n_data integer w integer w_vec(n_max) integer y integer y_vec(n_max) save d_vec save m_vec save w_vec save y_vec data d_vec / & 30, & 8, & 26, & 3, & 7, & 18, & 7, & 19, & 14, & 18, & 16, & 3, & 26, & 20, & 4, & 25, & 31, & 9, & 24, & 10, & 30, & 24, & 19, & 2, & 27, & 19, & 25, & 29, & 19, & 7, & 17, & 25, & 10, & 18 / data m_vec / & 7, & 12, & 9, & 10, & 1, & 5, & 11, & 4, & 10, & 5, & 3, & 3, & 3, & 4, & 6, & 1, & 3, & 9, & 2, & 6, & 6, & 7, & 6, & 8, & 3, & 4, & 8, & 9, & 4, & 10, & 3, & 2, & 11, & 7 / data w_vec / & 1, & 4, & 4, & 1, & 4, & 2, & 7, & 1, & 7, & 1, & 6, & 7, & 6, & 1, & 1, & 4, & 7, & 7, & 7, & 4, & 1, & 6, & 1, & 2, & 4, & 1, & 1, & 2, & 2, & 5, & 3, & 1, & 4, & 1 / data y_vec / & - 587, & - 169, & 70, & 135, & 470, & 576, & 694, & 1013, & 1066, & 1096, & 1190, & 1240, & 1288, & 1298, & 1391, & 1436, & 1492, & 1553, & 1560, & 1648, & 1680, & 1716, & 1768, & 1819, & 1839, & 1903, & 1929, & 1941, & 1943, & 1943, & 1992, & 1996, & 2038, & 2094 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 y = 0 m = 0 d = 0 w = 0 else y = y_vec(n_data) m = m_vec(n_data) d = d_vec(n_data) w = w_vec(n_data) end if return end subroutine weibull_cdf_values ( n_data, alpha, beta, x, fx ) c*********************************************************************72 c cc WEIBULL_CDF_VALUES returns some values of the Weibull CDF. c c Discussion: c c In Mathematica, the function can be evaluated by: c c Needs["Statistics`ContinuousDistributions`"] c dist = WeibullDistribution [ alpha, beta ] c CDF [ dist, x ] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, double precision ALPHA, the first parameter of the distribution. c c Output, double precision BETA, the second parameter of the distribution. c c Output, double precision X, the argument of the function. c c Output, double precision FX, the value of the function. c implicit none integer n_max parameter ( n_max = 12 ) double precision alpha double precision alpha_vec(n_max) double precision beta double precision beta_vec(n_max) double precision fx double precision fx_vec(n_max) integer n_data double precision x double precision x_vec(n_max) save alpha_vec save beta_vec save fx_vec save x_vec data alpha_vec / & 0.1000000000000000D+01, & 0.1000000000000000D+01, & 0.1000000000000000D+01, & 0.1000000000000000D+01, & 0.1000000000000000D+01, & 0.1000000000000000D+01, & 0.1000000000000000D+01, & 0.1000000000000000D+01, & 0.2000000000000000D+01, & 0.3000000000000000D+01, & 0.4000000000000000D+01, & 0.5000000000000000D+01 / data beta_vec / & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.2000000000000000D+01, & 0.3000000000000000D+01, & 0.4000000000000000D+01, & 0.5000000000000000D+01, & 0.2000000000000000D+01, & 0.2000000000000000D+01, & 0.2000000000000000D+01, & 0.2000000000000000D+01 / data fx_vec / & 0.8646647167633873D+00, & 0.9816843611112658D+00, & 0.9975212478233336D+00, & 0.9996645373720975D+00, & 0.6321205588285577D+00, & 0.4865828809674080D+00, & 0.3934693402873666D+00, & 0.3296799539643607D+00, & 0.8946007754381357D+00, & 0.9657818816883340D+00, & 0.9936702845725143D+00, & 0.9994964109502630D+00 / data x_vec / & 0.1000000000000000D+01, & 0.2000000000000000D+01, & 0.3000000000000000D+01, & 0.4000000000000000D+01, & 0.2000000000000000D+01, & 0.2000000000000000D+01, & 0.2000000000000000D+01, & 0.2000000000000000D+01, & 0.3000000000000000D+01, & 0.3000000000000000D+01, & 0.3000000000000000D+01, & 0.3000000000000000D+01 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 alpha = 0.0D+00 beta = 0.0D+00 x = 0.0D+00 fx = 0.0D+00 else alpha = alpha_vec(n_data) beta = beta_vec(n_data) x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine zeta_values ( n_data, n, zeta ) c*********************************************************************72 c cc ZETA_VALUES returns some values of the Riemann Zeta function. c c Discussion: c c ZETA(N) = sum ( 1 <= I < +oo ) 1 / I^N c c In Mathematica, the function can be evaluated by: c c Zeta[n] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 March 2007 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Parameters: c c Input/output, integer N_DATA. The user sets N_DATA to 0 before the c first call. On each call, the routine increments N_DATA by 1, and c returns the corresponding data; when there is no more data, the c output value of N_DATA will be 0 again. c c Output, integer N, the argument of the Zeta function. c c Output, double precision ZETA, the value of the Zeta function. c implicit none integer n_max parameter ( n_max = 15 ) integer n integer n_data integer n_vec(n_max) double precision zeta double precision zeta_vec(n_max) save n_vec save zeta_vec data n_vec / & 2, & 3, & 4, & 5, & 6, & 7, & 8, & 9, & 10, & 11, & 12, & 16, & 20, & 30, & 40 / data zeta_vec / & 0.164493406684822643647D+01, & 0.120205690315959428540D+01, & 0.108232323371113819152D+01, & 0.103692775514336992633D+01, & 0.101734306198444913971D+01, & 0.100834927738192282684D+01, & 0.100407735619794433939D+01, & 0.100200839292608221442D+01, & 0.100099457512781808534D+01, & 0.100049418860411946456D+01, & 0.100024608655330804830D+01, & 0.100001528225940865187D+01, & 0.100000095396203387280D+01, & 0.100000000093132743242D+01, & 0.100000000000090949478D+01 / if ( n_data .lt. 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max .lt. n_data ) then n_data = 0 n = 0 zeta = 0.0D+00 else n = n_vec(n_data) zeta = zeta_vec(n_data) end if return end