#! /usr/bin/env python3 # def test_values_test ( ): #*****************************************************************************80 # ## test_values_test() tests test_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 January 2024 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'test_values_test():' ) print ( ' python version: ' + platform.python_version ( ) ) print ( ' numpy version: ' + np.version.version ) print ( ' Test test_values().' ) abram0_values_test ( ) abram1_values_test ( ) abram2_values_test ( ) agm_values_test ( ) airy_ai_values_test ( ) airy_ai_int_values_test ( ) airy_ai_prime_values_test ( ) airy_bi_values_test ( ) airy_bi_int_values_test ( ) airy_bi_prime_values_test ( ) airy_cai_values_test ( ) airy_cbi_values_test ( ) airy_gi_values_test ( ) arccos_values_test ( ) arccosh_values_test ( ) arcsin_values_test ( ) arcsinh_values_test ( ) arctan_values_test ( ) arctan_int_values_test ( ) arctan2_values_test ( ) arctanh_values_test ( ) bei0_values_test ( ) bei1_values_test ( ) bell_values_test ( ) ber0_values_test ( ) ber1_values_test ( ) bernoulli_number_values_test ( ) bernoulli_poly_values_test ( ) bernstein_poly_01_values_test ( ) bessel_i0_values_test ( ) bessel_i0_int_values_test ( ) bessel_i0_spherical_values_test ( ) bessel_i1_values_test ( ) bessel_i1_spherical_values_test ( ) bessel_in_values_test ( ) bessel_ix_values_test ( ) bessel_j_spherical_values_test ( ) bessel_j0_values_test ( ) bessel_j0_int_values_test ( ) bessel_j0_spherical_values_test ( ) bessel_j0_zero_values_test ( ) bessel_j1_values_test ( ) bessel_j1_spherical_values_test ( ) bessel_jn_values_test ( ) bessel_jx_values_test ( ) bessel_k0_values_test ( ) bessel_k0_int_values_test ( ) bessel_k1_values_test ( ) bessel_kn_values_test ( ) bessel_kx_values_test ( ) bessel_y0_values_test ( ) bessel_y0_int_values_test ( ) bessel_y0_spherical_values_test ( ) bessel_y0_zero_values_test ( ) bessel_y1_values_test ( ) bessel_y1_spherical_values_test ( ) bessel_yn_values_test ( ) bessel_yx_values_test ( ) beta_values_test ( ) beta_cdf_values_test ( ) beta_inc_values_test ( ) beta_log_values_test ( ) beta_noncentral_cdf_values_test ( ) beta_pdf_values_test ( ) binomial_values_test ( ) binomial_cdf_values_test ( ) binomial_pdf_values_test ( ) bivariate_normal_cdf_values_test ( ) c8_log_values_test ( ) catalan_values_test ( ) cauchy_cdf_values_test ( ) cbrt_values_test ( ) cheby_t_poly_values_test ( ) cheby_t01_poly_values_test ( ) cheby_u_poly_values_test ( ) cheby_u01_poly_values_test ( ) cheby_v_poly_values_test ( ) cheby_v01_poly_values_test ( ) cheby_w_poly_values_test ( ) cheby_w01_poly_values_test ( ) chi_values_test ( ) chi_square_cdf_values_test ( ) chi_square_pdf_values_test ( ) chi_square_noncentral_cdf_values_test ( ) ci_values_test ( ) cin_values_test ( ) cinh_values_test ( ) clausen_values_test ( ) clebsch_gordan_values_test ( ) collatz_count_values_test ( ) cos_values_test ( ) cos_degree_values_test ( ) cos_power_int_values_test ( ) cosh_values_test ( ) cot_values_test ( ) cp_values_test ( ) datenum_values_test ( ) dawson_values_test ( ) debye1_values_test ( ) debye2_values_test ( ) debye3_values_test ( ) debye4_values_test ( ) dedekind_sum_values_test ( ) dielectric_values_test ( ) dilogarithm_values_test ( ) dixon_elliptic_values_test ( ) e1_values_test ( ) easter_gregorian_values_test ( ) easter_julian_values_test ( ) ei_values_test ( ) elliptic_ea_values_test ( ) elliptic_ek_values_test ( ) elliptic_em_values_test ( ) elliptic_fa_values_test ( ) elliptic_fk_values_test ( ) elliptic_fm_values_test ( ) elliptic_pia_values_test ( ) elliptic_pik_values_test ( ) elliptic_pim_values_test ( ) elliptic_inc_ea_values_test ( ) elliptic_inc_ek_values_test ( ) elliptic_inc_em_values_test ( ) elliptic_inc_fa_values_test ( ) elliptic_inc_fk_values_test ( ) elliptic_inc_fm_values_test ( ) elliptic_inc_pia_values_test ( ) elliptic_inc_pik_values_test ( ) elliptic_inc_pim_values_test ( ) erf_values_test ( ) erfc_values_test ( ) euler_number_values_test ( ) euler_poly_values_test ( ) exp_values_test ( ) exp3_int_values_test ( ) exponential_01_pdf_values_test ( ) exponential_cdf_values_test ( ) exponential_pdf_values_test ( ) extreme_values_cdf_values_test ( ) f_cdf_values_test ( ) f_noncentral_cdf_values_test ( ) factorial_values_test ( ) factorial2_values_test ( ) fresnel_cos_values_test ( ) fresnel_sin_values_test ( ) frobenius_number_data_values_test ( ) frobenius_number_order_values_test ( ) frobenius_number_order2_values_test ( ) gamma_values_test ( ) gamma_01_pdf_values_test ( ) gamma_inc_values_test ( ) gamma_inc_p_values_test ( ) gamma_inc_q_values_test ( ) gamma_log_values_test ( ) gamma_pdf_values_test ( ) gaussian_prime_complex_values_test ( ) gaussian_prime_real_values_test ( ) gcd_values_test ( ) gegenbauer_polynomial_values_test ( ) geometric_cdf_values_test ( ) goodwin_values_test ( ) gud_values_test ( ) harmonic_values_test ( ) hermite_function_values_test ( ) hermite_poly_phys_values_test ( ) hermite_poly_prob_values_test ( ) hyper_1f1_values_test ( ) hyper_2f1_values_test ( ) hyper_2f1_complex_values_test ( ) hypergeometric_cdf_values_test ( ) hypergeometric_pdf_values_test ( ) hypergeometric_u_values_test ( ) i0ml0_values_test ( ) i1ml1_values_test ( ) i4_fall_values_test ( ) i4_gpf_values_test ( ) i4_rise_values_test ( ) int_values_test ( ) inverse_chi_square_pdf_values_test ( ) inverse_gamma_pdf_values_test ( ) is_gaussian_prime_values_test ( ) is_prime_values_test ( ) jacobi_cn_values_test ( ) jacobi_dn_values_test ( ) jacobi_poly_values_test ( ) jacobi_sn_values_test ( ) jed_ce_values_test ( ) jed_mjd_values_test ( ) jed_rd_values_test ( ) jed_weekday_values_test ( ) kei0_values_test ( ) kei1_values_test ( ) ker0_values_test ( ) ker1_values_test ( ) laguerre_associated_values_test ( ) laguerre_general_values_test ( ) laguerre_polynomial_values_test ( ) lambert_w_values_test ( ) laplace_cdf_values_test ( ) legendre_associated_values_test ( ) legendre_associated_normalized_values_test ( ) legendre_associated_normalized_sphere_values_test ( ) legendre_function_q_values_test ( ) legendre_normalized_polynomial_values_test ( ) legendre_polynomial_values_test ( ) legendre_shifted_polynomial_values_test ( ) lerch_values_test ( ) lobachevsky_values_test ( ) lobatto_polynomial_values_test ( ) lobatto_polynomial_derivative_values_test ( ) log_values_test ( ) log_normal_cdf_values_test ( ) log_series_cdf_values_test ( ) log10_values_test ( ) logarithmic_integral_values_test ( ) logistic_cdf_values_test ( ) mcnugget_number_values_test ( ) mersenne_prime_values_test ( ) mertens_values_test ( ) mittag_leffler_ea_values_test ( ) moebius_values_test ( ) multinomial_pdf_values_test ( ) negative_binomial_cdf_values_test ( ) nine_j_values_test ( ) normal_01_cdf_values_test ( ) normal_cdf_values_test ( ) normal_pdf_values_test ( ) omega_values_test ( ) owen_values_test ( ) partition_count_values_test ( ) partition_distinct_count_values_test ( ) phi_values_test ( ) pi_values_test ( ) poisson_cdf_values_test ( ) polylogarithm_values_test ( ) polynomial_resultant_values_test ( ) polyomino_chiral_count_values_test ( ) polyomino_fixed_count_values_test ( ) polyomino_free_count_values_test ( ) prandtl_values_test ( ) prime_values_test ( ) psat_values_test ( ) psi_values_test ( ) r8_factorial_values_test ( ) r8_factorial_log_values_test ( ) r8_factorial2_values_test ( ) r8_fall_values_test ( ) r8_rise_values_test ( ) rayleigh_cdf_values_test ( ) scaled_inverse_chi_square_pdf_values_test ( ) secvir_values_test ( ) shi_values_test ( ) si_values_test ( ) sigma_values_test ( ) sin_values_test ( ) sin_degree_values_test ( ) sin_power_int_values_test ( ) sinh_values_test ( ) six_j_values_test ( ) sound_values_test ( ) sphere_unit_area_values_test ( ) sphere_unit_volume_values_test ( ) spherical_harmonic_values_test ( ) sqrt_values_test ( ) stirling1_values_test ( ) stirling2_values_test ( ) stromgen_values_test ( ) struve_h0_values_test ( ) struve_h1_values_test ( ) struve_l0_values_test ( ) struve_l1_values_test ( ) student_cdf_values_test ( ) student_noncentral_cdf_values_test ( ) subfactorial_values_test ( ) surten_values_test ( ) synch1_values_test ( ) synch2_values_test ( ) tan_values_test ( ) tanh_values_test ( ) tau_values_test ( ) thercon_values_test ( ) three_j_values_test ( ) tran02_values_test ( ) tran03_values_test ( ) tran04_values_test ( ) tran05_values_test ( ) tran06_values_test ( ) tran07_values_test ( ) tran08_values_test ( ) tran09_values_test ( ) trigamma_values_test ( ) truncated_normal_ab_cdf_values_test ( ) truncated_normal_ab_pdf_values_test ( ) truncated_normal_a_cdf_values_test ( ) truncated_normal_a_pdf_values_test ( ) truncated_normal_b_cdf_values_test ( ) truncated_normal_b_pdf_values_test ( ) tsat_values_test ( ) van_der_corput_values_test ( ) viscosity_values_test ( ) von_mises_cdf_values_test ( ) weekday_values_test ( ) weibull_cdf_values_test ( ) wright_omega_values_test ( ) zeta_values_test ( ) zeta_m1_values_test ( ) # # Terminate. # print ( '' ) print ( 'test_values_test():' ) print ( ' Normal end of execution.' ) return def abram0_values ( n_data ): #*****************************************************************************80 # ## abram0_values() returns some values of the Abramowitz0 function. # # Discussion: # # The function is defined by: # # abram0(X) = integral ( 0 <= T < infinity ) exp ( -T * T - X / T ) dT # # The data was reported by McLeod. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 17 November 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Allan McLeod, # Algorithm 757, MISCFUN: A software package to compute uncommon # special functions, # ACM transactions on Mathematical Software, # Volume 22, Number 3, September 1996, pages 288-301. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def abram0_values_test ( ): #*****************************************************************************80 # ## abram0_values_test() tests abram0_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 13 April 2007 # # Author: # # John Burkardt # print ( '' ) print ( 'abram0_values_test():' ) print ( ' abram0_values() stores values of' ) print ( ' the Abramowitz function of order 0.' ) print ( '' ) print ( ' X FX' ) print ( '' ) n_data = 0 while ( True ): [ n_data, x, fx ] = abram0_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16f' % ( x, fx ) ) return def abram1_values ( n_data ): #*****************************************************************************80 # ## abram1_values() returns some values of the Abramowitz1 function. # # Discussion: # # The function is defined by: # # abram1(x) = integral ( 0 <= t < infinity ) t * exp ( -t^2 - x / t ) dt # # The data was reported by McLeod. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 15 September 2004 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Allan McLeod, # Algorithm 757, MISCFUN: A software package to compute uncommon # special functions, # ACM transactions on Mathematical Software, # Volume 22, Number 3, September 1996, pages 288-301. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def abram1_values_test ( ): #*****************************************************************************80 # ## abram1_values_test() tests abram1_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 13 April 2007 # # Author: # # John Burkardt # print ( '' ) print ( 'abram1_values_test():' ) print ( ' abram1_values() stores values of' ) print ( ' the Abramowitz function of order 1.' ) print ( '' ) print ( ' X FX' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = abram1_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16f' % ( x, fx ) ) return def abram2_values ( n_data ): #*****************************************************************************80 # ## abram2_values() returns some values of the Abramowitz2 function. # # Discussion: # # The function is defined by: # # abram2(x) = Integral ( 0 <= t < infinity ) t^2 * exp( -t^2 - x / t ) dt # # The data was reported by McLeod. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 November 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Allan McLeod, # Algorithm 757, MISCFUN: A software package to compute uncommon # special functions, # ACM transactions on Mathematical Software, # Volume 22, Number 3, September 1996, pages 288-301. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def abram2_values_test ( ): #*****************************************************************************80 # ## abram2_values_test() tests abram2_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 13 April 2007 # # Author: # # John Burkardt # print ( '' ) print ( 'abram2_values_test():' ) print ( ' abram2_values() stores values of' ) print ( ' the Abramowitz function of order 2.' ) print ( '' ) print ( ' X FX' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = abram2_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16f' % ( x, fx ) ) return def agm_values ( n_data ): #*****************************************************************************80 # ## 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: # # 19 November 2014 # # Author: # # John Burkardt # # Reference: # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real A, B, the argument ofs the function. # # real FX, the value of the function. # import numpy as np n_max = 14 a_vec = np.array ( ( \ 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 ) ) b_vec = np.array ( ( \ 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 ) ) fx_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 a = 0.0 b = 0.0 fx = 0.0 else: a = a_vec[n_data] b = b_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, a, b, fx def agm_values_test ( ): #*****************************************************************************80 # ## agm_values_test() tests agm_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 09 February 2008 # # Author: # # John Burkardt # print ( '' ) print ( 'agm_values_test():' ) print ( ' agm_values() stores values of' ) print ( ' the arithmetic geometric mean function.' ) print ( '' ) print ( ' A B AGM(A,B)' ) print ( '' ) n_data = 0 while ( True ): n_data, a, b, fx = agm_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %24.16f' % ( a, b, fx ) ) return def airy_ai_int_values ( n_data ): #*****************************************************************************80 # ## 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: # # 05 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Allan McLeod, # Algorithm 757, MISCFUN: A software package to compute uncommon # special functions, # ACM transactions on Mathematical Software, # Volume 22, Number 3, September 1996, pages 288-301. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ -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 ) ) x_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data =n_data + 1 return n_data, x, fx def airy_ai_int_values_test ( ): #*****************************************************************************80 # ## airy_ai_int_values_test() tests airy_ai_int_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 05 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'airy_ai_int_values_test():' ) print ( ' airy_ai_int_values() stores values of' ) print ( ' the integral of the Airy Ai function.' ) print ( '' ) print ( ' X FX' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = airy_ai_int_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16f' % ( x, fx ) ) return def airy_ai_prime_values ( n_data ): #*****************************************************************************80 # ## 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: # # 15 September 2004 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real AIP, the derivative of the Airy AI function. # import numpy as np n_max = 11 fx_vec = np.array ( ( \ -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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data =n_data + 1 return n_data, x, fx def airy_ai_prime_values_test ( ): #*****************************************************************************80 # ## airy_ai_prime_values_test() tests airy_ai_prime_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'airy_ai_prime_values_test():' ) print ( ' airy_ai_prime_values() stores values of' ) print ( ' the derivative of the Airy Ai function.' ) print ( '' ) print ( ' X FX' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = airy_ai_prime_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16f' % ( x, fx ) ) return def airy_ai_values ( n_data ): #*****************************************************************************80 # ## 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: # # 15 September 2004 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real AI, the value of the Airy AI function. # import numpy as np n_max = 11 ai_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 ai = 0.0 else: x = x_vec[n_data] ai = ai_vec[n_data] n_data =n_data + 1 return n_data, x, ai def airy_ai_values_test ( ): #*****************************************************************************80 # ## airy_ai_values_test() tests airy_ai_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 13 April 2007 # # Author: # # John Burkardt # print ( '' ) print ( 'airy_ai_values_test():' ) print ( ' airy_ai_values() stores values of' ) print ( ' the Airy function Ai(x).' ) print ( '' ) print ( ' X Ai(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, ai = airy_ai_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16f' % ( x, ai ) ) return def airy_bi_int_values ( n_data ): #*****************************************************************************80 # ## 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: # # 14 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Allan McLeod, # Algorithm 757, MISCFUN: A software package to compute uncommon # special functions, # ACM transactions on Mathematical Software, # Volume 22, Number 3, September 1996, pages 288-301. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data =n_data + 1 return n_data, x, fx def airy_bi_int_values_test ( ): #*****************************************************************************80 # ## airy_bi_int_values_test() tests airy_bi_int_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 14 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'airy_bi_int_values_test():' ) print ( ' airy_bi_int_values() stores values of' ) print ( ' the integral of the Airy Bi function.' ) print ( '' ) print ( ' X FX' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = airy_bi_int_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def airy_bi_prime_values ( n_data ): #*****************************************************************************80 # ## 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: # # 16 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the derivative of the Airy BI function. # import numpy as np n_max = 11 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data =n_data + 1 return n_data, x, fx def airy_bi_prime_values_test ( ): #*****************************************************************************80 # ## airy_bi_prime_values_test() tests airy_bi_prime_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 16 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'airy_bi_prime_values_test():' ) print ( ' airy_bi_prime_values() stores values of' ) print ( ' the derivative of the Airy Bi function.' ) print ( '' ) print ( ' X FX' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = airy_bi_prime_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16f' % ( x, fx ) ) return def airy_bi_values ( n_data ): #*****************************************************************************80 # ## 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: # # 14 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the Airy BI function. # import numpy as np n_max = 11 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def airy_bi_values_test ( ): #*****************************************************************************80 # ## airy_bi_values_test() tests airy_bi_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 14 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'airy_bi_values_test():' ) print ( ' airy_bi_values() stores values of' ) print ( ' the Airy function Bi(x).' ) print ( '' ) print ( ' X Bi(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, ai = airy_bi_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16f' % ( x, ai ) ) return def airy_cai_values ( n_data ): #*****************************************************************************80 # ## airy_cai_values() returns some values of the Airy Ai(x) with 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: # # 13 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # complex X, the argument of the function. # # complex CAI, the value of the Airy AI function. # import numpy as np n_max = 10 cai_vec = np.array ( [ \ 0.1352924163128814 + 0.0000000000000000j, \ 0.1433824486882056 - 0.1092193342707378j, \ 0.2215404472324631 - 0.2588711788891803j, \ 0.4763929771766866 - 0.3036484220291284j, \ 0.5983692170633874 - 0.08154602160771214j, \ 0.5355608832923521 + 0.00000000000000000j, \ 0.5983692170633874 + 0.08154602160771214j, \ 0.4763929771766866 + 0.3036484220291284j, \ 0.2215404472324631 + 0.2588711788891803j, \ 0.1433824486882056 + 0.1092193342707378j ] ) x_vec = np.array ( [ \ 1.0000000000000000 + 0.0000000000000000j, \ 0.8090169943749474 + 0.5877852522924731j, \ 0.3090169943749474 + 0.9510565162951536j, \ -0.3090169943749474 + 0.9510565162951536j, \ -0.8090169943749474 + 0.5877852522924731j, \ -1.0000000000000000 + 0.0000000000000000j, \ -0.8090169943749474 - 0.5877852522924731j, \ -0.3090169943749474 - 0.9510565162951536j, \ 0.3090169943749474 - 0.9510565162951536j, \ 0.8090169943749474 - 0.5877852522924731j ] ) if ( n_data < 0 ): n_data = 0 if ( n_max <= n_data ): n_data = 0 x = 0.0 cai = 0.0 else: x = x_vec[n_data] cai = cai_vec[n_data] n_data = n_data + 1 return n_data, x, cai def airy_cai_values_test ( ): #*****************************************************************************80 # ## airy_cai_values_test() tests airy_cai_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 10 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'airy_cai_values_test():' ) print ( ' airy_cai_values() stores values of' ) print ( ' the complex Airy function Ai(x).' ) print ( '' ) print ( ' X.real X.imag Ai.real Ai.imag' ) print ( '' ) n_data = 0 while ( True ): n_data, x, cai = airy_cai_values ( n_data ) if ( n_data == 0 ): break print ( ' %12g %12g %24.16g %24.16g' \ % ( x.real, x.imag, cai.real, cai.imag ) ) return def airy_cbi_values ( n_data ): #*****************************************************************************80 # ## airy_cbi_values() returns some values of the Airy Bi(x) with 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: # # AiryBi[x] # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 13 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # complex X, the argument of the function. # # complex CBI, the value of the Airy BI function. # import numpy as np n_max = 10 cbi_vec = np.array ( ( \ 1.207423594952871 + 0.0000000000000000j, \ 0.9127160108293936 + 0.3800456133135556j, \ 0.6824453575635721 + 0.3343047153635002j, \ 0.5726265660086474 + 0.3988641086982559j, \ 0.2511841251049547 + 0.3401447690712719j, \ 0.1039973894969446 + 0.0000000000000000j, \ 0.2511841251049547 - 0.3401447690712719j, \ 0.5726265660086474 - 0.3988641086982559j, \ 0.6824453575635721 - 0.3343047153635002j, \ 0.9127160108293936 - 0.3800456133135556j ) ) x_vec = np.array ( ( \ 1.0000000000000000 + 0.0000000000000000j, \ 0.8090169943749474 + 0.5877852522924731j, \ 0.3090169943749474 + 0.9510565162951536j, \ -0.3090169943749474 + 0.9510565162951536j, \ -0.8090169943749474 + 0.5877852522924731j, \ -1.0000000000000000 + 0.0000000000000000j, \ -0.8090169943749474 - 0.5877852522924731j, \ -0.3090169943749474 - 0.9510565162951536j, \ 0.3090169943749474 - 0.9510565162951536j, \ 0.8090169943749474 - 0.5877852522924731j ) ) if ( n_data < 0 ): n_data = 0 if ( n_max <= n_data ): n_data = 0 x = 0.0 cbi = 0.0 else: x = x_vec[n_data] cbi = cbi_vec[n_data] n_data = n_data + 1 return n_data, x, cbi def airy_cbi_values_test ( ): #*****************************************************************************80 # ## airy_cbi_values_test() tests airy_cbi_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 10 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'airy_cbi_values_test():' ) print ( ' airy_cbi_values() stores values of' ) print ( ' the complex Airy function Bi(x).' ) print ( '' ) print ( ' X.real X.imag Bi.real(X) Bi.imag(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, cbi = airy_cbi_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %24.16f %24.16g' \ % ( x.real, x.imag, cbi.real, cbi.imag ) ) return def airy_gi_values ( n_data ): #*****************************************************************************80 # ## airy_gi_values() returns some values of the Airy Gi(x) function. # # Discussion: # # The function is defined by: # # airy_gi(x) = Integral ( 0 <= t < infinity ) sin ( x*t+t^3/3) dt / pi # # The data was reported by McLeod. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 14 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Allan McLeod, # Algorithm 757, MISCFUN: A software package to compute uncommon # special functions, # ACM transactions on Mathematical Software, # Volume 22, Number 3, September 1996, pages 288-301. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def airy_gi_values_test ( ): #*****************************************************************************80 # ## airy_gi_values_test() tests airy_gi_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'airy_gi_values_test():' ) print ( ' airy_gi_values() stores values of' ) print ( ' the Airy function Gi(x).' ) print ( '' ) print ( ' X Gi(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, ai = airy_gi_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, ai ) ) return def airy_hi_values ( n_data ): #*****************************************************************************80 # ## airy_hi_values() returns some values of the Airy Hi(x) function. # # Discussion: # # The function is defined by: # # airy_hi(x) = Integral ( 0 <= t < infinity ) exp ( x*t-t^3/3) dt / pi # # The data was reported by McLeod. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Allan McLeod, # Algorithm 757, MISCFUN: A software package to compute uncommon # special functions, # ACM transactions on Mathematical Software, # Volume 22, Number 3, September 1996, pages 288-301. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def airy_hi_values_test ( ): #*****************************************************************************80 # ## airy_hi_values_test() tests airy_hi_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'airy_hi_values_test():' ) print ( ' airy_hi_values() stores values of' ) print ( ' the Airy function Hi(x).' ) print ( '' ) print ( ' X Hi(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, ai = airy_hi_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, ai ) ) return def arccosh_values ( n_data ): #*****************************************************************************80 # ## 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: # # 21 December 2014 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 15 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def arccosh_values_test ( ): #*****************************************************************************80 # ## arccosh_values_test() tests arccosh_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'arccosh_values_test():' ) print ( ' arccosh_values() stores values of' ) print ( ' the hyperbolic arc cosine function.' ) print ( '' ) print ( ' X F(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = arccosh_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def arccos_values ( n_data ): #*****************************************************************************80 # ## 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: # # 22 December 2014 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 12 fx_vec = np.array ( ( \ 1.6709637479564564156, \ 1.5707963267948966192, \ 1.4706289056333368229, \ 1.3694384060045658278, \ 1.2661036727794991113, \ 1.1592794807274085998, \ 1.0471975511965977462, \ 0.92729521800161223243, \ 0.79539883018414355549, \ 0.64350110879328438680, \ 0.45102681179626243254, \ 0.00000000000000000000 ) ) x_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def arccos_values_test ( ): #*****************************************************************************80 # ## arccos_values_test() tests arccos_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'arccos_values_test():' ) print ( ' arccos_values() stores values of' ) print ( ' the arc cosine function.' ) print ( '' ) print ( ' X F(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = arccos_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def arcsinh_values ( n_data ): #*****************************************************************************80 # ## 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ -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 ) ) x_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def arcsinh_values_test ( ): #*****************************************************************************80 # ## arcsinh_values_test() tests arcsinh_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'arcsinh_values_test():' ) print ( ' arcsinh_values() stores values of' ) print ( ' the hyperbolic arc sine function.' ) print ( '' ) print ( ' X F(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = arcsinh_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def arcsin_values ( n_data ): #*****************************************************************************80 # ## 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: # # 22 December 2014 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 12 fx_vec = np.array ( ( \ -0.10016742116155979635, \ 0.00000000000000000000, \ 0.10016742116155979635, \ 0.20135792079033079146, \ 0.30469265401539750797, \ 0.41151684606748801938, \ 0.52359877559829887308, \ 0.64350110879328438680, \ 0.77539749661075306374, \ 0.92729521800161223243, \ 1.1197695149986341867, \ 1.5707963267948966192 ) ) x_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def arcsin_values_test ( ): #*****************************************************************************80 # ## arcsin_values_test() tests arcsin_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'arcsin_values_test():' ) print ( ' arcsin_values() stores values of' ) print ( ' the arc sine function.' ) print ( '' ) print ( ' X F(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = arcsin_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def arctan2_values ( n_data ): #*****************************************************************************80 # ## 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: # # 23 December 2014 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, Y, the arguments of the function. # # real F, the value of the function. # import numpy as np n_max = 19 f_vec = np.array ( ( \ -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 ) ) x_vec = np.array ( ( \ 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 ) ) y_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 x = 0.0 y = 0.0 f = 0.0 else: x = x_vec[n_data] y = y_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, x, y, f def arctan2_values_test ( ): #*****************************************************************************80 # ## arctan2_values_test() tests arctan2_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'arctan2_values_test():' ) print ( ' arctan2_values() stores values of' ) print ( ' the arc tangent function.' ) print ( '' ) print ( ' X Y F(X,Y)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, y, f = arctan2_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %24.16g' % ( x, y, f ) ) return def arctanh_values ( n_data ): #*****************************************************************************80 # ## 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 December 2014 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 15 fx_vec = np.array ( ( \ -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 ) ) x_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def arctanh_values_test ( ): #*****************************************************************************80 # ## arctanh_values_test() tests arctanh_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'arctanh_values_test():' ) print ( ' arctanh_values() stores values of' ) print ( ' the hyperbolic arc tangent function.' ) print ( '' ) print ( ' X F(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = arctanh_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def arctan_int_values ( n_data ): #*****************************************************************************80 # ## 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: # # 23 December 2014 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def arctan_int_values_test ( ): #*****************************************************************************80 # ## arctan_int_values_test() tests arctan_int_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'arctan_int_values_test():' ) print ( ' arctan_int_values() stores values of' ) print ( ' the arc tangent integral.' ) print ( '' ) print ( ' X F(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = arctan_int_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def arctan_values ( n_data ): #*****************************************************************************80 # ## 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: # # 23 December 2014 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 11 fx_vec = np.array ( ( \ 0.00000000000000000000, \ 0.24497866312686415417, \ 0.32175055439664219340, \ 0.46364760900080611621, \ 0.78539816339744830962, \ 1.1071487177940905030, \ 1.2490457723982544258, \ 1.3258176636680324651, \ 1.3734007669450158609, \ 1.4711276743037345919, \ 1.5208379310729538578 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def arctan_values_test ( ): #*****************************************************************************80 # ## arctan_values_test() tests arctan_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'arctan_values_test():' ) print ( ' arctan_values() stores values of' ) print ( ' the arc tangent function.' ) print ( '' ) print ( ' X F(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = arctan_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bei0_values ( n_data ): #*****************************************************************************80 # ## 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, # LC: QA47.A34, # ISBN: 0-486-61272-4. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Cambridge University Press, 1999, # LC: QA76.95.W65, # ISBN: 0-521-64314-7. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 11 fx_vec = np.array ( ( \ 0.0000000000000000, \ 0.06249321838219946, \ 0.2495660400366597, \ 0.5575600623030867, \ 0.9722916273066612, \ 1.457182044159804, \ 1.937586785266043, \ 2.283249966853915, \ 2.292690322699300, \ 1.686017203632139, \ 0.1160343815502004 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bei0_values_test ( ): #*****************************************************************************80 # ## bei0_values_test() tests bei0_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 24 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bei0_values_test():' ) print ( ' bei0_values() stores values of the Kelvin BEI function. of order 0.' ) print ( '' ) print ( ' X BEI(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bei0_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bei1_values ( n_data ): #*****************************************************************************80 # ## bei0_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: # # 25 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz, Irene Stegun, # Handbook of Mathematical Functions, # National Bureau of Standards, 1964, # LC: QA47.A34, # ISBN: 0-486-61272-4. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Cambridge University Press, 1999, # LC: QA76.95.W65, # ISBN: 0-521-64314-7. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 11 fx_vec = np.array ( ( \ 0.0000000000000000, \ 0.1711951797170153, \ 0.3075566313755366, \ 0.3678649890020899, \ 0.2997754370020335, \ 0.03866844396595048, \ -0.4874541770160708, \ -1.344042373111174, \ -2.563821688561078, \ -4.105685408400878, \ -5.797907901792625 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bei1_values_test ( ): #*****************************************************************************80 # ## bei1_values_test() tests bei1_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bei1_values_test():' ) print ( ' bei1_values() stores values of the Kelvin BEI function of order 1.' ) print ( '' ) print ( ' X BEI(1,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bei1_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bell_values ( n_data ): #*****************************************************************************80 # ## 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}] # # 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. # # Example: # # 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: # # 23 November 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer N, the order of the Bell number. # # integer C, the value of the Bell number. # import numpy as np n_max = 11 c_vec = np.array ( ( 1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975 ) ) n_vec = np.array ( ( 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ) ) if ( n_data < 0 ): n_data = 0 if ( n_max <= n_data ): n_data = 0 n = 0 c = 0 else: n = n_vec[n_data] c = c_vec[n_data] n_data = n_data + 1 return n_data, n, c def bell_values_test ( ): #*****************************************************************************80 # ## bell_values_test() tests bell_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 November 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bell_values_test():' ) print ( ' bell_values() returns values of' ) print ( ' the Bell numbers.' ) print ( '' ) print ( ' N BELL(N)' ) print ( '' ) n_data = 0 while ( True ): n_data, n, c = bell_values ( n_data ) if ( n_data == 0 ): break print ( '%6d %10d' % ( n, c ) ) return def ber0_values ( n_data ): #*****************************************************************************80 # ## 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, # LC: QA47.A34, # ISBN: 0-486-61272-4. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Cambridge University Press, 1999, # LC: QA76.95.W65, # ISBN: 0-521-64314-7. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 11 fx_vec = np.array ( ( \ 1.0000000000000000, \ 0.9990234639908383, \ 0.9843817812130869, \ 0.9210721835462558, \ 0.7517341827138082, \ 0.3999684171295313, \ -0.2213802495986939, \ -1.193598179589928, \ -2.563416557258580, \ -4.299086551599756, \ -6.230082478666358 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def ber0_values_test ( ): #*****************************************************************************80 # ## ber0_values_test() tests ber0_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'ber0_values_test():' ) print ( ' ber0_values() stores values of the Kelvin BER function. of order 0.' ) print ( '' ) print ( ' X BER(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = ber0_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def ber1_values ( n_data ): #*****************************************************************************80 # ## 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: # # 25 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz, Irene Stegun, # Handbook of Mathematical Functions, # National Bureau of Standards, 1964, # LC: QA47.A34, # ISBN: 0-486-61272-4. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Cambridge University Press, 1999, # LC: QA76.95.W65, # ISBN: 0-521-64314-7. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 11 fx_vec = np.array ( ( \ 0.0000000000000000, \ -0.1822431237551121, \ -0.3958682610197114, \ -0.6648654179597691, \ -0.9970776519264285, \ -1.373096897645111, \ -1.732644221128481, \ -1.959644131289749, \ -1.869248459031899, \ -1.202821631480086, \ 0.3597766667766728 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def ber1_values_test ( ): #*****************************************************************************80 # ## ber1_values_test() tests ber1_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 25 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'ber1_values_test():' ) print ( ' ber1_values() stores values of the Kelvin BER function of order 1.' ) print ( '' ) print ( ' X BER(1,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = ber1_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bernoulli_number_values ( n_data ): #*****************************************************************************80 # ## 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: # # 22 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer N, the order of the Bernoulli number. # # real C, the value of the Bernoulli number. # import numpy as np n_max = 10 c_vec = np.array ( ( 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 ) ) n_vec = np.array ( ( 0, 1, 2, 3, 4, 6, 8, 10, 20, 30 ) ) if ( n_data < 0 ): n_data = 0 if ( n_max <= n_data ): n_data = 0 n = 0 c = 0 else: n = n_vec[n_data] c = c_vec[n_data] n_data = n_data + 1 return n_data, n, c def bernoulli_number_values_test ( ): #*****************************************************************************80 # ## bernoulli_number_values_test() tests bernoulli_number_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bernoulli_number_values_test():' ) print ( ' bernoulli_number_values() returns values of' ) print ( ' the Bernoulli numbers.' ) print ( '' ) print ( ' N bernoulli_number(N)' ) print ( '' ) n_data = 0 while ( True ): n_data, n, c = bernoulli_number_values ( n_data ) if ( n_data == 0 ): break print ( '%6d %14.6g' % ( n, c ) ) return def bernoulli_poly_values ( n_data ): #*****************************************************************************80 # ## 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: # # 26 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer N, the order of the Bernoulli polynomial. # # real X, the argument of the Bernoulli polynomial. # # real FX, the value of the Bernoulli polynomial. # import numpy as np n_max = 27 fx_vec = np.array ( ( \ 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 ) ) n_vec = np.array ( ( \ 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 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 n = 0 x = 0.0 fx = 0.0 else: n = n_vec[n_data] x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, n, x, fx def bernoulli_poly_values_test ( ): #*****************************************************************************80 # ## bernoulli_poly_values_test() tests bernoulli_poly_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 26 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bernoulli_poly_values_test():' ) print ( ' bernoulli_poly_values() stores values of the Bernoulli polynomials.' ) print ( '' ) print ( ' N X FX' ) print ( '' ) n_data = 0 while ( True ): n_data, n, x, fx = bernoulli_poly_values ( n_data ) if ( n_data == 0 ): break print ( ' %6d %12f %24.16g' % ( n, x, fx ) ) return def bernstein_poly_01_values ( n_data ): #*****************************************************************************80 # ## 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: # # 26 December 2014 # # Author: # # John Burkardt # # Reference: # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer N, the degree of the polynomial. # # integer K, the index of the polynomial. # # real X, the argument of the polynomial. # # real F, the value of the polynomial B(N,K)(X). # import numpy as np n_max = 15 f_vec = np.array ( ( \ 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 ) ) k_vec = np.array ( ( \ 0, \ 0, 1, \ 0, 1, 2, \ 0, 1, 2, 3, \ 0, 1, 2, 3, 4 )) n_vec = np.array ( ( \ 0, \ 1, 1, \ 2, 2, 2, \ 3, 3, 3, 3, \ 4, 4, 4, 4, 4 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 n = 0 k = 0 x = 0.0 f = 0.0 else: n = n_vec[n_data] k = k_vec[n_data] x = x_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, n, k, x, f def bernstein_poly_01_values_test ( ): #*****************************************************************************80 # ## bernstein_poly_01_values_test() tests bernstein_poly_01_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 26 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bernstein_poly_01_values_test():' ) print ( ' bernstein_poly_01_values() stores values of Bernstein polynomials.' ) print ( '' ) print ( ' N K X F' ) print ( '' ) n_data = 0 while ( True ): n_data, n, k, x, f = bernstein_poly_01_values ( n_data ) if ( n_data == 0 ): break print ( ' %6d %6d %12f %24.16g' % ( n, k, x, f ) ) return def bessel_i0_int_values ( n_data ): #*****************************************************************************80 # ## 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: # # 30 December 2014 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_i0_int_values_test ( ): #*****************************************************************************80 # ## bessel_i0_int_values_test() tests bessel_i0_int_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_i0_int_values_test():' ) print ( ' bessel_i0_int_values() stores values of the Bessel I integral. of order 0.' ) print ( '' ) print ( ' X int_i(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_i0_int_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_i0_spherical_values ( n_data ): #*****************************************************************************80 # ## 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: # # 30 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz, Irene Stegun, # Handbook of Mathematical Functions, # National Bureau of Standards, 1964, # LC: QA47.A34, # ISBN: 0-486-61272-4. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 21 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_i0_spherical_values_test ( ): #*****************************************************************************80 # ## bessel_i0_spherical_values_test() tests bessel_i0_spherical_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_i0_spherical_values_test():' ) print ( ' bessel_i0_spherical_values: values of the spherical Bessel I function. of order 0.' ) print ( '' ) print ( ' X Spherical I(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_i0_spherical_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_i0_values ( n_data ): #*****************************************************************************80 # ## 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: # # 26 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_i0_values_test ( ): #*****************************************************************************80 # ## bessel_i0_values_test() tests bessel_i0_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 26 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_i0_values_test():' ) print ( ' bessel_i0_values() stores values of the Bessel I function. of order 0.' ) print ( '' ) print ( ' X I(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_i0_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_i1_spherical_values ( n_data ): #*****************************************************************************80 # ## 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: # # 31 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz, Irene Stegun, # Handbook of Mathematical Functions, # National Bureau of Standards, 1964, # LC: QA47.A34, # ISBN: 0-486-61272-4. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 21 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_i1_spherical_values_test ( ): #*****************************************************************************80 # ## bessel_i1_spherical_values_test() tests bessel_i1_spherical_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 31 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_i1_spherical_values_test():' ) print ( ' bessel_i1_spherical_values() stores values of the spherical Bessel I function. of order 1.' ) print ( '' ) print ( ' X Spherical I(1,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_i1_spherical_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_i1_values ( n_data ): #*****************************************************************************80 # ## 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: # # 27 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_i1_values_test ( ): #*****************************************************************************80 # ## bessel_i1_values_test() tests bessel_i1_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 27 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_i1_values_test():' ) print ( ' bessel_i1_values() stores values of the Bessel I function. of order 1.' ) print ( '' ) print ( ' X I(1,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_i1_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_in_values ( n_data ): #*****************************************************************************80 # ## 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: # # 02 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer NU, the order of the function. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 28 fx_vec = np.array ( ( \ 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 ) ) nu_vec = np.array ( ( \ 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 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 nu = 0 x = 0.0 fx = 0.0 else: nu = nu_vec[n_data] x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, nu, x, fx def bessel_in_values_test ( ): #*****************************************************************************80 # ## bessel_in_values_test() tests bessel_in_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 02 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_in_values_test():' ) print ( ' bessel_in_values() stores values of the Bessel I function. of order NU.' ) print ( '' ) print ( ' NU X I(NU,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, nu, x, fx = bessel_in_values ( n_data ) if ( n_data == 0 ): break print ( ' %4d %12f %24.16g' % ( nu, x, fx ) ) return def bessel_ix_values ( n_data ): #*****************************************************************************80 # ## 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: # # 06 January 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. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real NU, the order of the function. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 28 fx_vec = np.array ( ( \ 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 ) ) nu_vec = np.array ( ( \ 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 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 nu = 0 x = 0.0 fx = 0.0 else: nu = nu_vec[n_data] x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, nu, x, fx def bessel_ix_values_test ( ): #*****************************************************************************80 # ## bessel_ix_values_test() tests bessel_ix_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_ix_values_test():' ) print ( ' bessel_ix_values() stores values of the Bessel I function. of real order NU.' ) print ( '' ) print ( ' NU X I(NU,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, nu, x, fx = bessel_ix_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %24.16g' % ( nu, x, fx ) ) return def bessel_j0_int_values ( n_data ): #*****************************************************************************80 # ## 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: # # 11 January 2015 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_j0_int_values_test ( ): #*****************************************************************************80 # ## bessel_j0_int_values_test() tests bessel_j0_int_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_j0_int_values_test():' ) print ( ' bessel_j0_int_values() stores values of the Bessel J integral. of order 0.' ) print ( '' ) print ( ' X int_j(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_j0_int_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_j0_spherical_values ( n_data ): #*****************************************************************************80 # ## 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: # # 08 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 21 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_j0_spherical_values_test ( ): #*****************************************************************************80 # ## bessel_j0_spherical_values_test() tests bessel_j0_spherical_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 08 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_j0_spherical_values_test():' ) print ( ' bessel_j0_spherical_values() stores values of the spherical Bessel J function. of order 0.' ) print ( '' ) print ( ' X Spherical J(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_j0_spherical_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_j0_values ( n_data ): #*****************************************************************************80 # ## 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: # # 15 September 2004 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # n_max = 21; fx_vec = [ \ ]; x_vec = [ \ ]; import numpy as np n_max = 21 fx_vec = np.array ( ( \ -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 ) ) x_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_j0_values_test ( ): #*****************************************************************************80 # ## bessel_j0_values_test() tests bessel_j0_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 27 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_j0_values_test():' ) print ( ' bessel_j0_values() stores values of the Bessel J function. of order 0.' ) print ( '' ) print ( ' X J(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_j0_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_j0_zero_values ( n_data ): #*****************************************************************************80 # ## bessel_j0_values() returns some zeros of the J0 Bessel function. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 04 January 2016 # # Author: # # John Burkardt # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer K, the index of the zero. # # real FX, the value of the function. # import numpy as np n_max = 30 fx_vec = np.array ( [ \ 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 ] ) k_vec = np.array ( [ \ 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 n_data, k, fx def bessel_j0_zero_values_test ( ): #*****************************************************************************80 # ## bessel_j0_zero_values_test() tests bessel_j0_zero_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 04 January 2016 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_j0_zero_values_test():' ) print ( ' bessel_j0_zero_values stores zeros of' ) print ( ' the Bessel J0 function.' ) print ( '' ) print ( ' K X(K)' ) print ( '' ) n_data = 0 while ( True ): n_data, k, fx = bessel_j0_zero_values ( n_data ) if ( n_data == 0 ): break print ( ' %12d %24.16g' % ( k, fx ) ) return def bessel_j1_spherical_values ( n_data ): #*****************************************************************************80 # ## 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: # # 10 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 21 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_j1_spherical_values_test ( ): #*****************************************************************************80 # ## bessel_j1_spherical_values_test() tests bessel_j1_spherical_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 10 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_j1_spherical_values_test():' ) print ( ' bessel_j1_spherical_values() stores values of the spherical Bessel J function. of order 1.' ) print ( '' ) print ( ' X Spherical J(1,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_j1_spherical_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_j1_values ( n_data ): #*****************************************************************************80 # ## 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: # # 27 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 21 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_j1_values_test ( ): #*****************************************************************************80 # ## bessel_j1_values_test() tests bessel_j1_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 27 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_j1_values_test():' ) print ( ' bessel_j1_values() stores values of the Bessel J function. of order 1.' ) print ( '' ) print ( ' X J(1,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_j1_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_jn_values ( n_data ): #*****************************************************************************80 # ## 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: # # 10 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer NU, the order of the function. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) nu_vec = np.array ( ( \ 2, 2, 2, 2, \ 2, 5, 5, 5, \ 5, 5, 10, 10, \ 10, 10, 10, 20, \ 20, 20, 20, 20 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 nu = 0 x = 0.0 fx = 0.0 else: nu = nu_vec[n_data] x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, nu, x, fx def bessel_jn_values_test ( ): #*****************************************************************************80 # ## bessel_jn_values_test() tests bessel_jn_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 10 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_jn_values_test():' ) print ( ' bessel_jn_values() stores values of the Bessel J function. of order NU.' ) print ( '' ) print ( ' NU X J(NU,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, nu, x, fx = bessel_jn_values ( n_data ) if ( n_data == 0 ): break print ( ' %4d %12f %24.16g' % ( nu, x, fx ) ) return def bessel_j_spherical_values ( n_data ): #*****************************************************************************80 # ## 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 216 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz, Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer N, the index of the function. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 22 fx_vec = np.array ( [ \ 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 ] ) n_vec = np.array ( [ \ 0, \ 1, \ 2, \ 3, \ 4, \ 5, \ 6, \ 7, \ 8, \ 9, \ 10, \ 0, \ 1, \ 2, \ 3, \ 4, \ 5, \ 6, \ 7, \ 8, \ 9, \ 10 ] ) x_vec = np.array ( [ \ 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 if ( n_max <= n_data ): n_data = 0 n = 0 x = 0.0 fx = 0.0 else: n = n_vec[n_data] x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, n, x, fx def bessel_j_spherical_values_test ( ): #*****************************************************************************80 # ## bessel_j_spherical_values_test() tests bessel_j_spherical_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 16 January 2016 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_j_spherical_values_test():' ) print ( ' bessel_j_spherical_values() stores values of the spherical Bessel J function. of order N.' ) print ( '' ) print ( ' N X Spherical J(N,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, n, x, fx = bessel_j_spherical_values ( n_data ) if ( n_data == 0 ): break print ( ' %2d %12f %24.16g' % ( n, x, fx ) ) return def bessel_jx_values ( n_data ): #*****************************************************************************80 # ## 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: # # 11 January 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. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real NU, the order of the function. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 28 fx_vec = np.array ( ( \ 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 ) ) nu_vec = np.array ( ( \ 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 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 nu = 0 x = 0.0 fx = 0.0 else: nu = nu_vec[n_data] x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, nu, x, fx def bessel_jx_values_test ( ): #*****************************************************************************80 # ## bessel_jx_values_test() tests bessel_jx_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 11 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_jx_values_test():' ) print ( ' bessel_jx_values() stores values of the Bessel J function. of real order NU.' ) print ( '' ) print ( ' NU X J(NU,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, nu, x, fx = bessel_jx_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %24.16g' % ( nu, x, fx ) ) return def bessel_k0_int_values ( n_data ): #*****************************************************************************80 # ## 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: # # 11 January 2015 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_k0_int_values_test ( ): #*****************************************************************************80 # ## bessel_k0_int_values_test() tests bessel_k0_int_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 11 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_k0_int_values_test():' ) print ( ' bessel_k0_int_values() stores values of the Bessel K integral. of order 0.' ) print ( '' ) print ( ' X int_k(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_k0_int_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_k0_values ( n_data ): #*****************************************************************************80 # ## 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: # # 28 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_k0_values_test ( ): #*****************************************************************************80 # ## bessel_k0_values_test() tests bessel_k0_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_k0_values_test():' ) print ( ' bessel_k0_values() stores values of the Bessel K function. of order 0.' ) print ( '' ) print ( ' X K(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_k0_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_k1_values ( n_data ): #*****************************************************************************80 # ## 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: # # 28 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_k1_values_test ( ): #*****************************************************************************80 # ## bessel_k1_values_test() tests bessel_k1_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_k1_values_test():' ) print ( ' bessel_k1_values() stores values of the Bessel K function. of order 1.' ) print ( '' ) print ( ' X K(1,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_k1_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_kn_values ( n_data ): #*****************************************************************************80 # ## 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: # # 11 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer NU, the order of the function. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 28 fx_vec = np.array ( ( \ 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 ) ) nu_vec = np.array ( ( \ 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 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 nu = 0 x = 0.0 fx = 0.0 else: nu = nu_vec[n_data] x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, nu, x, fx def bessel_kn_values_test ( ): #*****************************************************************************80 # ## bessel_kn_values_test() tests bessel_kn_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 11 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_kn_values_test():' ) print ( ' bessel_kn_values() stores values of the Bessel K function. of order NU.' ) print ( '' ) print ( ' NU X K(NU,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, nu, x, fx = bessel_kn_values ( n_data ) if ( n_data == 0 ): break print ( ' %4d %12f %24.16g' % ( nu, x, fx ) ) return def bessel_kx_values ( n_data ): #*****************************************************************************80 # ## 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: # # 12 January 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. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real NU, the order of the function. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 28 fx_vec = np.array ( ( \ 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 ) ) nu_vec = np.array ( ( \ 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 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 nu = 0 x = 0.0 fx = 0.0 else: nu = nu_vec[n_data] x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, nu, x, fx def bessel_kx_values_test ( ): #*****************************************************************************80 # ## bessel_kx_values_test() tests bessel_kx_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_kx_values_test():' ) print ( ' bessel_kx_values() stores values of the Bessel K function. of real order NU.' ) print ( '' ) print ( ' NU X K(NU,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, nu, x, fx = bessel_kx_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %24.16g' % ( nu, x, fx ) ) return def bessel_y0_int_values ( n_data ): #*****************************************************************************80 # ## 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: # # 12 January 2015 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ -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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_y0_int_values_test ( ): #*****************************************************************************80 # ## bessel_y0_int_values_test() tests bessel_y0_int_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_y0_int_values_test():' ) print ( ' bessel_y0_int_values() stores values of the Bessel Y integral. of order 0.' ) print ( '' ) print ( ' X int_y(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_y0_int_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_y0_spherical_values ( n_data ): #*****************************************************************************80 # ## 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: # # 13 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz, Irene Stegun, # Handbook of Mathematical Functions, # National Bureau of Standards, 1964, # LC: QA47.A34, # ISBN: 0-486-61272-4. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 21 fx_vec = np.array ( ( \ -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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_y0_spherical_values_test ( ): #*****************************************************************************80 # ## bessel_y0_spherical_values_test() tests bessel_y0_spherical_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 13 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_y0_spherical_values_test():' ) print ( ' bessel_y0_spherical_values() stores values of the spherical Bessel Y function. of order 0.' ) print ( '' ) print ( ' X Spherical Y(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_y0_spherical_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_y0_values ( n_data ): #*****************************************************************************80 # ## 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: # # 29 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 16 fx_vec = np.array ( ( \ -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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_y0_values_test ( ): #*****************************************************************************80 # ## bessel_y0_values_test() tests bessel_y0_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 29 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_y0_values_test():' ) print ( ' bessel_y0_values() stores values of the Bessel Y function. of order 0.' ) print ( '' ) print ( ' X Y(0,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_y0_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_y0_zero_values ( n_data ): #*****************************************************************************80 # ## bessel_y0_values() returns some zeros of the Y0 Bessel function. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 March 2024 # # Author: # # John Burkardt # # Input: # # integer n_data: the user sets n_data to 0 before the first call. # # Output: # # integer n_data: on each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer k: the index of the zero. # # real fx: the value of the function. # import numpy as np n_max = 32 fx_vec = np.array ( [ \ 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 ] ) k_vec = np.array ( [ \ 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 n_data, k, fx def bessel_y0_zero_values_test ( ): #*****************************************************************************80 # ## bessel_y0_zero_values_test() tests bessel_y0_zero_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 March 2024 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_y0_zero_values_test():' ) print ( ' bessel_y0_zero_values stores zeros of' ) print ( ' the Bessel Y0 function.' ) print ( '' ) print ( ' K X(K)' ) print ( '' ) n_data = 0 while ( True ): n_data, k, fx = bessel_y0_zero_values ( n_data ) if ( n_data == 0 ): break print ( ' %12d %24.16g' % ( k, fx ) ) return def bessel_y1_spherical_values ( n_data ): #*****************************************************************************80 # ## 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: # # 13 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz, Irene Stegun, # Handbook of Mathematical Functions, # National Bureau of Standards, 1964, # LC: QA47.A34, # ISBN: 0-486-61272-4. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 21 fx_vec = np.array ( ( \ -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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_y1_spherical_values_test ( ): #*****************************************************************************80 # ## bessel_y1_spherical_values_test() tests bessel_y1_spherical_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 13 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_y1_spherical_values_test():' ) print ( ' bessel_y1_spherical_values() stores values of the spherical Bessel Y function of order 1.' ) print ( '' ) print ( ' X Spherical Y(1,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_y1_spherical_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_y1_values ( n_data ): #*****************************************************************************80 # ## 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: # # 29 December 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 16 fx_vec = np.array ( ( \ -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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def bessel_y1_values_test ( ): #*****************************************************************************80 # ## bessel_y1_values_test() tests bessel_y1_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 29 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_y1_values_test():' ) print ( ' bessel_y1_values() stores values of the Bessel Y function. of order 1.' ) print ( '' ) print ( ' X Y(1,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = bessel_y1_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16g' % ( x, fx ) ) return def bessel_yn_values ( n_data ): #*****************************************************************************80 # ## bessel_yn_values() returns some values of the Kn Bessel function. # # Discussion: # # In Mathematica, the function can be evaluated by: # # BesselY[n,x] # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 14 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer NU, the order of the function. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 20 fx_vec = np.array ( ( \ -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 ) ) nu_vec = np.array ( ( \ 2, 2, 2, 2, \ 2, 5, 5, 5, \ 5, 5, 10, 10, \ 10, 10, 10, 20, \ 20, 20, 20, 20 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 nu = 0 x = 0.0 fx = 0.0 else: nu = nu_vec[n_data] x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, nu, x, fx def bessel_yn_values_test ( ): #*****************************************************************************80 # ## bessel_yn_values_test() tests bessel_yn_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 14 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_yn_values_test():' ) print ( ' bessel_yn_values() stores values of the Bessel Y function. of order NU.' ) print ( '' ) print ( ' NU X Y(NU,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, nu, x, fx = bessel_yn_values ( n_data ) if ( n_data == 0 ): break print ( ' %4d %12f %24.16g' % ( nu, x, fx ) ) return def bessel_yx_values ( n_data ): #*****************************************************************************80 # ## 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 Kn 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: # # 14 January 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. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real NU, the order of the function. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 28 fx_vec = np.array ( ( \ -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 ) ) nu_vec = np.array ( ( \ 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 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 nu = 0 x = 0.0 fx = 0.0 else: nu = nu_vec[n_data] x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, nu, x, fx def bessel_yx_values_test ( ): #*****************************************************************************80 # ## bessel_yx_values_test() tests bessel_yx_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 14 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'bessel_yx_values_test():' ) print ( ' bessel_yx_values() stores values of the Bessel Y function. of real order NU.' ) print ( '' ) print ( ' NU X Y(NU,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, nu, x, fx = bessel_yx_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %24.16g' % ( nu, x, fx ) ) return def beta_cdf_values ( n_data ): #*****************************************************************************80 # ## 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 to X) T^(A-1) * (1-T)^(B-1) dT # / Integral (0 to 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: # # 15 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Karl Pearson, # Tables of the Incomplete Beta Function, # Cambridge University Press, 1968. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real A, B, the parameters of the function. # # real X, the argument of the function. # # real F, the value of the function. # import numpy as np n_max = 45 a_vec = np.array ( ( \ 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 )) b_vec = np.array ( ( \ 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 )) f_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 a = 0.0 b = 0.0 x = 0.0 f = 0.0 else: a = a_vec[n_data] b = b_vec[n_data] x = x_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, a, b, x, f def beta_cdf_values_test ( ): #*****************************************************************************80 # ## beta_cdf_values_test() tests beta_cdf_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 15 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'beta_cdf_values_test():' ) print ( ' beta_cdf_values() stores values of the BETA function.' ) print ( '' ) print ( ' A B X beta_cdf(A,B,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, a, b, x, f = beta_cdf_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %12f %24.16g' % ( a, b, x, f ) ) return def beta_inc_values ( n_data ): #*****************************************************************************80 # ## 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 to X) T^(A-1) * (1-T)^(B-1) dT # / Integral (0 to 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: # # 15 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Karl Pearson, # Tables of the Incomplete Beta Function, # Cambridge University Press, 1968. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real A, B, the parameters of the function. # # real X, the argument of the function. # # real F, the value of the function. # import numpy as np n_max = 45 a_vec = np.array ( ( \ 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 )) b_vec = np.array ( ( \ 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 )) f_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 a = 0.0 b = 0.0 x = 0.0 f = 0.0 else: a = a_vec[n_data] b = b_vec[n_data] x = x_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, a, b, x, f def beta_inc_values_test ( ): #*****************************************************************************80 # ## beta_inc_values_test() tests beta_inc_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 15 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'beta_inc_values_test():' ) print ( ' beta_inc_values() stores values of the BETA function.' ) print ( '' ) print ( ' A B X beta_inc(A,B,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, a, b, x, f = beta_inc_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %12f %24.16g' % ( a, b, x, f ) ) return def beta_log_values ( n_data ): #*****************************************************************************80 # ## 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: # # 15 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, Y, the arguments of the function. # # real FXY, the value of the function. # import numpy as np n_max = 17 f_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 ) ) y_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 y = 0.0 f = 0.0 else: x = x_vec[n_data] y = y_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, x, y, f def beta_log_values_test ( ): #*****************************************************************************80 # ## beta_log_values_test() tests beta_log_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 15 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'beta_log_values_test():' ) print ( ' beta_log_values() stores values of the Log(BETA) function.' ) print ( '' ) print ( ' X Y Log(BETA(X,Y))' ) print ( '' ) n_data = 0 while ( True ): n_data, x, y, f = beta_log_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %24.16g' % ( x, y, f ) ) return def beta_noncentral_cdf_values ( n_data ): #*****************************************************************************80 # ## 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: # # 15 January 2015 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # rea A, B, the shape parameters. # # real LAMDA, the noncentrality parameter. # It is necessary to misspell LAMBDA since Python uses it as a keyword. # # real X, the argument of the function. # # real FX, the value of the function. # import numpy as np n_max = 25 a_vec = np.array ( ( \ 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 )) b_vec = np.array ( ( \ 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 )) f_vec = np.array ( ( \ 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 ) ) lamda_vec = np.array ( ( \ 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 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 a = 0.0 b = 0.0 lamda = 0.0 x = 0.0 f = 0.0 else: a = a_vec[n_data] b = b_vec[n_data] lamda = lamda_vec[n_data] x = x_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, a, b, lamda, x, f def beta_noncentral_cdf_values_test ( ): #*****************************************************************************80 # ## beta_noncentral_cdf_values_test() tests beta_noncentral_cdf_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 15 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'beta_noncentral_cdf_values_test():' ) print ( ' beta_noncentral_cdf_values() stores values of the noncentral BETA CDF.' ) print ( '' ) print ( ' A B LAMDA X beta_noncentral_cdf(A,B,LAMDA,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, a, b, lamda, x, f = beta_noncentral_cdf_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %12f %12f %24.16g' % ( a, b, lamda, x, f ) ) return def beta_pdf_values ( n_data ): #*****************************************************************************80 # ## 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real ALPHA, BETA, the parameters of the function. # # real X, the argument of the function. # # real F, the value of the function. # import numpy as np n_max = 10 alpha_vec = np.array ( ( \ 1.092091484911879, \ 2.808477213834471, \ 1.287888961910225, \ 3.169828561512062, \ 2.006531407488083, \ 0.009191855792026001, \ 0.472723751058401, \ 4.204237253278341, \ 1.301514988836825, \ 1.758143299519481 )) beta_vec = np.array ( ( \ 4.781587882544648, \ 2.076535407379806, \ 0.549783967662353, \ 0.3086361453280091, \ 3.773367432107051, \ 4.487520304498656, \ 0.06808445791730976, \ 0.6155195788227712, \ 4.562418534907164, \ 4.114436583429598 )) f_vec = np.array ( ( \ 0.002826137156803199, \ 0.04208950342768649, \ 0.2184064957817208, \ 0.1335142301445414, \ 0.1070571849830009, \ 0.005796394377470491, \ 0.5518796772414584, \ 0.0, \ 2.87907465409348, \ 2.126992854611924 ) ) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 alpha = 0.0 beta = 0.0 x = 0.0 f = 0.0 else: alpha = alpha_vec[n_data] beta = beta_vec[n_data] x = x_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, alpha, beta, x, f def beta_pdf_values_test ( ): #*****************************************************************************80 # ## beta_pdf_values_test() tests beta_pdf_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 30 July 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'beta_pdf_values_test():' ) print ( ' beta_pdf_values() stores values of the BETA function.' ) print ( '' ) print ( ' ALPHA BETA X PDF()' ) print ( '' ) n_data = 0 while ( True ): n_data, alpha, beta, x, f = beta_pdf_values ( n_data ) if ( n_data == 0 ): break print ( ' %12g %12g %12g %24.16g' % ( alpha, beta, x, f ) ) return def beta_values ( n_data ): #*****************************************************************************80 # ## 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: # # 01 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real X, Y, the arguments of the function. # # real F, the value of the function. # import numpy as np n_max = 17 f_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 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 ) ) y_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 x = 0.0 y = 0.0 f = 0.0 else: x = x_vec[n_data] y = y_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, x, y, f def beta_values_test ( ): #*****************************************************************************80 # ## beta_values_test() tests beta_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 24 December 2014 # # Author: # # John Burkardt # print ( '' ) print ( 'beta_values_test():' ) print ( ' beta_values() stores values of the BETA function.' ) print ( '' ) print ( ' X Y BETA(X,Y)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, y, f = beta_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %24.16g' % ( x, y, f ) ) return def binomial_cdf_values ( n_data ): #*****************************************************************************80 # ## 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: # # 17 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Daniel Zwillinger, # CRC Standard Mathematical Tables and Formulae, # 30th Edition, CRC Press, 1996, pages 651-652. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer A, a parameter of the function. # # real B, a parameter of the function. # # integer X, the argument of the function. # # real F, the value of the function. # import numpy as np n_max = 17 a_vec = np.array ( ( \ 2, 2, 2, 2, \ 2, 4, 4, 4, \ 4, 10, 10, 10, \ 10, 10, 10, 10, \ 10 )) b_vec = np.array ( ( \ 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 )) f_vec = np.array ( ( \ 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 ) ) x_vec = np.array ( ( \ 0, 1, 2, 0, \ 1, 0, 1, 2, \ 3, 4, 4, 4, \ 4, 4, 4, 4, \ 4 ) ) if ( n_data < 0 ): n_data = 0 if ( n_max <= n_data ): n_data = 0 a = 0.0 b = 0.0 x = 0.0 f = 0.0 else: a = a_vec[n_data] b = b_vec[n_data] x = x_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, a, b, x, f def binomial_cdf_values_test ( ): #*****************************************************************************80 # ## binomial_cdf_values_test() tests binomial_cdf_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 17 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'binomial_cdf_values_test():' ) print ( ' binomial_cdf_values() stores values of the BINOMIAL CDF.' ) print ( '' ) print ( ' A B X binomial_cdf(A,B,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, a, b, x, f = binomial_cdf_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %12f %12d %24.16g' % ( a, b, x, f ) ) return def binomial_pdf_values ( n_data ): #*****************************************************************************80 # ## 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 and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Daniel Zwillinger, # CRC Standard Mathematical Tables and Formulae, # 30th Edition, CRC Press, 1996, pages 651-652. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer A, a parameter of the function. # # real B, a parameter of the function. # # integer X, the argument of the function. # # real F, the value of the function. # import numpy as np n_max = 10 a_vec = np.array ( ( \ 5, 12, 6, 13, 9, \ 1, 2, 17, 6, 8 )) b_vec = np.array ( ( \ 0.8295092339327006, \ 0.06611873491603133, \ 0.0438289977791071, \ 0.4495389603071763, \ 0.7972869541062562, \ 0.3507523379805466, \ 0.8590968552798568, \ 0.007512364073964213, \ 0.1136640464424993, \ 0.2671322702601793 )) f_vec = np.array ( ( \ 0.3927408939646697, \ 0.0006199968732461383, \ 0.764211224733124, \ 0.0004260353334364943, \ 0.302948289145794, \ 0.3507523379805466, \ 0.01985369619202562, \ 0.006854388879646552, \ 0.000002156446446382985, \ 0.0005691150511772053 ) ) x_vec = np.array ( ( \ 5, 5, 0, 0, 7, \ 1, 0, 2, 6, 7 ) ) if ( n_data < 0 ): n_data = 0 if ( n_max <= n_data ): n_data = 0 a = 0.0 b = 0.0 x = 0.0 f = 0.0 else: a = a_vec[n_data] b = b_vec[n_data] x = x_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, a, b, x, f def binomial_pdf_values_test ( ): #*****************************************************************************80 # ## binomial_pdf_values_test() tests binomial_pdf_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 27 July 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'binomial_pdf_values_test():' ) print ( ' binomial_pdf_values() stores values of the BINOMIAL PDF.' ) print ( '' ) print ( ' A B X binomial_pdf(A,B,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, a, b, x, f = binomial_pdf_values ( n_data ) if ( n_data == 0 ): break print ( ' %2d %24.16g %2d %24.16g' % ( a, b, x, f ) ) return def binomial_values ( n_data ): #*****************************************************************************80 # ## 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: # # 17 January 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer A, B, the arguments of the function. # # integer F, the value of the function. # import numpy as np n_max = 20 a_vec = np.array ( ( \ 1, 6, 6, 6, 15, \ 15, 15, 15, 15, 15, \ 15, 25, 25, 25, 25, \ 25, 25, 25, 25, 25 ) ) b_vec = np.array ( ( \ 0, 1, 3, 5, 1, \ 3, 5, 7, 9, 11, \ 13, 1, 3, 5, 7, \ 9, 11, 13, 15, 17 ) ) f_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 a = 0 b = 0 f = 0 else: a = a_vec[n_data] b = b_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, a, b, f def binomial_values_test ( ): #*****************************************************************************80 # ## binomial_values_test() tests binomial_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 17 January 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'binomial_values_test():' ) print ( ' binomial_values() stores values of the BINOMIAL function.' ) print ( '' ) print ( ' A B BINOMIAL(A,B)' ) print ( '' ) n_data = 0 while ( True ): n_data, a, b, f = binomial_values ( n_data ) if ( n_data == 0 ): break print ( ' %12d %12d %12d' % ( a, b, f ) ) return def bivariate_normal_cdf_values ( n_data ): #*****************************************************************************80 # ## 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.0;625 # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 February 2015 # # Author: # # John Burkardt # # Reference: # # J H Halton, # On the efficiency of certain quasi-random sequences of points # in evaluating multi-dimensional integrals, # Numerische Mathematik, # Volume 2, pages 84-90, 1960. # # J M 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer BASE, the base of the sequence. # # integer SEED, the index of the element of the sequence. # # real VALUE, the value of the SEED-th element of the # van der Corput sequence in base BASE. # import numpy as np n_max = 75 base_vec = np.array ( ( \ 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 )) seed_vec = np.array ( ( \ 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 )) value_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 base = 0 seed = 0 value = 0.0 else: base = base_vec[n_data] seed = seed_vec[n_data] value = value_vec[n_data] n_data = n_data + 1 return n_data, base, seed, value def van_der_corput_values_test ( ): #*****************************************************************************80 # ## van_der_corput_values_test() tests van_der_corput_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 February 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'van_der_corput_values_test():' ) print ( ' van_der_corput_values() stores values of the van_der_corput function.' ) print ( '' ) print ( ' BASE SEED VALUE' ) print ( '' ) n_data = 0 while ( True ): n_data, base, seed, value = van_der_corput_values ( n_data ) if ( n_data == 0 ): break print ( ' %12d %12d %24.12g' % ( base, seed, value ) ) return def viscosity_values ( n_data ): #*****************************************************************************80 # ## viscosity_values() returns some values of the viscosity function. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 February 2015 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real TC, the temperature, in degrees Celsius. # # real P, the pressure, in bar. # # real ETA, the viscosity, in MegaPascal seconds. # import numpy as np n_max = 34 eta_vec = np.array ( ( \ 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 )) p_vec = np.array ( ( \ 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 )) tc_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 tc = 0.0 p = 0.0 eta = 0.0 else: tc = tc_vec[n_data] p = p_vec[n_data] eta = eta_vec[n_data] n_data = n_data + 1 return n_data, tc, p, eta def viscosity_values_test ( ): #*****************************************************************************80 # ## viscosity_values_test() tests viscosity_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 February 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'viscosity_values_test():' ) print ( ' viscosity_values() stores values of the VISCOSITY function.' ) print ( '' ) print ( ' TC P ETA(TC,P)' ) print ( '' ) n_data = 0 while ( True ): n_data, tc, p, eta = viscosity_values ( n_data ) if ( n_data == 0 ): break print ( ' %12g %12g %24.16g' % ( tc, p, eta ) ) return def von_mises_cdf_values ( n_data ): #*****************************************************************************80 # ## von_mises_cdf_values() returns some values of the von Mises CDF. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 February 2015 # # Author: # # John Burkardt # # Reference: # # Kanti Mardia and Peter Jupp, # Directional Statistics, # Wiley, 2000, QA276.M335 # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real A, B, the parameters of the function. # # real X, the argument of the function. # # real F, the value of the function. # import numpy as np n_max = 23 a_vec = np.array ( ( \ 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 )) b_vec = np.array ( ( \ 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 )) f_vec = np.array ( ( \ 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 )) x_vec = np.array ( ( \ -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 if ( n_max <= n_data ): n_data = 0 a = 0.0 b = 0.0 x = 0.0 f = 0.0 else: a = a_vec[n_data] b = b_vec[n_data] x = x_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, a, b, x, f def von_mises_cdf_values_test ( ): #*****************************************************************************80 # ## von_mises_cdf_values_test() tests von_mises_cdf_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 February 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'von_mises_cdf_values_test():' ) print ( ' von_mises_cdf_values() stores values of the von Mises CDF.' ) print ( '' ) print ( ' A B X CDF(A,B,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, a, b, x, f = von_mises_cdf_values ( n_data ) if ( n_data == 0 ): break print ( ' %12g %12g %12g %24.16g' % ( a, b, x, f ) ) return def weekday_values ( n_data ): #*****************************************************************************80 # ## 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: # # 22 February 2015 # # 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer Y, M, D, the Common Era date. # # integer W, the day of the week. Sunday = 1. # import numpy as np n_max = 34 d_vec = np.array ( ( \ 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 )) m_vec = np.array ( ( \ 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 )) w_vec = np.array ( ( \ 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 )) y_vec = np.array ( ( \ - 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 if ( n_max <= n_data ): n_data = 0 y = 0 m = 0 d = 0 w = 0 else: y = y_vec[n_data] m = m_vec[n_data] d = d_vec[n_data] w = w_vec[n_data] n_data = n_data + 1 return n_data, y, m, d, w def weekday_values_test ( ): #*****************************************************************************80 # ## weekday_values_test() tests weekday_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 February 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'weekday_values_test():' ) print ( ' weekday_values stores() values of' ) print ( ' the weekday for a given Y/M/D date' ) print ( '' ) print ( ' Y M D W' ) print ( '' ) n_data = 0 while ( True ): n_data, y, m, d, w = weekday_values ( n_data ) if ( n_data == 0 ): break print ( ' %6d %6d %6d %6d' % ( y, m, d, w ) ) return def weibull_cdf_values ( n_data ): #*****************************************************************************80 # ## 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: # # 22 February 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real ALPHA, the first parameter of the distribution. # # real BETA, the second parameter of the distribution. # # real X, the argument of the function. # # real F, the value of the function. # import numpy as np n_max = 12 alpha_vec = np.array ( ( \ 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 )) beta_vec = np.array ( ( \ 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 )) f_vec = np.array ( ( \ 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 )) x_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 alpha = 0.0 beta = 0.0 x = 0.0 f = 0.0 else: alpha = alpha_vec[n_data] beta = beta_vec[n_data] x = x_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, alpha, beta, x, f def weibull_cdf_values_test ( ): #*****************************************************************************80 # ## weibull_cdf_values_test() tests weibull_cdf_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 February 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'weibull_cdf_values_test():' ) print ( ' weibull_cdf_values() stores values of the von Mises CDF.' ) print ( '' ) print ( ' ALPHA BETA X CDF(ALPHA,BETA,X)' ) print ( '' ) n_data = 0 while ( True ): n_data, alpha, beta, x, f = weibull_cdf_values ( n_data ) if ( n_data == 0 ): break print ( ' %12g %12g %12g %24.16g' % ( alpha, beta, x, f ) ) return def wright_omega_values ( n_data ): #*****************************************************************************80 # ## wright_omega_values() returns values of the Wright Omega function. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 15 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: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # complex Z, the argument of the function. # # complex FZ, the value of the function. # import numpy as np n_max = 10 fz_vec = np.array ( [ \ +0.5671432904097838 +0.0000000000000000j, \ +1.000000000000000 +0.0000000000000000j, \ +2.718281828459045 +0.0000000000000000j, \ -1.000000000000000 +0.0000000000000000j, \ -1.000000000000000 +0.0000000000000000j, \ -2.000000000000000 +0.0000000000000000j, \ -0.40637573995996 +0.0000000000000000j, \ +0.000000000000000 +1.0000000000000000j, \ -0.3181315052047641 +1.337235701430689j, \ +0.9372082083733697 +0.5054213160131512j ] ) z_vec = np.array ( [ \ +0.000000000000000 +0.000000000000000j, \ +1.000000000000000 +0.000000000000000j, \ +3.718281828459045 +0.000000000000000j, \ -1.000000000000000 +3.141592653589793j, \ -1.000000000000000 -3.141592653589793j, \ -1.306852819440055 +3.141592653589793j, \ -1.306852819440055 -3.141592653589793j, \ +0.000000000000000 +2.570796326794897j, \ +0.000000000000000 +3.141592653589793j, \ +1.000000000000000 +1.000000000000000j ] ) if ( n_data < 0 ): n_data = 0 if ( n_max <= n_data ): n_data = 0 z = 0.0 + 0.0j fz = 0.0 + 0.0j else: z = z_vec[n_data] fz = fz_vec[n_data] n_data = n_data + 1 return n_data, z, fz def wright_omega_values_test ( ): #*****************************************************************************80 # ## wright_omega_values_test() tests wright_omega_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 15 May 2016 # # Author: # # John Burkardt # print ( '' ) print ( 'wright_omega_values_test():' ) print ( ' wright_omega_values() stores values of' ) print ( ' the Wright Omega function.' ) print ( '' ) print ( ' Z.real Z.imag FZ.real FZ.imag' ) print ( '' ) n_data = 0 while ( True ): n_data, z, fz = wright_omega_values ( n_data ) if ( n_data == 0 ): break print ( ' %12g %12g %24.16g %24.16g' \ % ( z.real, z.imag, fz.real, fz.imag ) ) return def zeta_values ( n_data ): #*****************************************************************************80 # ## zeta_values() returns some values of the Riemann Zeta function. # # Discussion: # # ZETA(N) = sum ( 1 <= I < Infinity ) 1 / I^N # # In Mathematica, the function can be evaluated by: # # Zeta[n] # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 February 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # integer N, the argument of the Zeta function. # # real F, the value of the Zeta function. # import numpy as np n_max = 15 n_vec = np.array ( ( \ 2, \ 3, \ 4, \ 5, \ 6, \ 7, \ 8, \ 9, \ 10, \ 11, \ 12, \ 16, \ 20, \ 30, \ 40 )) f_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 n = 0 f = 0.0 else: n = n_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, n, f def zeta_values_test ( ): #*****************************************************************************80 # ## zeta_values_test() tests zeta_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 February 2015 # # Author: # # John Burkardt # print ( '' ) print ( 'zeta_values_test():' ) print ( ' zeta_values() stores values of the ZETA function.' ) print ( '' ) print ( ' N ZETA(N)' ) print ( '' ) n_data = 0 while ( True ): n_data, n, f = zeta_values ( n_data ) if ( n_data == 0 ): break print ( ' %6d %24.16f' % ( n, f ) ) return def zeta_m1_values ( n_data ): #*****************************************************************************80 # ## zeta_m1_values() returns some values of the Riemann Zeta function minus 1. # # Discussion: # # zeta_m1(N) = ZETA(N) - 1 # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 15 January 2016 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Input: # # integer n_data. The user sets n_data to 0 before the first call. # # Output: # # integer n_data. On each call, the routine increments n_data by 1, and # returns the corresponding data; when there is no more data, the # output value of n_data will be 0 again. # # real P, the argument. # # real F, the value. # import numpy as np n_max = 17 p_vec = np.array ( ( \ 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 )) f_vec = np.array ( ( \ 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 if ( n_max <= n_data ): n_data = 0 p = 0.0 f = 0.0 else: p = p_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, p, f def zeta_m1_values_test ( ): #*****************************************************************************80 # ## zeta_m1_values_test() tests zeta_m1_values(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 15 January 2017 # # Author: # # John Burkardt # print ( '' ) print ( 'zeta_m1_values_test():' ) print ( ' zeta_m1_values() stores values of the zeta_MINUS_ONE function.' ) print ( '' ) print ( ' N zeta_m1(N)' ) print ( '' ) n_data = 0 while ( True ): n_data, p, f = zeta_m1_values ( n_data ) if ( n_data == 0 ): break print ( ' %8f %24.16e' % ( p, f ) ) return if ( __name__ == '__main__' ): timestamp ( ) test_values_test ( ) timestamp ( )