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