function abram0 ( xvalue ) !*****************************************************************************80 ! !! ABRAM0 evaluates the Abramowitz function of order 0. ! ! Discussion: ! ! The function is defined by: ! ! ABRAM0(x) = Integral ( 0 <= t < infinity ) exp ( -t^2 - x / t ) dt ! ! The code uses Chebyshev expansions with the coefficients ! given to an accuracy of 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 07 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) ABRAM0, the value of the function. ! implicit none real ( kind = 8 ), dimension(0:8) :: ab0f = (/ & -0.68121927093549469816d0, & -0.78867919816149252495d0, & 0.5121581776818819543d-1, & -0.71092352894541296d-3, & 0.368681808504287d-5, & -0.917832337237d-8, & 0.1270202563d-10, & -0.1076888d-13, & 0.599d-17 /) real ( kind = 8 ), dimension(0:8) :: ab0g = (/ & -0.60506039430868273190d0, & -0.41950398163201779803d0, & 0.1703265125190370333d-1, & -0.16938917842491397d-3, & 0.67638089519710d-6, & -0.135723636255d-8, & 0.156297065d-11, & -0.112887d-14, & 0.55d-18 /) real ( kind = 8 ), dimension(0:8) :: ab0h = (/ & 1.38202655230574989705d0, & -0.30097929073974904355d0, & 0.794288809364887241d-2, & -0.6431910276847563d-4, & 0.22549830684374d-6, & -0.41220966195d-9, & 0.44185282d-12, & -0.30123d-15, & 0.14d-18 /) real ( kind = 8 ), dimension(0:27) :: ab0as = (/ & 1.97755499723693067407d+0, & -0.1046024792004819485d-1, & 0.69680790253625366d-3, & -0.5898298299996599d-4, & 0.577164455305320d-5, & -0.61523013365756d-6, & 0.6785396884767d-7, & -0.723062537907d-8, & 0.63306627365d-9, & -0.989453793d-11, & -0.1681980530d-10, & 0.673799551d-11, & -0.200997939d-11, & 0.54055903d-12, & -0.13816679d-12, & 0.3422205d-13, & -0.826686d-14, & 0.194566d-14, & -0.44268d-15, & 0.9562d-16, & -0.1883d-16, & 0.301d-17, & -0.19d-18, & -0.14d-18, & 0.11d-18, & -0.4d-19, & 0.2d-19, & -0.1d-19 /) real ( kind = 8 ) abram0 real ( kind = 8 ) asln real ( kind = 8 ) asval real ( kind = 8 ) cheval real ( kind = 8 ) fval real ( kind = 8 ) gval real ( kind = 8 ), parameter :: gval0 = 0.13417650264770070909D+00 real ( kind = 8 ), parameter :: half = 0.5D+00 real ( kind = 8 ) hval real ( kind = 8 ), parameter :: lnxmin = -708.3964D+00 integer, parameter :: nterma = 22 integer, parameter :: ntermf = 8 integer, parameter :: ntermg = 8 integer, parameter :: ntermh = 8 real ( kind = 8 ), parameter :: onerpi = 0.56418958354775628695D+00 real ( kind = 8 ), parameter :: rt3bpi = 0.97720502380583984317D+00 real ( kind = 8 ), parameter :: rtpib2 = 0.88622692545275801365D+00 real ( kind = 8 ), parameter :: six = 6.0D+00 real ( kind = 8 ) t real ( kind = 8 ), parameter :: three = 3.0D+00 real ( kind = 8 ), parameter :: two = 2.0D+00 real ( kind = 8 ) v real ( kind = 8 ) x real ( kind = 8 ), parameter :: xlow1 = 1.490116D-08 real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 x = xvalue if ( x < zero ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ABRAM0 - Fatal error!' write ( *, '(a)' ) ' Argument X < 0.' abram0 = zero else if ( x == zero ) then abram0 = rtpib2 else if ( x < xlow1 ) then abram0 = rtpib2 + x * ( log ( x ) - gval0 ) else if ( x <= two ) then t = ( x * x / two - half ) - half fval = cheval ( ntermf, ab0f, t ) gval = cheval ( ntermg, ab0g, t ) hval = cheval ( ntermh, ab0h, t ) abram0 = fval / onerpi + x * ( log ( x ) * hval - gval ) else v = three * ( ( x / two ) ** ( two / three ) ) t = ( six / v - half ) - half asval = cheval ( nterma, ab0as, t ) asln = log ( asval / rt3bpi ) - v if ( asln < lnxmin ) then abram0 = zero else abram0 = exp ( asln ) end if end if return end subroutine abram0_values ( n_data, x, fx ) !*****************************************************************************80 ! !! ABRAM0_VALUES returns some values of the Abramowitz0 function. ! ! Discussion: ! ! The function is defined by: ! ! ABRAM0(x) = Integral ( 0 <= t < infinity ) exp ( -t^2 - x / t ) dt ! ! The data was reported by McLeod. ! ! Modified: ! ! 21 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.87377726306985360531D+00, & 0.84721859650456925922D+00, & 0.77288934483988301615D+00, & 0.59684345853450151603D+00, & 0.29871735283675888392D+00, & 0.15004596450516388138D+00, & 0.11114662419157955096D+00, & 0.83909567153151897766D-01, & 0.56552321717943417515D-01, & 0.49876496603033790206D-01, & 0.44100889219762791328D-01, & 0.19738535180254062496D-01, & 0.86193088287161479900D-02, & 0.40224788162540127227D-02, & 0.19718658458164884826D-02, & 0.10045868340133538505D-02, & 0.15726917263304498649D-03, & 0.10352666912350263437D-04, & 0.91229759190956745069D-06, & 0.25628287737952698742D-09 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0019531250D+00, & 0.0078125000D+00, & 0.0312500000D+00, & 0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 1.2500000000D+00, & 1.5000000000D+00, & 1.8750000000D+00, & 2.0000000000D+00, & 2.1250000000D+00, & 3.0000000000D+00, & 4.0000000000D+00, & 5.0000000000D+00, & 6.0000000000D+00, & 7.0000000000D+00, & 10.0000000000D+00, & 15.0000000000D+00, & 20.0000000000D+00, & 40.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function abram1 ( xvalue ) !*****************************************************************************80 ! !! ABRAM1 evaluates the Abramowitz function of order 1. ! ! Discussion: ! ! The function is defined by: ! ! ABRAM1(x) = Integral ( 0 <= t < infinity ) t * exp ( -t^2 - x / t ) dt ! ! The code uses Chebyshev expansions with the coefficients ! given to an accuracy of 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 07 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) ABRAM1, the value of the function. ! implicit none real ( kind = 8 ) ab1as(0:27) real ( kind = 8 ) ab1f(0:9) real ( kind = 8 ) ab1g(0:8) real ( kind = 8 ) ab1h(0:8) real ( kind = 8 ) abram1 real ( kind = 8 ) asln real ( kind = 8 ) asval real ( kind = 8 ) cheval real ( kind = 8 ) fval real ( kind = 8 ) gval real ( kind = 8 ), parameter :: half = 0.5D+00 real ( kind = 8 ) hval real ( kind = 8 ) lnxmin integer, parameter :: nterma = 23 integer, parameter :: ntermf = 9 integer, parameter :: ntermg = 8 integer, parameter :: ntermh = 8 real ( kind = 8 ), parameter :: one = 1.0D+00 real ( kind = 8 ) onerpi real ( kind = 8 ), parameter :: rt3bpi = 0.97720502380583984317D+00 real ( kind = 8 ), parameter :: six = 6.0D+00 real ( kind = 8 ) t real ( kind = 8 ), parameter :: three = 3.0D+00 real ( kind = 8 ), parameter :: two = 2.0D+00 real ( kind = 8 ) v real ( kind = 8 ) x real ( kind = 8 ) xlow real ( kind = 8 ) xlow1 real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 data ab1f/1.47285192577978807369d0, & 0.10903497570168956257d0, & -0.12430675360056569753d0, & 0.306197946853493315d-2, & -0.2218410323076511d-4, & 0.6989978834451d-7, & -0.11597076444d-9, & 0.11389776d-12, & -0.7173d-16, & 0.3d-19/ data ab1g/0.39791277949054503528d0, & -0.29045285226454720849d0, & 0.1048784695465363504d-1, & -0.10249869522691336d-3, & 0.41150279399110d-6, & -0.83652638940d-9, & 0.97862595d-12, & -0.71868d-15, & 0.35d-18/ data ab1h/0.84150292152274947030d0, & -0.7790050698774143395d-1, & 0.133992455878390993d-2, & -0.808503907152788d-5, & 0.2261858281728d-7, & -0.3441395838d-10, & 0.3159858d-13, & -0.1884d-16, & 0.1d-19/ data ab1as(0)/ 2.13013643429065549448d0/ data ab1as(1)/ 0.6371526795218539933d-1/ data ab1as(2)/ -0.129334917477510647d-2/ data ab1as(3)/ 0.5678328753228265d-4/ data ab1as(4)/ -0.279434939177646d-5/ data ab1as(5)/ 0.5600214736787d-7/ data ab1as(6)/ 0.2392009242798d-7/ data ab1as(7)/ -0.750984865009d-8/ data ab1as(8)/ 0.173015330776d-8/ data ab1as(9)/ -0.36648877955d-9/ data ab1as(10)/ 0.7520758307d-10/ data ab1as(11)/-0.1517990208d-10/ data ab1as(12)/ 0.301713710d-11/ data ab1as(13)/-0.58596718d-12/ data ab1as(14)/ 0.10914455d-12/ data ab1as(15)/-0.1870536d-13/ data ab1as(16)/ 0.262542d-14/ data ab1as(17)/-0.14627d-15/ data ab1as(18)/-0.9500d-16/ data ab1as(19)/ 0.5873d-16/ data ab1as(20)/-0.2420d-16/ data ab1as(21)/ 0.868d-17/ data ab1as(22)/-0.290d-17/ data ab1as(23)/ 0.93d-18/ data ab1as(24)/-0.29d-18/ data ab1as(25)/ 0.9d-19/ data ab1as(26)/-0.3d-19/ data ab1as(27)/ 0.1d-19/ data onerpi/ 0.56418958354775628695d0/ ! ! Machine-dependent constants (suitable for IEEE machines) ! data xlow,xlow1,lnxmin/1.11023d-16,1.490116d-8,-708.3964d0/ x = xvalue if ( x < zero ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ABRAM1 - Fatal error!' write ( *, '(a)' ) ' Argument X < 0.' abram1 = zero else if ( x == zero ) then abram1 = half else if ( x < xlow ) then abram1 = half else if ( x < xlow1 ) then abram1 = ( one - x / onerpi - x * x * log ( x ) ) * half else if ( x <= two ) then t = ( x * x / two - half ) - half fval = cheval ( ntermf, ab1f, t ) gval = cheval ( ntermg, ab1g, t ) hval = cheval ( ntermh, ab1h, t ) abram1 = fval - x * ( gval / onerpi + x * log ( x ) * hval ) else v = three * ( ( x / two ) ** ( two / three ) ) t = ( six / v - half ) - half asval = cheval ( nterma, ab1as, t ) asln = log ( asval * sqrt ( v / three ) / rt3bpi ) - v if ( asln < lnxmin ) then abram1 = zero else abram1 = exp ( asln ) end if end if return end subroutine abram1_values ( n_data, x, fx ) !*****************************************************************************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. ! ! Modified: ! ! 21 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.49828219848799921792D+00, & 0.49324391773047288556D+00, & 0.47431612784691234649D+00, & 0.41095983258760410149D+00, & 0.25317617388227035867D+00, & 0.14656338138597777543D+00, & 0.11421547056018366587D+00, & 0.90026307383483764795D-01, & 0.64088214170742303375D-01, & 0.57446614314166191085D-01, & 0.51581624564800730959D-01, & 0.25263719555776416016D-01, & 0.11930803330196594536D-01, & 0.59270542280915272465D-02, & 0.30609215358017829567D-02, & 0.16307382136979552833D-02, & 0.28371851916959455295D-03, & 0.21122150121323238154D-04, & 0.20344578892601627337D-05, & 0.71116517236209642290D-09 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0019531250D+00, & 0.0078125000D+00, & 0.0312500000D+00, & 0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 1.2500000000D+00, & 1.5000000000D+00, & 1.8750000000D+00, & 2.0000000000D+00, & 2.1250000000D+00, & 3.0000000000D+00, & 4.0000000000D+00, & 5.0000000000D+00, & 6.0000000000D+00, & 7.0000000000D+00, & 10.0000000000D+00, & 15.0000000000D+00, & 20.0000000000D+00, & 40.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function abram2 ( xvalue ) !*****************************************************************************80 ! !! ABRAM2 evaluates the Abramowitz function of order 2. ! ! Discussion: ! ! The function is defined by: ! ! ABRAM2(x) = Integral ( 0 <= t < infinity ) t^2 * exp ( -t^2 - x / t ) dt ! ! The code uses Chebyshev expansions with the coefficients ! given to an accuracy of 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 07 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) ABRAM2, the value of the function. ! implicit none real ( kind = 8 ) abram2 real ( kind = 8 ) cheval real ( kind = 8 ), parameter :: half = 0.5D+00 integer, parameter :: nterma = 23 integer, parameter :: ntermf = 9 integer, parameter :: ntermg = 8 integer, parameter :: ntermh = 7 real ( kind = 8 ), parameter :: six = 6.0D+00 real ( kind = 8 ), parameter :: three = 3.0D+00 real ( kind = 8 ), parameter :: two = 2.0D+00 real ( kind = 8 ) x real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 real ( kind = 8 ) ab2f(0:9),ab2g(0:8),ab2h(0:7),ab2as(0:26), & asln,asval,fval,gval,hval,lnxmin, & onerpi,rtpib4,rt3bpi,t, & v,xlow,xlow1 data ab2f/1.03612162804243713846d0, & 0.19371246626794570012d0, & -0.7258758839233007378d-1, & 0.174790590864327399d-2, & -0.1281223233756549d-4, & 0.4115018153651d-7, & -0.6971047256d-10, & 0.6990183d-13, & -0.4492d-16, & 0.2d-19/ data ab2g/1.46290157198630741150d0, & 0.20189466883154014317d0, & -0.2908292087997129022d-1, & 0.47061049035270050d-3, & -0.257922080359333d-5, & 0.656133712946d-8, & -0.914110203d-11, & 0.774276d-14, & -0.429d-17/ data ab2h/0.30117225010910488881d0, & -0.1588667818317623783d-1, & 0.19295936935584526d-3, & -0.90199587849300d-6, & 0.206105041837d-8, & -0.265111806d-11, & 0.210864d-14, & -0.111d-17/ data ab2as(0)/ 2.46492325304334856893d0/ data ab2as(1)/ 0.23142797422248905432d0/ data ab2as(2)/ -0.94068173010085773d-3/ data ab2as(3)/ 0.8290270038089733d-4/ data ab2as(4)/ -0.883894704245866d-5/ data ab2as(5)/ 0.106638543567985d-5/ data ab2as(6)/ -0.13991128538529d-6/ data ab2as(7)/ 0.1939793208445d-7/ data ab2as(8)/ -0.277049938375d-8/ data ab2as(9)/ 0.39590687186d-9/ data ab2as(10)/-0.5408354342d-10/ data ab2as(11)/ 0.635546076d-11/ data ab2as(12)/-0.38461613d-12/ data ab2as(13)/-0.11696067d-12/ data ab2as(14)/ 0.6896671d-13/ data ab2as(15)/-0.2503113d-13/ data ab2as(16)/ 0.785586d-14/ data ab2as(17)/-0.230334d-14/ data ab2as(18)/ 0.64914d-15/ data ab2as(19)/-0.17797d-15/ data ab2as(20)/ 0.4766d-16/ data ab2as(21)/-0.1246d-16/ data ab2as(22)/ 0.316d-17/ data ab2as(23)/-0.77d-18/ data ab2as(24)/ 0.18d-18/ data ab2as(25)/-0.4d-19/ data ab2as(26)/ 0.1d-19/ data rt3bpi/ 0.97720502380583984317d0/ data rtpib4/ 0.44311346272637900682d0/ data onerpi/ 0.56418958354775628695d0/ ! ! Machine-dependent constants (suitable for IEEE machines) ! data xlow,xlow1,lnxmin/2.22045d-16,1.490116d-8,-708.3964d0/ x = xvalue if ( x < zero ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ABRAM2 - Fatal error!' write ( *, '(a)' ) ' Argument X < 0.' abram2 = zero else if ( x == zero ) then abram2 = rtpib4 else if ( x < xlow ) then abram2 = rtpib4 else if ( x < xlow1 ) then abram2 = rtpib4 - half * x + x * x * x * log ( x ) / six else if ( x <= 2.0D+00 ) then t = ( x * x / two - half ) - half fval = cheval ( ntermf, ab2f, t ) gval = cheval ( ntermg, ab2g, t ) hval = cheval ( ntermh, ab2h, t ) abram2 = fval / onerpi + x * ( x * x * log ( x ) * hval - gval ) else v = three * ( ( x / two ) ** ( two / three ) ) t = ( six / v - half ) - half asval = cheval ( nterma, ab2as, t ) asln = log ( asval / rt3bpi ) + log ( v / three ) - v if ( asln < lnxmin ) then abram2 = zero else abram2 = exp ( asln ) end if end if return end subroutine abram2_values ( n_data, x, fx ) !*****************************************************************************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. ! ! Modified: ! ! 22 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.44213858162107913430D+00, & 0.43923379545684026308D+00, & 0.42789857297092602234D+00, & 0.38652825661854504406D+00, & 0.26538204413231368110D+00, & 0.16848734838334595000D+00, & 0.13609200032513227112D+00, & 0.11070330027727917352D+00, & 0.82126019995530382267D-01, & 0.74538781999594581763D-01, & 0.67732034377612811390D-01, & 0.35641808698811851022D-01, & 0.17956589956618269083D-01, & 0.94058737143575370625D-02, & 0.50809356204299213556D-02, & 0.28149565414209719359D-02, & 0.53808696422559303431D-03, & 0.44821756380146327259D-04, & 0.46890678427324100410D-05, & 0.20161544850996420504D-08 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0019531250D+00, & 0.0078125000D+00, & 0.0312500000D+00, & 0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 1.2500000000D+00, & 1.5000000000D+00, & 1.8750000000D+00, & 2.0000000000D+00, & 2.1250000000D+00, & 3.0000000000D+00, & 4.0000000000D+00, & 5.0000000000D+00, & 6.0000000000D+00, & 7.0000000000D+00, & 10.0000000000D+00, & 15.0000000000D+00, & 20.0000000000D+00, & 40.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function airy_ai_int ( xvalue ) !*****************************************************************************80 ! !! AIRY_AI_INT calculates the integral of the Airy function Ai. ! ! Discussion: ! ! The function is defined by: ! ! AIRY_AI_INT(x) = Integral ( 0 <= t <= x ) Ai(t) dt ! ! The program uses Chebyshev expansions, the coefficients of which ! are given to 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 07 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) AIRY_AI_INT, the value of the function. ! implicit none real ( kind = 8 ) aaint1(0:25) real ( kind = 8 ) aaint2(0:21) real ( kind = 8 ) aaint3(0:40) real ( kind = 8 ) aaint4(0:17) real ( kind = 8 ) aaint5(0:17) real ( kind = 8 ) airy_ai_int real ( kind = 8 ) airzer real ( kind = 8 ) arg real ( kind = 8 ) cheval real ( kind = 8 ), parameter :: eight = 8.0D+00 real ( kind = 8 ) forty1 real ( kind = 8 ), parameter :: four = 4.0D+00 real ( kind = 8 ) fr996 real ( kind = 8 ) gval real ( kind = 8 ) hval real ( kind = 8 ) nine real ( kind = 8 ) ninhun integer, parameter :: nterm1 = 22 integer, parameter :: nterm2 = 17 integer, parameter :: nterm3 = 37 integer nterm4 integer nterm5 real ( kind = 8 ), parameter :: one = 1.0D+00 real ( kind = 8 ) piby4 real ( kind = 8 ) pitim6 real ( kind = 8 ) rt2b3p real ( kind = 8 ) t real ( kind = 8 ) temp real ( kind = 8 ), parameter :: three = 3.0D+00 real ( kind = 8 ), parameter :: two = 2.0D+00 real ( kind = 8 ) x real ( kind = 8 ) xhigh1 real ( kind = 8 ) xlow1 real ( kind = 8 ) xneg1 real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 real ( kind = 8 ) z data aaint1(0)/ 0.37713517694683695526d0/ data aaint1(1)/ -0.13318868432407947431d0/ data aaint1(2)/ 0.3152497374782884809d-1/ data aaint1(3)/ -0.318543076436574077d-2/ data aaint1(4)/ -0.87398764698621915d-3/ data aaint1(5)/ 0.46699497655396971d-3/ data aaint1(6)/ -0.9544936738983692d-4/ data aaint1(7)/ 0.542705687156716d-5/ data aaint1(8)/ 0.239496406252188d-5/ data aaint1(9)/ -0.75690270205649d-6/ data aaint1(10)/ 0.9050138584518d-7/ data aaint1(11)/ 0.320529456043d-8/ data aaint1(12)/-0.303825536444d-8/ data aaint1(13)/ 0.48900118596d-9/ data aaint1(14)/-0.1839820572d-10/ data aaint1(15)/-0.711247519d-11/ data aaint1(16)/ 0.151774419d-11/ data aaint1(17)/-0.10801922d-12/ data aaint1(18)/-0.963542d-14/ data aaint1(19)/ 0.313425d-14/ data aaint1(20)/-0.29446d-15/ data aaint1(21)/-0.477d-17/ data aaint1(22)/ 0.461d-17/ data aaint1(23)/-0.53d-18/ data aaint1(24)/ 0.1d-19/ data aaint1(25)/ 0.1d-19/ data aaint2(0)/ 1.92002524081984009769d0/ data aaint2(1)/ -0.4220049417256287021d-1/ data aaint2(2)/ -0.239457722965939223d-2/ data aaint2(3)/ -0.19564070483352971d-3/ data aaint2(4)/ -0.1547252891056112d-4/ data aaint2(5)/ -0.140490186137889d-5/ data aaint2(6)/ -0.12128014271367d-6/ data aaint2(7)/ -0.1179186050192d-7/ data aaint2(8)/ -0.104315578788d-8/ data aaint2(9)/ -0.10908209293d-9/ data aaint2(10)/-0.929633045d-11/ data aaint2(11)/-0.110946520d-11/ data aaint2(12)/-0.7816483d-13/ data aaint2(13)/-0.1319661d-13/ data aaint2(14)/-0.36823d-15/ data aaint2(15)/-0.21505d-15/ data aaint2(16)/ 0.1238d-16/ data aaint2(17)/-0.557d-17/ data aaint2(18)/ 0.84d-18/ data aaint2(19)/-0.21d-18/ data aaint2(20)/ 0.4d-19/ data aaint2(21)/-0.1d-19/ data aaint3(0)/ 0.47985893264791052053d0/ data aaint3(1)/ -0.19272375126169608863d0/ data aaint3(2)/ 0.2051154129525428189d-1/ data aaint3(3)/ 0.6332000070732488786d-1/ data aaint3(4)/ -0.5093322261845754082d-1/ data aaint3(5)/ 0.1284424078661663016d-1/ data aaint3(6)/ 0.2760137088989479413d-1/ data aaint3(7)/ -0.1547066673866649507d-1/ data aaint3(8)/ -0.1496864655389316026d-1/ data aaint3(9)/ 0.336617614173574541d-2/ data aaint3(10)/ 0.530851163518892985d-2/ data aaint3(11)/ 0.41371226458555081d-3/ data aaint3(12)/-0.102490579926726266d-2/ data aaint3(13)/-0.32508221672025853d-3/ data aaint3(14)/ 0.8608660957169213d-4/ data aaint3(15)/ 0.6671367298120775d-4/ data aaint3(16)/ 0.449205999318095d-5/ data aaint3(17)/-0.670427230958249d-5/ data aaint3(18)/-0.196636570085009d-5/ data aaint3(19)/ 0.22229677407226d-6/ data aaint3(20)/ 0.22332222949137d-6/ data aaint3(21)/ 0.2803313766457d-7/ data aaint3(22)/-0.1155651663619d-7/ data aaint3(23)/-0.433069821736d-8/ data aaint3(24)/-0.6227777938d-10/ data aaint3(25)/ 0.26432664903d-9/ data aaint3(26)/ 0.5333881114d-10/ data aaint3(27)/-0.522957269d-11/ data aaint3(28)/-0.382229283d-11/ data aaint3(29)/-0.40958233d-12/ data aaint3(30)/ 0.11515622d-12/ data aaint3(31)/ 0.3875766d-13/ data aaint3(32)/ 0.140283d-14/ data aaint3(33)/-0.141526d-14/ data aaint3(34)/-0.28746d-15/ data aaint3(35)/ 0.923d-17/ data aaint3(36)/ 0.1224d-16/ data aaint3(37)/ 0.157d-17/ data aaint3(38)/-0.19d-18/ data aaint3(39)/-0.8d-19/ data aaint3(40)/-0.1d-19/ data aaint4/1.99653305828522730048d0, & -0.187541177605417759d-2, & -0.15377536280305750d-3, & -0.1283112967682349d-4, & -0.108128481964162d-5, & -0.9182131174057d-7, & -0.784160590960d-8, & -0.67292453878d-9, & -0.5796325198d-10, & -0.501040991d-11, & -0.43420222d-12, & -0.3774305d-13, & -0.328473d-14, & -0.28700d-15, & -0.2502d-16, & -0.220d-17, & -0.19d-18, & -0.2d-19/ data aaint5/1.13024602034465716133d0, & -0.464718064639872334d-2, & -0.35137413382693203d-3, & -0.2768117872545185d-4, & -0.222057452558107d-5, & -0.18089142365974d-6, & -0.1487613383373d-7, & -0.123515388168d-8, & -0.10310104257d-9, & -0.867493013d-11, & -0.73080054d-12, & -0.6223561d-13, & -0.525128d-14, & -0.45677d-15, & -0.3748d-16, & -0.356d-17, & -0.23d-18, & -0.4d-19/ data nine,forty1/ 9.0d0, 41.0d0/ data ninhun,fr996/ 900.0d0, 4996.0d0 / data piby4/0.78539816339744830962d0/ data pitim6/18.84955592153875943078d0/ data rt2b3p/0.46065886596178063902d0/ data airzer/0.35502805388781723926d0/ ! ! Machine-dependant constants (suitable for IEEE machines) ! data nterm4,nterm5/15,15/ data xlow1,xhigh1,xneg1/2.22045d-16,14.480884d0,-2.727134d10/ x = xvalue if ( x < xneg1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'AIRY_AI_INT - Fatal error!' write ( *, '(a)' ) ' X too negative for accurate computation.' airy_ai_int = -two / three return else if ( x < -eight ) then z = - ( x + x ) * sqrt ( -x ) / three arg = z + piby4 temp = nine * z * z t = ( fr996 - temp ) / ( ninhun + temp ) gval = cheval ( nterm4, aaint4, t ) hval = cheval ( nterm5, aaint5, t ) temp = gval * cos ( arg ) + hval * sin ( arg ) / z airy_ai_int = rt2b3p * temp / sqrt ( z ) - two / three else if ( x <= -xlow1 )then t = -x / four - one airy_ai_int = x * cheval ( nterm3, aaint3, t ) else if ( x < xlow1 ) then airy_ai_int = airzer * x else if ( x <= four ) then t = x / two - one airy_ai_int = cheval ( nterm1, aaint1, t ) * x else if ( x <= xhigh1 ) then z = ( x + x ) * sqrt ( x ) / three temp = three * z t = ( forty1 - temp ) / ( nine + temp ) temp = exp ( -z ) * cheval ( nterm2, aaint2, t ) / sqrt ( pitim6 * z ) airy_ai_int = one / three - temp else airy_ai_int = one / three end if return end subroutine airy_ai_int_values ( n_data, x, fx ) !*****************************************************************************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. ! ! Modified: ! ! 22 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & -0.75228838916610124300D+00, & -0.57348350185854889466D+00, & -0.76569840313421291743D+00, & -0.65181015505382467421D+00, & -0.55881974894471876922D+00, & -0.56902352870716815309D+00, & -0.47800749642926168100D+00, & -0.46567398346706861416D+00, & -0.96783140945618013679D-01, & -0.34683049857035607494D-03, & 0.34658366917927930790D-03, & 0.27657581846051227124D-02, & 0.14595330491185717833D+00, & 0.23631734191710977960D+00, & 0.33289264538612212697D+00, & 0.33318759129779422976D+00, & 0.33332945170523851439D+00, & 0.33333331724248357420D+00, & 0.33333333329916901594D+00, & 0.33333333333329380187D+00 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & -12.0000000000D+00, & -11.0000000000D+00, & -10.0000000000D+00, & -9.5000000000D+00, & -9.0000000000D+00, & -6.5000000000D+00, & -4.0000000000D+00, & -1.0000000000D+00, & -0.2500000000D+00, & -0.0009765625D+00, & 0.0009765625D+00, & 0.0078125000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 4.0000000000D+00, & 4.5000000000D+00, & 6.0000000000D+00, & 8.0000000000D+00, & 10.0000000000D+00, & 12.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function airy_bi_int ( xvalue ) !*****************************************************************************80 ! !! AIRY_BI_INT calculates the integral of the Airy function Bi. ! ! Discussion: ! ! The function is defined by: ! ! AIRY_BI_INT(x) = Integral ( 0 <= t <= x ) Bi(t) dt ! ! The program uses Chebyshev expansions, the coefficients of which ! are given to 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 07 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) AIRY_BI_INT, the value of the function. ! implicit none real ( kind = 8 ) abint1(0:36) real ( kind = 8 ) abint2(0:37) real ( kind = 8 ) abint3(0:37) real ( kind = 8 ) abint4(0:20) real ( kind = 8 ) abint5(0:20) real ( kind = 8 ) airy_bi_int real ( kind = 8 ) cheval real ( kind = 8 ), parameter :: eight = 8.0D+00 real ( kind = 8 ), parameter :: four = 4.0D+00 integer, parameter :: nterm1 = 33 integer, parameter :: nterm2 = 30 integer, parameter :: nterm3 = 34 integer nterm4 integer nterm5 real ( kind = 8 ), parameter :: one = 1.0D+00 real ( kind = 8 ), parameter :: three = 3.0D+00 real ( kind = 8 ) x real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 real ( kind = 8 ) arg,birzer,f1,f2,nine,ninhun, & onept5,piby4,rt2b3p,sixten,seven,t,temp, & thr644,xlow1,xhigh1,xmax,xneg1, & z data abint1(0)/ 0.38683352445038543350d0/ data abint1(1)/ -0.8823213550888908821d-1/ data abint1(2)/ 0.21463937440355429239d0/ data abint1(3)/ -0.4205347375891315126d-1/ data abint1(4)/ 0.5932422547496086771d-1/ data abint1(5)/ -0.840787081124270210d-2/ data abint1(6)/ 0.871824772778487955d-2/ data abint1(7)/ -0.12191600199613455d-3/ data abint1(8)/ 0.44024821786023234d-3/ data abint1(9)/ 0.27894686666386678d-3/ data abint1(10)/-0.7052804689785537d-4/ data abint1(11)/ 0.5901080066770100d-4/ data abint1(12)/-0.1370862587982142d-4/ data abint1(13)/ 0.505962573749073d-5/ data abint1(14)/-0.51598837766735d-6/ data abint1(15)/ 0.397511312349d-8/ data abint1(16)/ 0.9524985978055d-7/ data abint1(17)/-0.3681435887321d-7/ data abint1(18)/ 0.1248391688136d-7/ data abint1(19)/-0.249097619137d-8/ data abint1(20)/ 0.31775245551d-9/ data abint1(21)/ 0.5434365270d-10/ data abint1(22)/-0.4024566915d-10/ data abint1(23)/ 0.1393855527d-10/ data abint1(24)/-0.303817509d-11/ data abint1(25)/ 0.40809511d-12/ data abint1(26)/ 0.1634116d-13/ data abint1(27)/-0.2683809d-13/ data abint1(28)/ 0.896641d-14/ data abint1(29)/-0.183089d-14/ data abint1(30)/ 0.21333d-15/ data abint1(31)/ 0.1108d-16/ data abint1(32)/-0.1276d-16/ data abint1(33)/ 0.363d-17/ data abint1(34)/-0.62d-18/ data abint1(35)/ 0.5d-19/ data abint1(36)/ 0.1d-19/ data abint2(0)/ 2.04122078602516135181d0/ data abint2(1)/ 0.2124133918621221230d-1/ data abint2(2)/ 0.66617599766706276d-3/ data abint2(3)/ 0.3842047982808254d-4/ data abint2(4)/ 0.362310366020439d-5/ data abint2(5)/ 0.50351990115074d-6/ data abint2(6)/ 0.7961648702253d-7/ data abint2(7)/ 0.717808442336d-8/ data abint2(8)/ -0.267770159104d-8/ data abint2(9)/ -0.168489514699d-8/ data abint2(10)/-0.36811757255d-9/ data abint2(11)/ 0.4757128727d-10/ data abint2(12)/ 0.5263621945d-10/ data abint2(13)/ 0.778973500d-11/ data abint2(14)/-0.460546143d-11/ data abint2(15)/-0.183433736d-11/ data abint2(16)/ 0.32191249d-12/ data abint2(17)/ 0.29352060d-12/ data abint2(18)/-0.1657935d-13/ data abint2(19)/-0.4483808d-13/ data abint2(20)/ 0.27907d-15/ data abint2(21)/ 0.711921d-14/ data abint2(22)/-0.1042d-16/ data abint2(23)/-0.119591d-14/ data abint2(24)/ 0.4606d-16/ data abint2(25)/ 0.20884d-15/ data abint2(26)/-0.2416d-16/ data abint2(27)/-0.3638d-16/ data abint2(28)/ 0.863d-17/ data abint2(29)/ 0.591d-17/ data abint2(30)/-0.256d-17/ data abint2(31)/-0.77d-18/ data abint2(32)/ 0.66d-18/ data abint2(33)/ 0.3d-19/ data abint2(34)/-0.15d-18/ data abint2(35)/ 0.2d-19/ data abint2(36)/ 0.3d-19/ data abint2(37)/-0.1d-19/ data abint3(0)/ 0.31076961598640349251d0/ data abint3(1)/ -0.27528845887452542718d0/ data abint3(2)/ 0.17355965706136543928d0/ data abint3(3)/ -0.5544017909492843130d-1/ data abint3(4)/ -0.2251265478295950941d-1/ data abint3(5)/ 0.4107347447812521894d-1/ data abint3(6)/ 0.984761275464262480d-2/ data abint3(7)/ -0.1555618141666041932d-1/ data abint3(8)/ -0.560871870730279234d-2/ data abint3(9)/ 0.246017783322230475d-2/ data abint3(10)/ 0.165740392292336978d-2/ data abint3(11)/-0.3277587501435402d-4/ data abint3(12)/-0.24434680860514925d-3/ data abint3(13)/-0.5035305196152321d-4/ data abint3(14)/ 0.1630264722247854d-4/ data abint3(15)/ 0.851914057780934d-5/ data abint3(16)/ 0.29790363004664d-6/ data abint3(17)/-0.64389707896401d-6/ data abint3(18)/-0.15046988145803d-6/ data abint3(19)/ 0.1587013535823d-7/ data abint3(20)/ 0.1276766299622d-7/ data abint3(21)/ 0.140578534199d-8/ data abint3(22)/-0.46564739741d-9/ data abint3(23)/-0.15682748791d-9/ data abint3(24)/-0.403893560d-11/ data abint3(25)/ 0.666708192d-11/ data abint3(26)/ 0.128869380d-11/ data abint3(27)/-0.6968663d-13/ data abint3(28)/-0.6254319d-13/ data abint3(29)/-0.718392d-14/ data abint3(30)/ 0.115296d-14/ data abint3(31)/ 0.42276d-15/ data abint3(32)/ 0.2493d-16/ data abint3(33)/-0.971d-17/ data abint3(34)/-0.216d-17/ data abint3(35)/-0.2d-19/ data abint3(36)/ 0.6d-19/ data abint3(37)/ 0.1d-19/ data abint4(0)/ 1.99507959313352047614d0/ data abint4(1)/ -0.273736375970692738d-2/ data abint4(2)/ -0.30897113081285850d-3/ data abint4(3)/ -0.3550101982798577d-4/ data abint4(4)/ -0.412179271520133d-5/ data abint4(5)/ -0.48235892316833d-6/ data abint4(6)/ -0.5678730727927d-7/ data abint4(7)/ -0.671874810365d-8/ data abint4(8)/ -0.79811649857d-9/ data abint4(9)/ -0.9514271478d-10/ data abint4(10)/-0.1137468966d-10/ data abint4(11)/-0.136359969d-11/ data abint4(12)/-0.16381418d-12/ data abint4(13)/-0.1972575d-13/ data abint4(14)/-0.237844d-14/ data abint4(15)/-0.28752d-15/ data abint4(16)/-0.3475d-16/ data abint4(17)/-0.422d-17/ data abint4(18)/-0.51d-18/ data abint4(19)/-0.6d-19/ data abint4(20)/-0.1d-19/ data abint5(0)/ 1.12672081961782566017d0/ data abint5(1)/ -0.671405567525561198d-2/ data abint5(2)/ -0.69812918017832969d-3/ data abint5(3)/ -0.7561689886425276d-4/ data abint5(4)/ -0.834985574510207d-5/ data abint5(5)/ -0.93630298232480d-6/ data abint5(6)/ -0.10608556296250d-6/ data abint5(7)/ -0.1213128916741d-7/ data abint5(8)/ -0.139631129765d-8/ data abint5(9)/ -0.16178918054d-9/ data abint5(10)/-0.1882307907d-10/ data abint5(11)/-0.220272985d-11/ data abint5(12)/-0.25816189d-12/ data abint5(13)/-0.3047964d-13/ data abint5(14)/-0.358370d-14/ data abint5(15)/-0.42831d-15/ data abint5(16)/-0.4993d-16/ data abint5(17)/-0.617d-17/ data abint5(18)/-0.68d-18/ data abint5(19)/-0.10d-18/ data abint5(20)/-0.1d-19/ data onept5/ 1.5d0 / data seven/ 7.0d0 / data nine,sixten/ 9.0d0 , 16.0d0 / data ninhun,thr644/900.0d0 , 3644.0d0 / data piby4/0.78539816339744830962d0/ data rt2b3p/0.46065886596178063902d0/ data birzer/0.61492662744600073515d0/ ! ! Machine-dependent parameters (suitable for IEEE machines) ! data nterm4,nterm5/17,17/ data xlow1,xhigh1/2.22044604d-16,104.587632d0/ data xneg1,xmax/-2.727134d10,1.79d308/ x = xvalue if ( x < xneg1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'AIRY_BI_INT - Warning!' write ( *, '(a)' ) ' Argument is too negative for accurate computation.' airy_bi_int = zero else if ( x < -seven ) then z = - ( x + x ) * sqrt ( -x ) / three arg = z + piby4 temp = nine * z * z t = ( thr644 - temp ) / ( ninhun + temp ) f1 = cheval ( nterm4, abint4, t ) * sin ( arg ) f2 = cheval ( nterm5, abint5, t ) * cos ( arg ) / z airy_bi_int = ( f2 - f1 ) * rt2b3p / sqrt ( z ) else if ( x <= -xlow1 ) then t = - ( x + x ) / seven - one airy_bi_int = x * cheval ( nterm3, abint3, t ) else if ( x < xlow1 ) then airy_bi_int = birzer * x else if ( x <= eight ) then t = x / four - one airy_bi_int = x * exp ( onept5 * x ) * cheval ( nterm1, abint1, t ) else if ( x <= xhigh1 ) then t = sixten * sqrt ( eight / x ) / x - one z = ( x + x ) * sqrt ( x ) / three temp = rt2b3p * cheval ( nterm2, abint2, t ) / sqrt ( z ) temp = z + log ( temp ) airy_bi_int = exp ( temp ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'AIRY_BI_INT - Warning!' write ( *, '(a)' ) ' Argument is too large for accurate computation.' airy_bi_int = xmax end if return end subroutine airy_bi_int_values ( n_data, x, fx ) !*****************************************************************************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. ! ! Modified: ! ! 23 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.17660819031554631869D-01, & -0.15040424806140020451D-01, & 0.14756446293227661920D-01, & -0.11847304264848446271D+00, & -0.64916741266165856037D-01, & 0.97260832464381044540D-01, & 0.50760058495287539119D-01, & -0.37300500963429492179D+00, & -0.13962988442666578531D+00, & -0.12001735266723296160D-02, & 0.12018836117890354598D-02, & 0.36533846550952011043D+00, & 0.87276911673800812196D+00, & 0.48219475263803429675D+02, & 0.44006525804904178439D+06, & 0.17608153976228301458D+07, & 0.73779211705220007228D+07, & 0.14780980310740671617D+09, & 0.97037614223613433849D+11, & 0.11632737638809878460D+15 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & -12.0000000000D+00, & -10.0000000000D+00, & -8.0000000000D+00, & -7.5000000000D+00, & -7.0000000000D+00, & -6.5000000000D+00, & -4.0000000000D+00, & -1.0000000000D+00, & -0.2500000000D+00, & -0.0019531250D+00, & 0.0019531250D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 4.0000000000D+00, & 8.0000000000D+00, & 8.5000000000D+00, & 9.0000000000D+00, & 10.0000000000D+00, & 12.0000000000D+00, & 14.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function airy_gi ( xvalue ) !*****************************************************************************80 ! !! AIRY_GI computes the modified Airy function Gi(x). ! ! Discussion: ! ! The function is defined by: ! ! AIRY_GI(x) = Integral ( 0 <= t < infinity ) sin ( x*t+t^3/3) dt / pi ! ! The approximation uses Chebyshev expansions with the coefficients ! given to 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 07 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) AIRY_GI, the value of the function. ! implicit none real ( kind = 8 ) airy_gi real ( kind = 8 ) cheval real ( kind = 8 ), parameter :: four = 4.0D+00 integer, parameter :: nterm1 = 28 integer, parameter :: nterm2 = 23 integer, parameter :: nterm3 = 39 integer nterm4 integer nterm5 integer nterm6 real ( kind = 8 ), parameter :: one = 1.0D+00 real ( kind = 8 ), parameter :: three = 3.0D+00 real ( kind = 8 ) x real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 real ( kind = 8 ) argip1(0:30),argip2(0:29),argin1(0:42), & arbin1(0:10),arbin2(0:11),arhin1(0:15), & bi,cheb1,cheb2,cosz,five,five14, & gizero,minate,nine,onebpi,one76,one024,piby4, & rtpiin,seven,seven2,sinz,t,temp,twelhu,twent8, & xcube,xhigh1,xhigh2,xhigh3,xlow1,xminus, & zeta data argip1(0)/ 0.26585770795022745082d0/ data argip1(1)/ -0.10500333097501922907d0/ data argip1(2)/ 0.841347475328454492d-2/ data argip1(3)/ 0.2021067387813439541d-1/ data argip1(4)/ -0.1559576113863552234d-1/ data argip1(5)/ 0.564342939043256481d-2/ data argip1(6)/ -0.59776844826655809d-3/ data argip1(7)/ -0.42833850264867728d-3/ data argip1(8)/ 0.22605662380909027d-3/ data argip1(9)/ -0.3608332945592260d-4/ data argip1(10)/-0.785518988788901d-5/ data argip1(11)/ 0.473252480746370d-5/ data argip1(12)/-0.59743513977694d-6/ data argip1(13)/-0.15917609165602d-6/ data argip1(14)/ 0.6336129065570d-7/ data argip1(15)/-0.276090232648d-8/ data argip1(16)/-0.256064154085d-8/ data argip1(17)/ 0.47798676856d-9/ data argip1(18)/ 0.4488131863d-10/ data argip1(19)/-0.2346508882d-10/ data argip1(20)/ 0.76839085d-12/ data argip1(21)/ 0.73227985d-12/ data argip1(22)/-0.8513687d-13/ data argip1(23)/-0.1630201d-13/ data argip1(24)/ 0.356769d-14/ data argip1(25)/ 0.25001d-15/ data argip1(26)/-0.10859d-15/ data argip1(27)/-0.158d-17/ data argip1(28)/ 0.275d-17/ data argip1(29)/-0.5d-19/ data argip1(30)/-0.6d-19/ data argip2(0)/ 2.00473712275801486391d0/ data argip2(1)/ 0.294184139364406724d-2/ data argip2(2)/ 0.71369249006340167d-3/ data argip2(3)/ 0.17526563430502267d-3/ data argip2(4)/ 0.4359182094029882d-4/ data argip2(5)/ 0.1092626947604307d-4/ data argip2(6)/ 0.272382418399029d-5/ data argip2(7)/ 0.66230900947687d-6/ data argip2(8)/ 0.15425323370315d-6/ data argip2(9)/ 0.3418465242306d-7/ data argip2(10)/ 0.728157724894d-8/ data argip2(11)/ 0.151588525452d-8/ data argip2(12)/ 0.30940048039d-9/ data argip2(13)/ 0.6149672614d-10/ data argip2(14)/ 0.1202877045d-10/ data argip2(15)/ 0.233690586d-11/ data argip2(16)/ 0.43778068d-12/ data argip2(17)/ 0.7996447d-13/ data argip2(18)/ 0.1494075d-13/ data argip2(19)/ 0.246790d-14/ data argip2(20)/ 0.37672d-15/ data argip2(21)/ 0.7701d-16/ data argip2(22)/ 0.354d-17/ data argip2(23)/-0.49d-18/ data argip2(24)/ 0.62d-18/ data argip2(25)/-0.40d-18/ data argip2(26)/-0.1d-19/ data argip2(27)/ 0.2d-19/ data argip2(28)/-0.3d-19/ data argip2(29)/ 0.1d-19/ data argin1(0)/ -0.20118965056732089130d0/ data argin1(1)/ -0.7244175303324530499d-1/ data argin1(2)/ 0.4505018923894780120d-1/ data argin1(3)/ -0.24221371122078791099d0/ data argin1(4)/ 0.2717884964361678294d-1/ data argin1(5)/ -0.5729321004818179697d-1/ data argin1(6)/ -0.18382107860337763587d0/ data argin1(7)/ 0.7751546082149475511d-1/ data argin1(8)/ 0.18386564733927560387d0/ data argin1(9)/ 0.2921504250185567173d-1/ data argin1(10)/-0.6142294846788018811d-1/ data argin1(11)/-0.2999312505794616238d-1/ data argin1(12)/ 0.585937118327706636d-2/ data argin1(13)/ 0.822221658497402529d-2/ data argin1(14)/ 0.132579817166846893d-2/ data argin1(15)/-0.96248310766565126d-3/ data argin1(16)/-0.45065515998211807d-3/ data argin1(17)/ 0.772423474325474d-5/ data argin1(18)/ 0.5481874134758052d-4/ data argin1(19)/ 0.1245898039742876d-4/ data argin1(20)/-0.246196891092083d-5/ data argin1(21)/-0.169154183545285d-5/ data argin1(22)/-0.16769153169442d-6/ data argin1(23)/ 0.9636509337672d-7/ data argin1(24)/ 0.3253314928030d-7/ data argin1(25)/ 0.5091804231d-10/ data argin1(26)/-0.209180453553d-8/ data argin1(27)/-0.41237387870d-9/ data argin1(28)/ 0.4163338253d-10/ data argin1(29)/ 0.3032532117d-10/ data argin1(30)/ 0.340580529d-11/ data argin1(31)/-0.88444592d-12/ data argin1(32)/-0.31639612d-12/ data argin1(33)/-0.1505076d-13/ data argin1(34)/ 0.1104148d-13/ data argin1(35)/ 0.246508d-14/ data argin1(36)/-0.3107d-16/ data argin1(37)/-0.9851d-16/ data argin1(38)/-0.1453d-16/ data argin1(39)/ 0.118d-17/ data argin1(40)/ 0.67d-18/ data argin1(41)/ 0.6d-19/ data argin1(42)/-0.1d-19/ data arbin1/1.99983763583586155980d0, & -0.8104660923669418d-4, & 0.13475665984689d-6, & -0.70855847143d-9, & 0.748184187d-11, & -0.12902774d-12, & 0.322504d-14, & -0.10809d-15, & 0.460d-17, & -0.24d-18, & 0.1d-19/ data arbin2/0.13872356453879120276d0, & -0.8239286225558228d-4, & 0.26720919509866d-6, & -0.207423685368d-8, & 0.2873392593d-10, & -0.60873521d-12, & 0.1792489d-13, & -0.68760d-15, & 0.3280d-16, & -0.188d-17, & 0.13d-18, & -0.1d-19/ data arhin1/1.99647720399779650525d0, & -0.187563779407173213d-2, & -0.12186470897787339d-3, & -0.814021609659287d-5, & -0.55050925953537d-6, & -0.3763008043303d-7, & -0.258858362365d-8, & -0.17931829265d-9, & -0.1245916873d-10, & -0.87171247d-12, & -0.6084943d-13, & -0.431178d-14, & -0.29787d-15, & -0.2210d-16, & -0.136d-17, & -0.14d-18/ data five,seven,minate/ 5.0d0, 7.0d0 , -8.0d0 / data nine,twent8,seven2/ 9.0d0, 28.0d0 , 72.0d0 / data one76,five14/ 176.0d0 , 514.0d0 / data one024,twelhu/ 1024.0d0, 1200.0d0 / data gizero/0.20497554248200024505d0/ data onebpi/0.31830988618379067154d0/ data piby4/0.78539816339744830962d0/ data rtpiin/0.56418958354775628695d0/ ! ! Machine-dependent constants (suitable for IEEE machines) ! data nterm4,nterm5,nterm6/9,10,14/ data xlow1,xhigh1/2.22045d-16,208063.8307d0/ data xhigh2,xhigh3/0.14274d308,-2097152.0d0/ x = xvalue if ( x < -xhigh1 * xhigh1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'AIRY_GI - Fatal error!' write ( *, '(a)' ) ' Argument too negative for accurate computation.' airy_gi = zero else if ( x <= xhigh3 ) then xminus = -x t = xminus * sqrt ( xminus ) zeta = ( t + t ) / three temp = rtpiin / sqrt ( sqrt ( xminus ) ) cosz = cos ( zeta + piby4 ) sinz = sin ( zeta + piby4 ) / zeta xcube = x * x * x bi = ( cosz + sinz * five / seven2 ) * temp t = ( xcube + twelhu ) / ( one76 - xcube ) airy_gi = bi + cheval ( nterm6, arhin1, t ) * onebpi / x else if ( x < minate ) then xminus = -x t = xminus * sqrt ( xminus ) zeta = ( t + t ) / three temp = rtpiin / sqrt ( sqrt ( xminus ) ) cosz = cos ( zeta + piby4 ) sinz = sin ( zeta + piby4 ) / zeta xcube = x * x * x t = - ( one024 / ( xcube ) + one ) cheb1 = cheval ( nterm4, arbin1, t ) cheb2 = cheval ( nterm5, arbin2, t ) bi = ( cosz * cheb1 + sinz * cheb2 ) * temp t = ( xcube + twelhu ) / ( one76 - xcube ) airy_gi = bi + cheval ( nterm6, arhin1, t ) * onebpi / x else if ( x <= -xlow1 ) then t = -( x + four ) / four airy_gi = cheval ( nterm3, argin1, t ) else if ( x < xlow1 ) then airy_gi = gizero else if ( x <= seven ) then t = ( nine * x - twent8 ) / ( x + twent8 ) airy_gi = cheval ( nterm1, argip1, t ) else if ( x <= xhigh1 ) then xcube = x * x * x t = ( twelhu - xcube ) / ( five14 + xcube ) airy_gi = onebpi * cheval ( nterm2, argip2, t ) / x else if ( x <= xhigh2 ) then airy_gi = onebpi / x else airy_gi = zero end if return end subroutine airy_gi_values ( n_data, x, fx ) !*****************************************************************************80 ! !! AIRY_GI_VALUES returns some values of the Airy Gi function. ! ! Discussion: ! ! The function is defined by: ! ! AIRY_GI(x) = Integral ( 0 <= t < infinity ) sin ( x*t+t^3/3) dt / pi ! ! The data was reported by McLeod. ! ! Modified: ! ! 24 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.20468308070040542435D+00, & 0.18374662832557904078D+00, & -0.11667221729601528265D+00, & 0.31466934902729557596D+00, & -0.37089040722426257729D+00, & -0.25293059772424019694D+00, & 0.28967410658692701936D+00, & -0.34644836492634090590D+00, & 0.28076035913873049496D+00, & 0.21814994508094865815D+00, & 0.20526679000810503329D+00, & 0.22123695363784773258D+00, & 0.23521843981043793760D+00, & 0.82834303363768729338D-01, & 0.45757385490989281893D-01, & 0.44150012014605159922D-01, & 0.39951133719508907541D-01, & 0.35467706833949671483D-01, & 0.31896005100679587981D-01, & 0.26556892713512410405D-01 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & -0.0019531250D+00, & -0.1250000000D+00, & -1.0000000000D+00, & -4.0000000000D+00, & -8.0000000000D+00, & -8.2500000000D+00, & -9.0000000000D+00, & -10.0000000000D+00, & -11.0000000000D+00, & -13.0000000000D+00, & 0.0019531250D+00, & 0.1250000000D+00, & 1.0000000000D+00, & 4.0000000000D+00, & 7.0000000000D+00, & 7.2500000000D+00, & 8.0000000000D+00, & 9.0000000000D+00, & 10.0000000000D+00, & 12.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function airy_hi ( xvalue ) !*****************************************************************************80 ! !! AIRY_HI computes the modified Airy function Hi(x). ! ! Discussion: ! ! The function is defined by: ! ! AIRY_HI(x) = Integral ( 0 <= t < infinity ) exp(x*t-t^3/3) dt / pi ! ! The approximation uses Chebyshev expansions with the coefficients ! given to 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 07 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) AIRY_HI, the value of the function. ! implicit none real ( kind = 8 ) airy_hi real ( kind = 8 ) cheval real ( kind = 8 ), parameter :: four = 4.0D+00 integer, parameter :: nterm1 = 29 integer, parameter :: nterm2 = 17 integer, parameter :: nterm3 = 22 integer nterm4 integer nterm5 real ( kind = 8 ), parameter :: one = 1.0D+00 real ( kind = 8 ), parameter :: three = 3.0D+00 real ( kind = 8 ), parameter :: two = 2.0D+00 real ( kind = 8 ) x real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 real ( kind = 8 ) arhip(0:31),arbip(0:23),argip1(0:29), & arhin1(0:21),arhin2(0:15), & bi,five14,gi,hizero,lnrtpi, & minate,onebpi,one76,seven,t,temp, & thre43,twelhu,twelve,xcube, & xhigh1,xlow1,xmax,xneg1,xneg2, & zeta data arhip(0)/ 1.24013562561762831114d0/ data arhip(1)/ 0.64856341973926535804d0/ data arhip(2)/ 0.55236252592114903246d0/ data arhip(3)/ 0.20975122073857566794d0/ data arhip(4)/ 0.12025669118052373568d0/ data arhip(5)/ 0.3768224931095393785d-1/ data arhip(6)/ 0.1651088671548071651d-1/ data arhip(7)/ 0.455922755211570993d-2/ data arhip(8)/ 0.161828480477635013d-2/ data arhip(9)/ 0.40841282508126663d-3/ data arhip(10)/0.12196479721394051d-3/ data arhip(11)/0.2865064098657610d-4/ data arhip(12)/0.742221556424344d-5/ data arhip(13)/0.163536231932831d-5/ data arhip(14)/0.37713908188749d-6/ data arhip(15)/0.7815800336008d-7/ data arhip(16)/0.1638447121370d-7/ data arhip(17)/0.319857665992d-8/ data arhip(18)/0.61933905307d-9/ data arhip(19)/0.11411161191d-9/ data arhip(20)/0.2064923454d-10/ data arhip(21)/0.360018664d-11/ data arhip(22)/0.61401849d-12/ data arhip(23)/0.10162125d-12/ data arhip(24)/0.1643701d-13/ data arhip(25)/0.259084d-14/ data arhip(26)/0.39931d-15/ data arhip(27)/0.6014d-16/ data arhip(28)/0.886d-17/ data arhip(29)/0.128d-17/ data arhip(30)/0.18d-18/ data arhip(31)/0.3d-19/ data arbip(0)/ 2.00582138209759064905d0/ data arbip(1)/ 0.294478449170441549d-2/ data arbip(2)/ 0.3489754514775355d-4/ data arbip(3)/ 0.83389733374343d-6/ data arbip(4)/ 0.3136215471813d-7/ data arbip(5)/ 0.167865306015d-8/ data arbip(6)/ 0.12217934059d-9/ data arbip(7)/ 0.1191584139d-10/ data arbip(8)/ 0.154142553d-11/ data arbip(9)/ 0.24844455d-12/ data arbip(10)/ 0.4213012d-13/ data arbip(11)/ 0.505293d-14/ data arbip(12)/-0.60032d-15/ data arbip(13)/-0.65474d-15/ data arbip(14)/-0.22364d-15/ data arbip(15)/-0.3015d-16/ data arbip(16)/ 0.959d-17/ data arbip(17)/ 0.616d-17/ data arbip(18)/ 0.97d-18/ data arbip(19)/-0.37d-18/ data arbip(20)/-0.21d-18/ data arbip(21)/-0.1d-19/ data arbip(22)/ 0.2d-19/ data arbip(23)/ 0.1d-19/ data argip1(0)/ 2.00473712275801486391d0/ data argip1(1)/ 0.294184139364406724d-2/ data argip1(2)/ 0.71369249006340167d-3/ data argip1(3)/ 0.17526563430502267d-3/ data argip1(4)/ 0.4359182094029882d-4/ data argip1(5)/ 0.1092626947604307d-4/ data argip1(6)/ 0.272382418399029d-5/ data argip1(7)/ 0.66230900947687d-6/ data argip1(8)/ 0.15425323370315d-6/ data argip1(9)/ 0.3418465242306d-7/ data argip1(10)/ 0.728157724894d-8/ data argip1(11)/ 0.151588525452d-8/ data argip1(12)/ 0.30940048039d-9/ data argip1(13)/ 0.6149672614d-10/ data argip1(14)/ 0.1202877045d-10/ data argip1(15)/ 0.233690586d-11/ data argip1(16)/ 0.43778068d-12/ data argip1(17)/ 0.7996447d-13/ data argip1(18)/ 0.1494075d-13/ data argip1(19)/ 0.246790d-14/ data argip1(20)/ 0.37672d-15/ data argip1(21)/ 0.7701d-16/ data argip1(22)/ 0.354d-17/ data argip1(23)/-0.49d-18/ data argip1(24)/ 0.62d-18/ data argip1(25)/-0.40d-18/ data argip1(26)/-0.1d-19/ data argip1(27)/ 0.2d-19/ data argip1(28)/-0.3d-19/ data argip1(29)/ 0.1d-19/ data arhin1(0)/ 0.31481017206423404116d0/ data arhin1(1)/ -0.16414499216588964341d0/ data arhin1(2)/ 0.6176651597730913071d-1/ data arhin1(3)/ -0.1971881185935933028d-1/ data arhin1(4)/ 0.536902830023331343d-2/ data arhin1(5)/ -0.124977068439663038d-2/ data arhin1(6)/ 0.24835515596994933d-3/ data arhin1(7)/ -0.4187024096746630d-4/ data arhin1(8)/ 0.590945437979124d-5/ data arhin1(9)/ -0.68063541184345d-6/ data arhin1(10)/ 0.6072897629164d-7/ data arhin1(11)/-0.367130349242d-8/ data arhin1(12)/ 0.7078017552d-10/ data arhin1(13)/ 0.1187894334d-10/ data arhin1(14)/-0.120898723d-11/ data arhin1(15)/ 0.1189656d-13/ data arhin1(16)/ 0.594128d-14/ data arhin1(17)/-0.32257d-15/ data arhin1(18)/-0.2290d-16/ data arhin1(19)/ 0.253d-17/ data arhin1(20)/ 0.9d-19/ data arhin1(21)/-0.2d-19/ data arhin2/1.99647720399779650525d0, & -0.187563779407173213d-2, & -0.12186470897787339d-3, & -0.814021609659287d-5, & -0.55050925953537d-6, & -0.3763008043303d-7, & -0.258858362365d-8, & -0.17931829265d-9, & -0.1245916873d-10, & -0.87171247d-12, & -0.6084943d-13, & -0.431178d-14, & -0.29787d-15, & -0.2210d-16, & -0.136d-17, & -0.14d-18/ data seven/ 7.0d0 / data minate,twelve,one76/ -8.0d0 , 12.0d0 , 176.0d0 / data thre43,five14,twelhu/ 343.0d0, 514.0d0, 1200.0d0 / data hizero/0.40995108496400049010d0/ data lnrtpi/0.57236494292470008707d0/ data onebpi/0.31830988618379067154d0/ ! ! Machine-dependent constants (suitable for IEEE machines) ! data nterm4,nterm5/19,14/ data xlow1,xhigh1/2.220446d-16,104.4175d0/ data xneg1,xneg2,xmax/-0.14274d308,-208063.831d0,1.79d308/ x = xvalue if ( xhigh1 < x ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'AIRY_HI - Fatal error!' write ( *, '(a)' ) ' Argument too large.' airy_hi = xmax return end if ! ! Code for x < 0.0 ! if ( x < zero ) then if ( x < minate ) then if ( x < xneg1 ) then airy_hi = zero else if ( x < xneg2 ) then temp = one airy_hi = - temp * onebpi / x else xcube = x * x * x t = ( xcube + twelhu ) / ( one76 - xcube ) temp = cheval ( nterm5, arhin2, t ) airy_hi = - temp * onebpi / x end if end if else if ( -xlow1 < x ) then airy_hi = hizero else t = ( four * x + twelve ) / ( x - twelve ) airy_hi = cheval ( nterm4, arhin1, t ) end if end if ! ! Code for x >= 0.0 ! else if ( x <= seven ) then if ( x < xlow1 ) then airy_hi = hizero else t = ( x + x ) / seven - one temp = ( x + x + x ) / two airy_hi = exp ( temp ) * cheval ( nterm1, arhip, t ) end if else xcube = x * x * x temp = sqrt ( xcube ) zeta = ( temp + temp ) / three t = two * ( sqrt ( thre43 / xcube ) ) - one temp = cheval ( nterm2, arbip, t ) temp = zeta + log ( temp ) - log ( x ) / four - lnrtpi bi = exp ( temp ) t = ( twelhu - xcube ) / ( xcube + five14 ) gi = cheval ( nterm3, argip1, t ) * onebpi / x airy_hi = bi - gi end if end if return end subroutine airy_hi_values ( n_data, x, fx ) !*****************************************************************************80 ! !! AIRY_HI_VALUES returns some values of the Airy Hi function. ! ! Discussion: ! ! The function is defined by: ! ! AIRY_HI(x) = Integral ( 0 <= t < infinity ) exp ( x * t - t^3 / 3 ) dt / pi ! ! The data was reported by McLeod. ! ! Modified: ! ! 24 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.40936798278458884024D+00, & 0.37495291608048868619D+00, & 0.22066960679295989454D+00, & 0.77565356679703713590D-01, & 0.39638826473124717315D-01, & 0.38450072575004151871D-01, & 0.35273216868317898556D-01, & 0.31768535282502272742D-01, & 0.28894408288051391369D-01, & 0.24463284011678541180D-01, & 0.41053540139998941517D+00, & 0.44993502381204990817D+00, & 0.97220515514243332184D+00, & 0.83764237105104371193D+02, & 0.80327744952044756016D+05, & 0.15514138847749108298D+06, & 0.11995859641733262114D+07, & 0.21472868855967642259D+08, & 0.45564115351632913590D+09, & 0.32980722582904761929D+12 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & -0.0019531250D+00, & -0.1250000000D+00, & -1.0000000000D+00, & -4.0000000000D+00, & -8.0000000000D+00, & -8.2500000000D+00, & -9.0000000000D+00, & -10.0000000000D+00, & -11.0000000000D+00, & -13.0000000000D+00, & 0.0019531250D+00, & 0.1250000000D+00, & 1.0000000000D+00, & 4.0000000000D+00, & 7.0000000000D+00, & 7.2500000000D+00, & 8.0000000000D+00, & 9.0000000000D+00, & 10.0000000000D+00, & 12.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function arctan_int ( xvalue ) !*****************************************************************************80 ! !! ARCTAN_INT calculates the inverse tangent integral. ! ! Discussion: ! ! The function is defined by: ! ! ARCTAN_INT(x) = Integral ( 0 <= t <= x ) arctan ( t ) / t dt ! ! The approximation uses Chebyshev series with the coefficients ! given to an accuracy of 20D. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 24 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) ARCTAN_INT, the value of the function. ! implicit none real ( kind = 8 ), dimension ( 0:22 ) :: atnina = (/ & 1.91040361296235937512d0, & -0.4176351437656746940d-1, & 0.275392550786367434d-2, & -0.25051809526248881d-3, & 0.2666981285121171d-4, & -0.311890514107001d-5, & 0.38833853132249d-6, & -0.5057274584964d-7, & 0.681225282949d-8, & -0.94212561654d-9, & 0.13307878816d-9, & -0.1912678075d-10, & 0.278912620d-11, & -0.41174820d-12, & 0.6142987d-13, & -0.924929d-14, & 0.140387d-14, & -0.21460d-15, & 0.3301d-16, & -0.511d-17, & 0.79d-18, & -0.12d-18, & 0.2d-19 /) real ( kind = 8 ) arctan_int real ( kind = 8 ) cheval real ( kind = 8 ), parameter :: half = 0.5D+00 integer ind integer, parameter :: nterms = 19 real ( kind = 8 ), parameter :: one = 1.0D+00 real ( kind = 8 ) t real ( kind = 8 ), parameter :: twobpi = 0.63661977236758134308D+00 real ( kind = 8 ) x real ( kind = 8 ), parameter :: xlow = 7.4505806D-09 real ( kind = 8 ), parameter :: xupper = 4.5036D+15 real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 ind = 1 x = xvalue if ( x < zero ) then x = -x ind = -1 end if if ( x < xlow ) then arctan_int = x else if ( x <= one ) then t = x * x t = ( t - half ) + ( t - half ) arctan_int = x * cheval ( nterms, atnina, t ) else if ( x <= xupper ) then t = one / ( x * x ) t = ( t - half ) + ( t - half ) arctan_int = log ( x ) / twobpi + cheval ( nterms, atnina, t ) / x else arctan_int = log ( x ) / twobpi end if if ( ind < 0 ) then arctan_int = -arctan_int end if return end subroutine arctan_int_values ( n_data, x, fx ) !*****************************************************************************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. ! ! Modified: ! ! 25 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.19531241721588483191D-02, & -0.39062433772980711281D-02, & 0.78124470192576499535D-02, & 0.15624576181996527280D-01, & -0.31246610349485401551D-01, & 0.62472911335014397321D-01, & 0.12478419717389654039D+00, & -0.24830175098230686908D+00, & 0.48722235829452235711D+00, & 0.91596559417721901505D+00, & 0.12749694484943800618D+01, & -0.15760154034463234224D+01, & 0.24258878412859089996D+01, & 0.33911633326292997361D+01, & 0.44176450919422186583D+01, & -0.47556713749547247774D+01, & 0.50961912150934111303D+01, & 0.53759175735714876256D+01, & -0.61649904785027487422D+01, & 0.72437843013083534973D+01 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0019531250D+00, & -0.0039062500D+00, & 0.0078125000D+00, & 0.0156250000D+00, & -0.0312500000D+00, & 0.0625000000D+00, & 0.1250000000D+00, & -0.2500000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 1.5000000000D+00, & -2.0000000000D+00, & 4.0000000000D+00, & 8.0000000000D+00, & 16.0000000000D+00, & -20.0000000000D+00, & 25.0000000000D+00, & 30.0000000000D+00, & -50.0000000000D+00, & 100.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function bessel_i0_int ( xvalue ) !*****************************************************************************80 ! !! BESSEL_I0_INT computes the integral of the modified Bessel function I0(X). ! ! Discussion: ! ! The function is defined by: ! ! I0_INT(x) = Integral ( 0 <= t <= x ) I0(t) dt ! ! The program uses Chebyshev expansions, the coefficients of ! which are given to 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 29 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) BESSEL_I0_INT, the value of the function. ! implicit none real ( kind = 8 ) ari01(0:28) real ( kind = 8 ) ari0a(0:33) real ( kind = 8 ), parameter :: ateen = 18.0D+00 real ( kind = 8 ) bessel_i0_int real ( kind = 8 ) cheval real ( kind = 8 ), parameter :: half = 0.5D+00 integer ind real ( kind = 8 ), parameter :: lnr2pi = 0.91893853320467274178D+00 integer, parameter :: nterm1 = 25 integer, parameter :: nterm2 = 27 real ( kind = 8 ) t real ( kind = 8 ) temp real ( kind = 8 ), parameter :: thirt6 = 36.0D+00 real ( kind = 8 ), parameter :: three = 3.0D+00 real ( kind = 8 ) x real ( kind = 8 ), parameter :: xhigh = 713.758339D+00 real ( kind = 8 ), parameter :: xlow = 0.5161914D-07 real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 data ari01(0)/ 0.41227906926781516801d0/ data ari01(1)/ -0.34336345150081519562d0/ data ari01(2)/ 0.22667588715751242585d0/ data ari01(3)/ -0.12608164718742260032d0/ data ari01(4)/ 0.6012484628777990271d-1/ data ari01(5)/ -0.2480120462913358248d-1/ data ari01(6)/ 0.892773389565563897d-2/ data ari01(7)/ -0.283253729936696605d-2/ data ari01(8)/ 0.79891339041712994d-3/ data ari01(9)/ -0.20053933660964890d-3/ data ari01(10)/ 0.4416816783014313d-4/ data ari01(11)/-0.822377042246068d-5/ data ari01(12)/ 0.120059794219015d-5/ data ari01(13)/-0.11350865004889d-6/ data ari01(14)/ 0.69606014466d-9/ data ari01(15)/ 0.180622772836d-8/ data ari01(16)/-0.26039481370d-9/ data ari01(17)/-0.166188103d-11/ data ari01(18)/ 0.510500232d-11/ data ari01(19)/-0.41515879d-12/ data ari01(20)/-0.7368138d-13/ data ari01(21)/ 0.1279323d-13/ data ari01(22)/ 0.103247d-14/ data ari01(23)/-0.30379d-15/ data ari01(24)/-0.1789d-16/ data ari01(25)/ 0.673d-17/ data ari01(26)/ 0.44d-18/ data ari01(27)/-0.14d-18/ data ari01(28)/-0.1d-19/ data ari0a(0)/ 2.03739654571143287070d0/ data ari0a(1)/ 0.1917631647503310248d-1/ data ari0a(2)/ 0.49923334519288147d-3/ data ari0a(3)/ 0.2263187103659815d-4/ data ari0a(4)/ 0.158682108285561d-5/ data ari0a(5)/ 0.16507855636318d-6/ data ari0a(6)/ 0.2385058373640d-7/ data ari0a(7)/ 0.392985182304d-8/ data ari0a(8)/ 0.46042714199d-9/ data ari0a(9)/ -0.7072558172d-10/ data ari0a(10)/-0.6747183961d-10/ data ari0a(11)/-0.2026962001d-10/ data ari0a(12)/-0.87320338d-12/ data ari0a(13)/ 0.175520014d-11/ data ari0a(14)/ 0.60383944d-12/ data ari0a(15)/-0.3977983d-13/ data ari0a(16)/-0.8049048d-13/ data ari0a(17)/-0.1158955d-13/ data ari0a(18)/ 0.827318d-14/ data ari0a(19)/ 0.282290d-14/ data ari0a(20)/-0.77667d-15/ data ari0a(21)/-0.48731d-15/ data ari0a(22)/ 0.7279d-16/ data ari0a(23)/ 0.7873d-16/ data ari0a(24)/-0.785d-17/ data ari0a(25)/-0.1281d-16/ data ari0a(26)/ 0.121d-17/ data ari0a(27)/ 0.214d-17/ data ari0a(28)/-0.27d-18/ data ari0a(29)/-0.36d-18/ data ari0a(30)/ 0.7d-19/ data ari0a(31)/ 0.6d-19/ data ari0a(32)/-0.2d-19/ data ari0a(33)/-0.1d-19/ ind = 1 x = xvalue if ( xvalue < zero ) then ind = -1 x = -x end if if ( x < xlow ) then bessel_i0_int = x else if ( x <= ateen ) then t = ( three * x - ateen ) / ( x + ateen ) bessel_i0_int = x * exp ( x ) * cheval ( nterm1, ari01, t ) else if ( x <= xhigh ) then t = ( thirt6 / x - half ) - half temp = x - half * log ( x ) - lnr2pi + log ( cheval ( nterm2, ari0a, t )) bessel_i0_int = exp ( temp ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BESSEL_I0_INT - Fatal error!' write ( *, '(a)' ) ' Argument magnitude too large.' bessel_i0_int = exp ( xhigh - lnr2pi - half * log ( xhigh ) ) end if if ( ind == -1 ) then bessel_i0_int = -bessel_i0_int end if return end subroutine bessel_i0_int_values ( n_data, x, fx ) !*****************************************************************************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. ! ! Modified: ! ! 29 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.19531256208818052282D-02, & -0.39062549670565734544D-02, & 0.62520348032546565850D-01, & 0.12516285581366971819D+00, & -0.51051480879740303760D+00, & 0.10865210970235898158D+01, & 0.27750019054282535299D+01, & -0.13775208868039716639D+02, & 0.46424372058106108576D+03, & 0.64111867658021584522D+07, & -0.10414860803175857953D+08, & 0.44758598913855743089D+08, & -0.11852985311558287888D+09, & 0.31430078220715992752D+09, & -0.83440212900794309620D+09, & 0.22175367579074298261D+10, & 0.58991731842803636487D+10, & -0.41857073244691522147D+11, & 0.79553885818472357663D+12, & 0.15089715082719201025D+17 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0019531250D+00, & -0.0039062500D+00, & 0.0625000000D+00, & 0.1250000000D+00, & -0.5000000000D+00, & 1.0000000000D+00, & 2.0000000000D+00, & -4.0000000000D+00, & 8.0000000000D+00, & 18.0000000000D+00, & -18.5000000000D+00, & 20.0000000000D+00, & -21.0000000000D+00, & 22.0000000000D+00, & -23.0000000000D+00, & 24.0000000000D+00, & 25.0000000000D+00, & -27.0000000000D+00, & 30.0000000000D+00, & 40.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function bessel_j0_int ( xvalue ) !*****************************************************************************80 ! !! BESSEL_J0_INT calculates the integral of the Bessel function J0. ! ! Discussion: ! ! The function is defined by: ! ! J0_INT(x) = Integral ( 0 <= t <= x ) J0(t) dt ! ! The code uses Chebyshev expansions whose coefficients are ! given to 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 07 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) BESSEL_J0_INT, the value of the function. ! implicit none real ( kind = 8 ) bessel_j0_int real ( kind = 8 ) cheval integer ind integer, parameter :: nterm1 = 22 integer, parameter :: nterm2 = 18 integer, parameter :: nterm3 = 16 real ( kind = 8 ), parameter :: one = 1.0D+00 real ( kind = 8 ) x real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 real ( kind = 8 ) arj01(0:23),arj0a1(0:21),arj0a2(0:18), & five12,one28,pib41,pib411,pib412, & pib42,rt2bpi,sixten,t,temp,xhigh,xlow, & xmpi4 data sixten/ 16.0d0 / data one28,five12/ 128.0d0 , 512d0 / data rt2bpi/0.79788456080286535588d0/ data pib411,pib412/ 201.0d0 , 256.0d0/ data pib42/0.24191339744830961566d-3/ data arj01(0)/ 0.38179279321690173518d0/ data arj01(1)/ -0.21275636350505321870d0/ data arj01(2)/ 0.16754213407215794187d0/ data arj01(3)/ -0.12853209772196398954d0/ data arj01(4)/ 0.10114405455778847013d0/ data arj01(5)/ -0.9100795343201568859d-1/ data arj01(6)/ 0.6401345264656873103d-1/ data arj01(7)/ -0.3066963029926754312d-1/ data arj01(8)/ 0.1030836525325064201d-1/ data arj01(9)/ -0.255670650399956918d-2/ data arj01(10)/ 0.48832755805798304d-3/ data arj01(11)/-0.7424935126036077d-4/ data arj01(12)/ 0.922260563730861d-5/ data arj01(13)/-0.95522828307083d-6/ data arj01(14)/ 0.8388355845986d-7/ data arj01(15)/-0.633184488858d-8/ data arj01(16)/ 0.41560504221d-9/ data arj01(17)/-0.2395529307d-10/ data arj01(18)/ 0.122286885d-11/ data arj01(19)/-0.5569711d-13/ data arj01(20)/ 0.227820d-14/ data arj01(21)/-0.8417d-16/ data arj01(22)/ 0.282d-17/ data arj01(23)/-0.9d-19/ data arj0a1(0)/ 1.24030133037518970827d0/ data arj0a1(1)/ -0.478125353632280693d-2/ data arj0a1(2)/ 0.6613148891706678d-4/ data arj0a1(3)/ -0.186042740486349d-5/ data arj0a1(4)/ 0.8362735565080d-7/ data arj0a1(5)/ -0.525857036731d-8/ data arj0a1(6)/ 0.42606363251d-9/ data arj0a1(7)/ -0.4211761024d-10/ data arj0a1(8)/ 0.488946426d-11/ data arj0a1(9)/ -0.64834929d-12/ data arj0a1(10)/ 0.9617234d-13/ data arj0a1(11)/-0.1570367d-13/ data arj0a1(12)/ 0.278712d-14/ data arj0a1(13)/-0.53222d-15/ data arj0a1(14)/ 0.10844d-15/ data arj0a1(15)/-0.2342d-16/ data arj0a1(16)/ 0.533d-17/ data arj0a1(17)/-0.127d-17/ data arj0a1(18)/ 0.32d-18/ data arj0a1(19)/-0.8d-19/ data arj0a1(20)/ 0.2d-19/ data arj0a1(21)/-0.1d-19/ data arj0a2(0)/ 1.99616096301341675339d0/ data arj0a2(1)/ -0.190379819246668161d-2/ data arj0a2(2)/ 0.1539710927044226d-4/ data arj0a2(3)/ -0.31145088328103d-6/ data arj0a2(4)/ 0.1110850971321d-7/ data arj0a2(5)/ -0.58666787123d-9/ data arj0a2(6)/ 0.4139926949d-10/ data arj0a2(7)/ -0.365398763d-11/ data arj0a2(8)/ 0.38557568d-12/ data arj0a2(9)/ -0.4709800d-13/ data arj0a2(10)/ 0.650220d-14/ data arj0a2(11)/-0.99624d-15/ data arj0a2(12)/ 0.16700d-15/ data arj0a2(13)/-0.3028d-16/ data arj0a2(14)/ 0.589d-17/ data arj0a2(15)/-0.122d-17/ data arj0a2(16)/ 0.27d-18/ data arj0a2(17)/-0.6d-19/ data arj0a2(18)/ 0.1d-19/ ! ! Machine-dependent constants (suitable for IEEE machines) ! data xlow,xhigh/3.650024d-8,9.0072d15/ x = xvalue ind = 1 if ( x < zero ) then x = -x ind = -1 end if if ( x < xlow ) then bessel_j0_int = x else if ( x <= sixten ) then t = x * x / one28 - one bessel_j0_int = x * cheval ( nterm1, arj01, t ) else if ( x <= xhigh ) then t = five12 / ( x * x ) - one pib41 = pib411 / pib412 xmpi4 = ( x - pib41 ) - pib42 temp = cos ( xmpi4 ) * cheval ( nterm2, arj0a1, t ) / x temp = temp - sin ( xmpi4) * cheval ( nterm3, arj0a2, t ) bessel_j0_int = one - rt2bpi * temp / sqrt ( x ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BESSEL_J0_INT - Fatal error!' write ( *, '(a)' ) ' Argument magnitude too large.' bessel_j0_int = one end if if ( ind == -1 ) then bessel_j0_int = -bessel_j0_int end if return end subroutine bessel_j0_int_values ( n_data, x, fx ) !*****************************************************************************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. ! ! Modified: ! ! 29 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.97656242238978822427D-03, & 0.39062450329491108875D-02, & -0.62479657927917933620D-01, & 0.12483733492120479139D+00, & -0.48968050664604505505D+00, & 0.91973041008976023931D+00, & -0.14257702931970265690D+01, & 0.10247341594606064818D+01, & -0.12107468348304501655D+01, & 0.11008652032736190799D+01, & -0.10060334829904124192D+01, & 0.81330572662485953519D+00, & -0.10583788214211277585D+01, & 0.87101492116545875169D+00, & -0.88424908882547488420D+00, & 0.11257761503599914603D+01, & -0.90141212258183461184D+00, & 0.91441344369647797803D+00, & -0.94482281938334394886D+00, & 0.92266255696016607257D+00 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0009765625D+00, & 0.0039062500D+00, & -0.0625000000D+00, & 0.1250000000D+00, & -0.5000000000D+00, & 1.0000000000D+00, & -2.0000000000D+00, & 4.0000000000D+00, & -8.0000000000D+00, & 16.0000000000D+00, & -16.5000000000D+00, & 18.0000000000D+00, & -20.0000000000D+00, & 25.0000000000D+00, & -30.0000000000D+00, & 40.0000000000D+00, & -50.0000000000D+00, & 75.0000000000D+00, & -80.0000000000D+00, & 100.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function bessel_k0_int ( xvalue ) !*****************************************************************************80 ! !! BESSEL_K0_INT calculates the integral of the modified Bessel function K0(X). ! ! Discussion: ! ! The function is defined by: ! ! K0_INT(x) = Integral ( 0 <= t <= x ) K0(t) dt ! ! The code uses Chebyshev expansions, whose coefficients are ! given to 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 29 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) BESSEL_K0_INT, the value of the function. ! implicit none real ( kind = 8 ) bessel_k0_int real ( kind = 8 ) cheval real ( kind = 8 ), parameter :: half = 0.5D+00 integer, parameter :: nterm1 = 14 integer, parameter :: nterm2 = 14 integer, parameter :: nterm3 = 23 real ( kind = 8 ), parameter :: six = 6.0D+00 real ( kind = 8 ) x real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 real ( kind = 8 ) ak0in1(0:15),ak0in2(0:15),ak0ina(0:27), & const1,const2,eightn,fval, & piby2,rt2bpi,t,temp,twelve, & xhigh,xlow data twelve,eightn / 12.0d0 , 18.0d0 / data const1/1.11593151565841244881d0/ data const2/-0.11593151565841244881d0/ data piby2/1.57079632679489661923d0/ data rt2bpi/0.79788456080286535588d0/ data ak0in1/16.79702714464710959477d0, & 9.79134687676889407070d0, & 2.80501316044337939300d0, & 0.45615620531888502068d0, & 0.4716224457074760784d-1, & 0.335265148269698289d-2, & 0.17335181193874727d-3, & 0.679951889364702d-5, & 0.20900268359924d-6, & 0.516603846976d-8, & 0.10485708331d-9, & 0.177829320d-11, & 0.2556844d-13, & 0.31557d-15, & 0.338d-17, & 0.3d-19/ data ak0in2/10.76266558227809174077d0, & 5.62333479849997511550d0, & 1.43543664879290867158d0, & 0.21250410143743896043d0, & 0.2036537393100009554d-1, & 0.136023584095623632d-2, & 0.6675388699209093d-4, & 0.250430035707337d-5, & 0.7406423741728d-7, & 0.176974704314d-8, & 0.3485775254d-10, & 0.57544785d-12, & 0.807481d-14, & 0.9747d-16, & 0.102d-17, & 0.1d-19/ data ak0ina(0)/ 1.91172065445060453895d0/ data ak0ina(1)/ -0.4183064565769581085d-1/ data ak0ina(2)/ 0.213352508068147486d-2/ data ak0ina(3)/ -0.15859497284504181d-3/ data ak0ina(4)/ 0.1497624699858351d-4/ data ak0ina(5)/ -0.167955955322241d-5/ data ak0ina(6)/ 0.21495472478804d-6/ data ak0ina(7)/ -0.3058356654790d-7/ data ak0ina(8)/ 0.474946413343d-8/ data ak0ina(9)/ -0.79424660432d-9/ data ak0ina(10)/ 0.14156555325d-9/ data ak0ina(11)/-0.2667825359d-10/ data ak0ina(12)/ 0.528149717d-11/ data ak0ina(13)/-0.109263199d-11/ data ak0ina(14)/ 0.23518838d-12/ data ak0ina(15)/-0.5247991d-13/ data ak0ina(16)/ 0.1210191d-13/ data ak0ina(17)/-0.287632d-14/ data ak0ina(18)/ 0.70297d-15/ data ak0ina(19)/-0.17631d-15/ data ak0ina(20)/ 0.4530d-16/ data ak0ina(21)/-0.1190d-16/ data ak0ina(22)/ 0.319d-17/ data ak0ina(23)/-0.87d-18/ data ak0ina(24)/ 0.24d-18/ data ak0ina(25)/-0.7d-19/ data ak0ina(26)/ 0.2d-19/ data ak0ina(27)/-0.1d-19/ ! ! Machine-dependent values (suitable for IEEE machines) ! data xlow,xhigh/4.47034836d-8,36.0436534d0/ x = xvalue if ( x < zero ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BESSEL_K0_INT - Fatal error!' write ( *, '(a)' ) ' Argument X < 0.' bessel_k0_int = zero else if ( x == zero ) then bessel_k0_int = zero else if ( x < xlow ) then bessel_k0_int = x * ( const1 - log ( x ) ) else if ( x <= six ) then t = ( ( x * x ) / eightn - half ) - half fval = ( const2 + log ( x ) ) * cheval ( nterm2, ak0in2, t ) bessel_k0_int = x * ( cheval ( nterm1, ak0in1, t ) - fval ) else if ( x < xhigh ) then fval = piby2 t = ( twelve / x - half ) - half temp = exp ( -x ) * cheval ( nterm3, ak0ina, t ) fval = fval - temp / ( sqrt ( x ) * rt2bpi ) bessel_k0_int = fval else bessel_k0_int = piby2 end if return end subroutine bessel_k0_int_values ( n_data, x, fx ) !*****************************************************************************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. ! ! Modified: ! ! 29 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.78587929563466784589D-02, & 0.26019991617330578111D-01, & 0.24311842237541167904D+00, & 0.39999633750480508861D+00, & 0.92710252093114907345D+00, & 0.12425098486237782662D+01, & 0.14736757343168286825D+01, & 0.15606495706051741364D+01, & 0.15673873907283660493D+01, & 0.15696345532693743714D+01, & 0.15701153443250786355D+01, & 0.15706574852894436220D+01, & 0.15707793116159788598D+01, & 0.15707942066465767196D+01, & 0.15707962315469192247D+01, & 0.15707963262340149876D+01, & 0.15707963267948756308D+01, & 0.15707963267948966192D+01, & 0.15707963267948966192D+01, & 0.15707963267948966192D+01 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0009765625D+00, & 0.0039062500D+00, & 0.0625000000D+00, & 0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 2.0000000000D+00, & 4.0000000000D+00, & 5.0000000000D+00, & 6.0000000000D+00, & 6.5000000000D+00, & 8.0000000000D+00, & 10.0000000000D+00, & 12.0000000000D+00, & 15.0000000000D+00, & 20.0000000000D+00, & 30.0000000000D+00, & 50.0000000000D+00, & 80.0000000000D+00, & 100.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function bessel_y0_int ( xvalue ) !*****************************************************************************80 ! !! BESSEL_Y0_INT calculates the integral of the Bessel function Y0. ! ! Discussion: ! ! The function is defined by: ! ! Y0_INT(x) = Integral ( 0 <= t <= x ) Y0(t) dt ! ! The code uses Chebyshev expansions whose coefficients are ! given to 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 23 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) BESSEL_Y0_INT, the value of the function. ! implicit none real ( kind = 8 ) bessel_y0_int real ( kind = 8 ) cheval real ( kind = 8 ), parameter :: nine = 9.0D+00 integer, parameter :: nterm1 = 22 integer, parameter :: nterm2 = 22 integer, parameter :: nterm3 = 17 integer, parameter :: nterm4 = 15 real ( kind = 8 ), parameter :: one = 1.0D+00 real ( kind = 8 ), parameter :: sixten = 16.0D+00 real ( kind = 8 ) x real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 real ( kind = 8 ) arj01(0:23),ary01(0:24),ary0a1(0:21), & ary0a2(0:18),five12,gal2m1,gamln2, & one28,pib41,pib411,pib412, & pib42,rt2bpi,t,temp,twobpi,xhigh, & xlow,xmpi4 data one28,five12/ 128.0d0 , 512.0d0 / data rt2bpi/0.79788456080286535588d0/ data pib411,pib412/ 201.0d0, 256.0d0/ data pib42/0.24191339744830961566d-3/ data twobpi/0.63661977236758134308d0/ data gal2m1/-1.11593151565841244881d0/ data gamln2/-0.11593151565841244881d0/ data arj01(0)/ 0.38179279321690173518d0/ data arj01(1)/ -0.21275636350505321870d0/ data arj01(2)/ 0.16754213407215794187d0/ data arj01(3)/ -0.12853209772196398954d0/ data arj01(4)/ 0.10114405455778847013d0/ data arj01(5)/ -0.9100795343201568859d-1/ data arj01(6)/ 0.6401345264656873103d-1/ data arj01(7)/ -0.3066963029926754312d-1/ data arj01(8)/ 0.1030836525325064201d-1/ data arj01(9)/ -0.255670650399956918d-2/ data arj01(10)/ 0.48832755805798304d-3/ data arj01(11)/-0.7424935126036077d-4/ data arj01(12)/ 0.922260563730861d-5/ data arj01(13)/-0.95522828307083d-6/ data arj01(14)/ 0.8388355845986d-7/ data arj01(15)/-0.633184488858d-8/ data arj01(16)/ 0.41560504221d-9/ data arj01(17)/-0.2395529307d-10/ data arj01(18)/ 0.122286885d-11/ data arj01(19)/-0.5569711d-13/ data arj01(20)/ 0.227820d-14/ data arj01(21)/-0.8417d-16/ data arj01(22)/ 0.282d-17/ data arj01(23)/-0.9d-19/ data ary01(0)/ 0.54492696302724365490d0/ data ary01(1)/ -0.14957323588684782157d0/ data ary01(2)/ 0.11085634486254842337d0/ data ary01(3)/ -0.9495330018683777109d-1/ data ary01(4)/ 0.6820817786991456963d-1/ data ary01(5)/ -0.10324653383368200408d0/ data ary01(6)/ 0.10625703287534425491d0/ data ary01(7)/ -0.6258367679961681990d-1/ data ary01(8)/ 0.2385645760338293285d-1/ data ary01(9)/ -0.644864913015404481d-2/ data ary01(10)/ 0.131287082891002331d-2/ data ary01(11)/-0.20988088174989640d-3/ data ary01(12)/ 0.2716042484138347d-4/ data ary01(13)/-0.291199114014694d-5/ data ary01(14)/ 0.26344333093795d-6/ data ary01(15)/-0.2041172069780d-7/ data ary01(16)/ 0.137124781317d-8/ data ary01(17)/-0.8070680792d-10/ data ary01(18)/ 0.419883057d-11/ data ary01(19)/-0.19459104d-12/ data ary01(20)/ 0.808782d-14/ data ary01(21)/-0.30329d-15/ data ary01(22)/ 0.1032d-16/ data ary01(23)/-0.32d-18/ data ary01(24)/ 0.1d-19/ data ary0a1(0)/ 1.24030133037518970827d0/ data ary0a1(1)/ -0.478125353632280693d-2/ data ary0a1(2)/ 0.6613148891706678d-4/ data ary0a1(3)/ -0.186042740486349d-5/ data ary0a1(4)/ 0.8362735565080d-7/ data ary0a1(5)/ -0.525857036731d-8/ data ary0a1(6)/ 0.42606363251d-9/ data ary0a1(7)/ -0.4211761024d-10/ data ary0a1(8)/ 0.488946426d-11/ data ary0a1(9)/ -0.64834929d-12/ data ary0a1(10)/ 0.9617234d-13/ data ary0a1(11)/-0.1570367d-13/ data ary0a1(12)/ 0.278712d-14/ data ary0a1(13)/-0.53222d-15/ data ary0a1(14)/ 0.10844d-15/ data ary0a1(15)/-0.2342d-16/ data ary0a1(16)/ 0.533d-17/ data ary0a1(17)/-0.127d-17/ data ary0a1(18)/ 0.32d-18/ data ary0a1(19)/-0.8d-19/ data ary0a1(20)/ 0.2d-19/ data ary0a1(21)/-0.1d-19/ data ary0a2(0)/ 1.99616096301341675339d0/ data ary0a2(1)/ -0.190379819246668161d-2/ data ary0a2(2)/ 0.1539710927044226d-4/ data ary0a2(3)/ -0.31145088328103d-6/ data ary0a2(4)/ 0.1110850971321d-7/ data ary0a2(5)/ -0.58666787123d-9/ data ary0a2(6)/ 0.4139926949d-10/ data ary0a2(7)/ -0.365398763d-11/ data ary0a2(8)/ 0.38557568d-12/ data ary0a2(9)/ -0.4709800d-13/ data ary0a2(10)/ 0.650220d-14/ data ary0a2(11)/-0.99624d-15/ data ary0a2(12)/ 0.16700d-15/ data ary0a2(13)/-0.3028d-16/ data ary0a2(14)/ 0.589d-17/ data ary0a2(15)/-0.122d-17/ data ary0a2(16)/ 0.27d-18/ data ary0a2(17)/-0.6d-19/ data ary0a2(18)/ 0.1d-19/ ! ! Machine-dependent constants (suitable for IEEE machines) ! data xlow,xhigh/3.16101364d-8,9.007199256d15/ x = xvalue if ( x < zero ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BESSEL_Y0_INT - Fatal error!' write ( *, '(a)' ) ' Argument X < 0.' bessel_y0_int = zero else if ( x == zero ) then bessel_y0_int = zero else if ( x < xlow ) then bessel_y0_int = ( log ( x ) + gal2m1 ) * twobpi * x else if ( x <= sixten ) then t = x * x / one28 - one temp = ( log ( x ) + gamln2 ) * cheval ( nterm1, arj01, t ) temp = temp - cheval ( nterm2, ary01, t ) bessel_y0_int = twobpi * x * temp else if ( x <= xhigh ) then t = five12 / ( x * x ) - one pib41 = pib411 / pib412 xmpi4 = ( x - pib41 ) - pib42 temp = sin ( xmpi4 ) * cheval ( nterm3, ary0a1, t ) / x temp = temp + cos ( xmpi4 ) * cheval ( nterm4, ary0a2, t ) bessel_y0_int = - rt2bpi * temp / sqrt ( x ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BESSEL_Y0_INT - Fatal error!' write ( *, '(a)' ) ' Argument too large.' bessel_y0_int = zero end if return end subroutine bessel_y0_int_values ( n_data, x, fx ) !*****************************************************************************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. ! ! Modified: ! ! 30 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & -0.91442642860172110926D-02, & -0.29682047390397591290D-01, & -0.25391431276585388961D+00, & -0.56179545591464028187D+00, & -0.63706937660742309754D+00, & -0.28219285008510084123D+00, & 0.38366964785312561103D+00, & -0.12595061285798929390D+00, & 0.24129031832266684828D+00, & 0.17138069757627037938D+00, & 0.18958142627134083732D+00, & 0.17203846136449706946D+00, & -0.16821597677215029611D+00, & -0.93607927351428988679D-01, & 0.88229711948036648408D-01, & -0.89324662736274161841D-02, & -0.54814071000063488284D-01, & -0.94958246003466381588D-01, & -0.19598064853404969850D-01, & -0.83084772357154773468D-02 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0019531250D+00, & 0.0078125000D+00, & 0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 2.0000000000D+00, & 4.0000000000D+00, & 6.0000000000D+00, & 10.0000000000D+00, & 16.0000000000D+00, & 16.2500000000D+00, & 17.0000000000D+00, & 20.0000000000D+00, & 25.0000000000D+00, & 30.0000000000D+00, & 40.0000000000D+00, & 50.0000000000D+00, & 70.0000000000D+00, & 100.0000000000D+00, & 125.0000000000D+00 /) ! if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function cheval ( n, a, t ) !*****************************************************************************80 ! !! CHEVAL evaluates a Chebyshev series. ! ! Discussion: ! ! This function evaluates a Chebyshev series, using the ! Clenshaw method with Reinsch modification, as analysed ! in the paper by Oliver. ! ! Modified: ! ! 07 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! J Oliver, ! An error analysis of the modified Clenshaw method for ! evaluating Chebyshev and Fourier series, ! Journal of the IMA, ! Volume 20, 1977, pages 379-391. ! ! Parameters: ! ! Input, integer N, the number of terms in the sequence. ! ! Input, real ( kind = 8 ) A(0:N), the coefficients of the Chebyshev series. ! ! Input, real ( kind = 8 ) T, the value at which the series is ! to be evaluated. ! ! Output, real ( kind = 8 ) CHEVAL, the value of the Chebyshev series at T. ! implicit none integer n real ( kind = 8 ) :: a(0:n) real ( kind = 8 ) :: cheval real ( kind = 8 ) :: d1 real ( kind = 8 ) :: d2 real ( kind = 8 ), parameter :: half = 0.5D+00 integer i real ( kind = 8 ) :: t real ( kind = 8 ), parameter :: test = 0.6D+00 real ( kind = 8 ) :: tt real ( kind = 8 ), parameter :: two = 2.0D+00 real ( kind = 8 ) :: u0 real ( kind = 8 ) :: u1 real ( kind = 8 ) :: u2 real ( kind = 8 ), parameter :: zero = 0.0D+00 u1 = zero ! ! T <= -0.6, Reinsch modification. ! if ( t <= -test ) then d1 = zero tt = ( t + half ) + half tt = tt + tt do i = n, 0, -1 d2 = d1 u2 = u1 d1 = tt * u2 + a(i) - d2 u1 = d1 - u2 end do cheval = ( d1 - d2 ) / two ! ! -0.6 < T < 0.6, Standard Clenshaw method. ! else if ( t < test ) then u0 = zero tt = t + t do i = n, 0, -1 u2 = u1 u1 = u0 u0 = tt * u1 + a(i) - u2 end do cheval = ( u0 - u2 ) / two ! ! 0.6 <= T, Reinsch modification. ! else d1 = zero tt = ( t - half ) - half tt = tt + tt do i = n, 0, -1 d2 = d1 u2 = u1 d1 = tt * u2 + a(i) + d2 u1 = d1 + u2 end do cheval = ( d1 + d2 ) / two end if return end function clausen ( xvalue ) !*****************************************************************************80 ! !! CLAUSEN calculates Clausen's integral. ! ! Discussion: ! ! The function is defined by: ! ! CLAUSEN(x) = Integral ( 0 <= t <= x ) -ln ( 2 * sin ( t / 2 ) ) dt ! ! The code uses Chebyshev expansions with the coefficients ! given to 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 07 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) CLAUSEN, the value of the function. ! implicit none real ( kind = 8 ) aclaus(0:15) real ( kind = 8 ) cheval real ( kind = 8 ) clausen real ( kind = 8 ), parameter :: half = 0.5D+00 integer indx integer, parameter :: nterms = 13 real ( kind = 8 ), parameter :: one = 1.0D+00 real ( kind = 8 ), parameter :: pi = 3.1415926535897932385D+00 real ( kind = 8 ), parameter :: pisq = 9.8696044010893586188D+00 real ( kind = 8 ) t real ( kind = 8 ) x real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 real ( kind = 8 ) twopi,twopia,twopib,xhigh,xsmall data twopi/6.2831853071795864769d0/ data twopia,twopib/6.28125d0, 0.19353071795864769253d-2/ data aclaus/2.14269436376668844709d0, & 0.7233242812212579245d-1, & 0.101642475021151164d-2, & 0.3245250328531645d-4, & 0.133315187571472d-5, & 0.6213240591653d-7, & 0.313004135337d-8, & 0.16635723056d-9, & 0.919659293d-11, & 0.52400462d-12, & 0.3058040d-13, & 0.181969d-14, & 0.11004d-15, & 0.675d-17, & 0.42d-18, & 0.3d-19/ ! ! Set machine-dependent constants (suitable for IEEE machines) ! data xsmall,xhigh/2.3406689d-8,4.5036d15/ x = xvalue if ( xhigh < abs ( x ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CLAUSEN - Warning!' write ( *, '(a)' ) ' Argument magnitude too large for accurate computation.' clausen = zero return end if indx = 1 if ( x < zero ) then x = -x indx = -1 end if ! ! Argument reduced using simulated extra precision ! if ( twopi < x ) then t = aint ( x / twopi ) x = ( x - t * twopia ) - t * twopib end if if ( pi < x ) then x = ( twopia - x ) + twopib indx = -indx end if if ( x == zero ) then clausen = zero else if ( x < xsmall ) then clausen = x * ( one - log ( x ) ) else t = ( x * x ) / pisq - half t = t + t if ( one < t ) then t = one end if clausen = x * cheval ( nterms, aclaus, t ) - x * log ( x ) end if if ( indx < 0 ) then clausen = -clausen end if return end subroutine clausen_values ( n_data, x, fx ) !*****************************************************************************80 ! !! CLAUSEN_VALUES returns some values of the Clausen's integral. ! ! Discussion: ! ! The function is defined by: ! ! CLAUSEN(x) = Integral ( 0 <= t <= x ) -ln ( 2 * sin ( t / 2 ) ) dt ! ! The data was reported by McLeod. ! ! Modified: ! ! 25 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.14137352886760576684D-01, & 0.13955467081981281934D+00, & -0.38495732156574238507D+00, & 0.84831187770367927099D+00, & 0.10139591323607685043D+01, & -0.93921859275409211003D+00, & 0.72714605086327924743D+00, & 0.43359820323553277936D+00, & -0.98026209391301421161D-01, & -0.56814394442986978080D+00, & -0.70969701784448921625D+00, & 0.99282013254695671871D+00, & -0.98127747477447367875D+00, & -0.64078266570172320959D+00, & 0.86027963733231192456D+00, & 0.39071647608680211043D+00, & 0.47574793926539191502D+00, & 0.10105014481412878253D+01, & 0.96332089044363075154D+00, & -0.61782699481929311757D+00 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0019531250D+00, & 0.0312500000D+00, & -0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & -1.5000000000D+00, & 2.0000000000D+00, & 2.5000000000D+00, & -3.0000000000D+00, & 4.0000000000D+00, & 4.2500000000D+00, & -5.0000000000D+00, & 5.5000000000D+00, & 6.0000000000D+00, & 8.0000000000D+00, & -10.0000000000D+00, & 15.0000000000D+00, & 20.0000000000D+00, & -30.0000000000D+00, & 50.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function debye1 ( xvalue ) !*****************************************************************************80 ! !! DEBYE1 calculates the Debye function of order 1. ! ! Discussion: ! ! The function is defined by: ! ! DEBYE1(x) = 1 / x * Integral ( 0 <= t <= x ) t / ( exp ( t ) - 1 ) dt ! ! The code uses Chebyshev series whose coefficients ! are given to 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 07 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) DEBYE1, the value of the function. ! implicit none real ( kind = 8 ) adeb1(0:18) real ( kind = 8 ) cheval real ( kind = 8 ) debye1 real ( kind = 8 ), parameter :: eight = 8.0D+00 real ( kind = 8 ), parameter :: four = 4.0D+00 real ( kind = 8 ), parameter :: half = 0.5D+00 integer i integer nexp integer, parameter :: nterms = 15 real ( kind = 8 ), parameter :: one = 1.0D+00 real ( kind = 8 ), parameter :: quart = 0.25D+00 real ( kind = 8 ) x real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 real ( kind = 8 ) debinf,expmx, & nine,rk,sum1,t,thirt6,xk,xlim,xlow, & xupper data nine,thirt6 /9.0d0, 36.0d0 / data debinf/0.60792710185402662866d0/ data adeb1/2.40065971903814101941d0, & 0.19372130421893600885d0, & -0.623291245548957703d-2, & 0.35111747702064800d-3, & -0.2282224667012310d-4, & 0.158054678750300d-5, & -0.11353781970719d-6, & 0.835833611875d-8, & -0.62644247872d-9, & 0.4760334890d-10, & -0.365741540d-11, & 0.28354310d-12, & -0.2214729d-13, & 0.174092d-14, & -0.13759d-15, & 0.1093d-16, & -0.87d-18, & 0.7d-19, & -0.1d-19/ ! ! Machine-dependent constants (suitable for IEEE machines) ! data xlow,xupper,xlim/0.298023d-7,35.35051d0,708.39642d0/ x = xvalue if ( x < zero ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DEBYE1 - Fatal error!' write ( *, '(a)' ) ' Argument X < 0.' debye1 = zero else if ( x < xlow ) then debye1 = ( ( x - nine ) * x + thirt6 ) / thirt6 else if ( x <= four ) then t = ( ( x * x / eight ) - half ) - half debye1 = cheval ( nterms, adeb1, t ) - quart * x else debye1 = one / ( x * debinf ) if ( x < xlim ) then expmx = exp ( -x ) if ( xupper < x ) then debye1 = debye1 - expmx * ( one + one / x ) else sum1 = zero rk = aint ( xlim / x ) nexp = int ( rk ) xk = rk * x do i = nexp, 1, -1 t = ( one + one / xk ) / rk sum1 = sum1 * expmx + t rk = rk - one xk = xk - x end do debye1 = debye1 - sum1 * expmx end if end if end if return end subroutine debye1_values ( n_data, x, fx ) !*****************************************************************************80 ! !! DEBYE1_VALUES returns some values of Debye's function of order 1. ! ! Discussion: ! ! The function is defined by: ! ! DEBYE1(x) = 1 / x * Integral ( 0 <= t <= x ) t / ( exp ( t ) - 1 ) dt ! ! The data was reported by McLeod. ! ! Modified: ! ! 27 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.99951182471380889183D+00, & 0.99221462647120597836D+00, & 0.96918395997895308324D+00, & 0.88192715679060552968D+00, & 0.77750463411224827642D+00, & 0.68614531078940204342D+00, & 0.60694728460981007205D+00, & 0.53878956907785587703D+00, & 0.48043521957304283829D+00, & 0.38814802129793784501D+00, & 0.36930802829242526815D+00, & 0.32087619770014612104D+00, & 0.29423996623154246701D+00, & 0.27126046678502189985D+00, & 0.20523930310221503723D+00, & 0.16444346567994602563D+00, & 0.10966194482735821276D+00, & 0.82246701178200016086D-01, & 0.54831135561510852445D-01, & 0.32898681336964528729D-01 /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0019531250D+00, & 0.0312500000D+00, & 0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & 1.5000000000D+00, & 2.0000000000D+00, & 2.5000000000D+00, & 3.0000000000D+00, & 4.0000000000D+00, & 4.2500000000D+00, & 5.0000000000D+00, & 5.5000000000D+00, & 6.0000000000D+00, & 8.0000000000D+00, & 10.0000000000D+00, & 15.0000000000D+00, & 20.0000000000D+00, & 30.0000000000D+00, & 50.0000000000D+00 /) ! if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function debye2 ( xvalue ) !*****************************************************************************80 ! !! DEBYE2 calculates the Debye function of order 2. ! ! Discussion: ! ! The function is defined by: ! ! DEBYE2(x) = 2 / x^2 * Integral ( 0 <= t <= x ) t^2 / ( exp ( t ) - 1 ) dt ! ! The code uses Chebyshev series whose coefficients ! are given to 20 decimal places. ! ! This subroutine is set up to work on IEEE machines. ! ! Modified: ! ! 24 August 2004 ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input, real ( kind = 8 ) XVALUE, the argument of the function. ! ! Output, real ( kind = 8 ) DEBYE2, the value of the function. ! implicit none real ( kind = 8 ) cheval real ( kind = 8 ) debye2 real ( kind = 8 ), parameter :: eight = 8.0D+00 real ( kind = 8 ), parameter :: four = 4.0D+00 real ( kind = 8 ), parameter :: half = 0.5D+00 integer i integer nexp integer, parameter :: nterms = 17 real ( kind = 8 ), parameter :: one = 1.0D+00 real ( kind = 8 ), parameter :: three = 3.0D+00 real ( kind = 8 ), parameter :: two = 2.0D+00 real ( kind = 8 ) x real ( kind = 8 ) xvalue real ( kind = 8 ), parameter :: zero = 0.0D+00 real ( kind = 8 ) adeb2(0:18),debinf,expmx, & rk,sum1,t,twent4,xk,xlim1, & xlim2,xlow,xupper data twent4/24.0d0/ data debinf/4.80822761263837714160d0/ data adeb2/2.59438102325707702826d0, & 0.28633572045307198337d0, & -0.1020626561580467129d-1, & 0.60491097753468435d-3, & -0.4052576589502104d-4, & 0.286338263288107d-5, & -0.20863943030651d-6, & 0.1552378758264d-7, & -0.117312800866d-8, & 0.8973585888d-10, & -0.693176137d-11, & 0.53980568d-12, & -0.4232405d-13, & 0.333778d-14, & -0.26455d-15, & 0.2106d-16, & -0.168d-17, & 0.13d-18, & -0.1d-19/ ! ! Machine-dependent constants ! data xlow,xupper/0.298023d-7,35.35051d0/ data xlim1,xlim2/708.39642d0,2.1572317d154/ ! x = xvalue if ( x < zero ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DEBYE2 - Fatal error!' write ( *, '(a)' ) ' Argument X < 0.' debye2 = zero else if ( x < xlow ) then debye2 = ( ( x - eight ) * x + twent4 ) / twent4 else if ( x <= four ) then t = ( ( x * x / eight ) - half ) - half debye2 = cheval ( nterms, adeb2, t ) - x / three else if ( x <= xupper ) then expmx = exp ( -x ) sum1 = zero rk = aint ( xlim1 / x ) nexp = int ( rk ) xk = rk * x do i = nexp, 1, -1 t = ( one + two / xk + two / ( xk * xk ) ) / rk sum1 = sum1 * expmx + t rk = rk - one xk = xk - x end do debye2 = debinf / ( x * x ) - two * sum1 * expmx else if ( x < xlim1 ) then expmx = exp ( -x ) sum1 = ( ( x + two ) * x + two ) / ( x * x ) debye2 = debinf / ( x * x ) - two * sum1 * expmx else if ( x <= xlim2 ) then debye2 = debinf / ( x * x ) else debye2 = zero end if return end subroutine debye2_values ( n_data, x, fx ) !*****************************************************************************80 ! !! DEBYE2_VALUES returns some values of Debye's function of order 2. ! ! Discussion: ! ! The function is defined by: ! ! DEBYE2(x) = 2 / x^2 * Integral ( 0 <= t <= x ) t^2 / ( exp ( t ) - 1 ) dt ! ! The data was reported by McLeod. ! ! Modified: ! ! 27 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &