program main !******************************************************************************* ! !! TOMS757_PRB tests the routines in TOMS757. ! ! Discussion: ! ! This program tests the 37 functions in the MISCFUN package. ! It is a fairly simple code with each function being tested ! at 20 different arguments. The code compares the value ! from the function with a pre-computed value, and produces ! the absolute and relative errors. ! ! Author: ! ! Allan McLeod, ! Department of Mathematics and Statistics, ! Paisley University, High Street, Paisley, Scotland, PA12BE ! macl_ms0@paisley.ac.uk ! implicit none call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOMS757_PRB:' write ( *, '(a)' ) ' Test the uncommon special function routines in TOMS757.' call test01 call test02 call test03 call test04 call test05 call test06 call test07 call test08 call test09 call test10 call test11 call test12 call test13 call test14 call test15 call test16 call test17 call test18 call test19 call test20 call test21 call test22 call test23 call test24 call test25 call test26 call test27 call test28 call test29 call test30 call test31 call test32 call test33 call test34 call test35 call test36 call test37 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOMS757_PRB:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) !stop end subroutine test01 !******************************************************************************* ! !! TEST01 tests ABRAM0. ! implicit none real ( kind = 8 ) abram0 real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST01' write ( *, '(a)' ) ' Testing function ABRAM0' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call abram0_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = abram0 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test02 !******************************************************************************* ! !! TEST02 tests ABRAM1. ! implicit none ! real ( kind = 8 ) abram1 real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST02' write ( *, '(a)' ) ' Testing function ABRAM1' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call abram1_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = abram1 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test03 !******************************************************************************* ! !! TEST03 tests ABRAM2. ! implicit none real ( kind = 8 ) abram2 real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03' write ( *, '(a)' ) ' Testing function ABRAM2' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call abram2_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = abram2 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test04 !******************************************************************************* ! !! TEST04 tests AIRY_AI_INT. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) airy_ai_int real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST04' write ( *, '(a)' ) ' Testing function AIRY_AI_INT' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call airy_ai_int_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = airy_ai_int ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test05 !******************************************************************************* ! !! TEST05 tests AIRY_BI_INT. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) airy_bi_int real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05' write ( *, '(a)' ) ' Testing function AIRY_BI_INT' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call airy_bi_int_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = airy_bi_int ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test06 !******************************************************************************* ! !! TEST06 tests AIRY_GI. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) airy_gi real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST06' write ( *, '(a)' ) ' Testing function AIRY_GI' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call airy_gi_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = airy_gi ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test07 !******************************************************************************* ! !! TEST07 tests AIRY_HI. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) airy_hi real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST07' write ( *, '(a)' ) ' Testing function AIRY_HI' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call airy_hi_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = airy_hi ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test08 !******************************************************************************* ! !! TEST08 tests ARCTAN_INT. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) arctan_int real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08' write ( *, '(a)' ) ' Testing function ARCTAN_INT' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call arctan_int_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = arctan_int ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test09 !******************************************************************************* ! !! TEST09 tests BESSEL_I0_INT. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) bessel_i0_int real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST09' write ( *, '(a)' ) ' Testing function BESSEL_I0_INT' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call bessel_i0_int_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = bessel_i0_int ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test10 !******************************************************************************* ! !! TEST10 tests BESSEL_J0_INT. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) bessel_j0_int real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST10' write ( *, '(a)' ) ' Testing function BESSEL_J0_INT' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call bessel_j0_int_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = bessel_j0_int ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test11 !******************************************************************************* ! !! TEST11 tests BESSEL_K0_INT. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) bessel_k0_int real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11' write ( *, '(a)' ) ' Testing function BESSEL_K0_INT' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call bessel_k0_int_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = bessel_k0_int ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test12 !******************************************************************************* ! !! TEST12 tests BESSEL_Y0_INT. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) bessel_y0_int real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST12' write ( *, '(a)' ) ' Testing function BESSEL_Y0_INT' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call bessel_y0_int_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = bessel_y0_int ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test13 !******************************************************************************* ! !! TEST13 tests CLAUSEN. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) clausen real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST13' write ( *, '(a)' ) ' Testing function CLAUSEN' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call clausen_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = clausen ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test14 !******************************************************************************* ! !! TEST14 tests DEBYE1. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) debye1 real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST14' write ( *, '(a)' ) ' Testing function DEBYE1' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call debye1_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = debye1 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test15 !******************************************************************************* ! !! TEST15 tests DEBYE2. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) debye2 real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST15' write ( *, '(a)' ) ' Testing function DEBYE2' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call debye2_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = debye2 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test16 !******************************************************************************* ! !! TEST16 tests DEBYE3. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) debye3 real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST16' write ( *, '(a)' ) ' Testing function DEBYE3' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call debye3_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = debye3 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test17 !******************************************************************************* ! !! TEST17 tests DEBYE4. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) debye4 real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17' write ( *, '(a)' ) ' Testing function DEBYE4' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call debye4_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = debye4 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test18 !******************************************************************************* ! !! TEST18 tests EXP3_INT. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) exp3_int real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST18' write ( *, '(a)' ) ' Testing function EXP3_INT' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call exp3_int_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = exp3_int ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test19 !******************************************************************************* ! !! TEST19 tests GOODWIN. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx real ( kind = 8 ) goodwin integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST19' write ( *, '(a)' ) ' Testing function GOODWIN' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call goodwin_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = goodwin ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test20 !******************************************************************************* ! !! TEST20 tests I0ML0. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx real ( kind = 8 ) i0ml0 integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST20' write ( *, '(a)' ) ' Testing function I0ML0' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call i0ml0_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = i0ml0 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test21 !******************************************************************************* ! !! TEST21 tests I1ML1. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx real ( kind = 8 ) i1ml1 integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST21' write ( *, '(a)' ) ' Testing function I1ML1' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call i1ml1_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = i1ml1 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test22 !******************************************************************************* ! !! TEST22 tests LOBACHEVSKY. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx real ( kind = 8 ) lobachevsky integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST22' write ( *, '(a)' ) ' Testing function LOBACHEVSKY' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call lobachevsky_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = lobachevsky ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test23 !******************************************************************************* ! !! TEST23 tests STROMGEN. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) stromgen real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST23' write ( *, '(a)' ) ' Testing function STROMGEN' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call stromgen_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = stromgen ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test24 !******************************************************************************* ! !! TEST24 tests STRUVE_H0. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) struve_h0 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST24' write ( *, '(a)' ) ' Testing function STRUVE_H0' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call struve_h0_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = struve_h0 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test25 !******************************************************************************* ! !! TEST25 tests STRUVE_H1. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) struve_h1 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST25' write ( *, '(a)' ) ' Testing function STRUVE_H1' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call struve_h1_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = struve_h1 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test26 !******************************************************************************* ! !! TEST26 tests STRUVE_L0. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) struve_l0 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST26' write ( *, '(a)' ) ' Testing function STRUVE_L0' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call struve_l0_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = struve_l0 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test27 !******************************************************************************* ! !! TEST27 tests STRUVE_L1. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) struve_l1 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST27' write ( *, '(a)' ) ' Testing function STRUVE_L1' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call struve_l1_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = struve_l1 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test28 !******************************************************************************* ! !! TEST28 tests SYNCH1. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) synch1 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST28' write ( *, '(a)' ) ' Testing function SYNCH1' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call synch1_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = synch1 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test29 !******************************************************************************* ! !! TEST29 tests SYNCH2. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) synch2 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST29' write ( *, '(a)' ) ' Testing function SYNCH2' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call synch2_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = synch2 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test30 !******************************************************************************* ! !! TEST30 tests TRAN02. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) tran02 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST30' write ( *, '(a)' ) ' Testing function TRAN02' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call tran02_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = tran02 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test31 !******************************************************************************* ! !! TEST31 tests TRAN03. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) tran03 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST31' write ( *, '(a)' ) ' Testing function TRAN03' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call tran03_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = tran03 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test32 !******************************************************************************* ! !! TEST32 tests TRAN04. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) tran04 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST32' write ( *, '(a)' ) ' Testing function TRAN04' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call tran04_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = tran04 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test33 !******************************************************************************* ! !! TEST33 tests TRAN05. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) tran05 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST33' write ( *, '(a)' ) ' Testing function TRAN05' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call tran05_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = tran05 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test34 !******************************************************************************* ! !! TEST34 tests TRAN06. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) tran06 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST34' write ( *, '(a)' ) ' Testing function TRAN06' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call tran06_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = tran06 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test35 !******************************************************************************* ! !! TEST35 tests TRAN07. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) tran07 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST35' write ( *, '(a)' ) ' Testing function TRAN07' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call tran07_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = tran07 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test36 !******************************************************************************* ! !! TEST36 tests TRAN08. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) tran08 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST36' write ( *, '(a)' ) ' Testing function TRAN08' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call tran08_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = tran08 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end subroutine test37 !******************************************************************************* ! !! TEST37 tests TRAN09. ! implicit none real ( kind = 8 ) abserr real ( kind = 8 ) comp real ( kind = 8 ) fx integer n_data real ( kind = 8 ) relerr real ( kind = 8 ) tran09 real ( kind = 8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST37' write ( *, '(a)' ) ' Testing function TRAN09' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Argument Abs. error Rel. error' write ( *, '(a)' ) ' ' n_data = 0 do call tran09_values ( n_data, x, fx ) if ( n_data <= 0 ) then exit end if comp = tran09 ( x ) abserr = abs ( fx - comp ) relerr = abserr / abs ( fx ) write ( *, '(2x,f15.10,2x,d15.5,8x,d15.5)' ) x, abserr, relerr end do return end function abram0 ( xvalue ) !******************************************************************************* ! !! 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 = (/ & real ( -0.68121927093549469816d0, kind = 8 ), & real ( -0.78867919816149252495d0, kind = 8 ), & real ( 0.5121581776818819543d-1, kind = 8 ), & real ( -0.71092352894541296d-3, kind = 8 ), & real ( 0.368681808504287d-5, kind = 8 ), & real ( -0.917832337237d-8, kind = 8 ), & real ( 0.1270202563d-10, kind = 8 ), & real ( -0.1076888d-13, kind = 8 ), & real ( 0.599d-17, kind = 8 ) /) real ( kind = 8 ), dimension(0:8) :: ab0g = (/ & real ( -0.60506039430868273190d0, kind = 8 ), & real ( -0.41950398163201779803d0, kind = 8 ), & real ( 0.1703265125190370333d-1, kind = 8 ), & real ( -0.16938917842491397d-3, kind = 8 ), & real ( 0.67638089519710d-6, kind = 8 ), & real ( -0.135723636255d-8, kind = 8 ), & real ( 0.156297065d-11, kind = 8 ), & real ( -0.112887d-14, kind = 8 ), & real ( 0.55d-18, kind = 8 ) /) real ( kind = 8 ), dimension(0:8) :: ab0h = (/ & real ( 1.38202655230574989705d0, kind = 8 ), & real ( -0.30097929073974904355d0, kind = 8 ), & real ( 0.794288809364887241d-2, kind = 8 ), & real ( -0.6431910276847563d-4, kind = 8 ), & real ( 0.22549830684374d-6, kind = 8 ), & real ( -0.41220966195d-9, kind = 8 ), & real ( 0.44185282d-12, kind = 8 ), & real ( -0.30123d-15, kind = 8 ), & real ( 0.14d-18, kind = 8 ) /) real ( kind = 8 ), dimension(0:27) :: ab0as = (/ & real ( 1.97755499723693067407d+0, kind = 8 ), & real ( -0.1046024792004819485d-1, kind = 8 ), & real ( 0.69680790253625366d-3, kind = 8 ), & real ( -0.5898298299996599d-4, kind = 8 ), & real ( 0.577164455305320d-5, kind = 8 ), & real ( -0.61523013365756d-6, kind = 8 ), & real ( 0.6785396884767d-7, kind = 8 ), & real ( -0.723062537907d-8, kind = 8 ), & real ( 0.63306627365d-9, kind = 8 ), & real ( -0.989453793d-11, kind = 8 ), & real ( -0.1681980530d-10, kind = 8 ), & real ( 0.673799551d-11, kind = 8 ), & real ( -0.200997939d-11, kind = 8 ), & real ( 0.54055903d-12, kind = 8 ), & real ( -0.13816679d-12, kind = 8 ), & real ( 0.3422205d-13, kind = 8 ), & real ( -0.826686d-14, kind = 8 ), & real ( 0.194566d-14, kind = 8 ), & real ( -0.44268d-15, kind = 8 ), & real ( 0.9562d-16, kind = 8 ), & real ( -0.1883d-16, kind = 8 ), & real ( 0.301d-17, kind = 8 ), & real ( -0.19d-18, kind = 8 ), & real ( -0.14d-18, kind = 8 ), & real ( 0.11d-18, kind = 8 ), & real ( -0.4d-19, kind = 8 ), & real ( 0.2d-19, kind = 8 ), & real ( -0.1d-19, kind = 8 ) /) 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 ) !******************************************************************************* ! !! ABRAM0_VALUES returns some values of the Abramowitz0 function. ! ! Discussion: ! ! The function is defined by: ! ! ABRAM0(x) = Integral ( 0 <= t < infinity ) exp ( -t^2 - x / t ) dt ! ! The data was reported by McLeod. ! ! Modified: ! ! 21 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz and Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & real ( 0.87377726306985360531D+00, kind = 8 ), & real ( 0.84721859650456925922D+00, kind = 8 ), & real ( 0.77288934483988301615D+00, kind = 8 ), & real ( 0.59684345853450151603D+00, kind = 8 ), & real ( 0.29871735283675888392D+00, kind = 8 ), & real ( 0.15004596450516388138D+00, kind = 8 ), & real ( 0.11114662419157955096D+00, kind = 8 ), & real ( 0.83909567153151897766D-01, kind = 8 ), & real ( 0.56552321717943417515D-01, kind = 8 ), & real ( 0.49876496603033790206D-01, kind = 8 ), & real ( 0.44100889219762791328D-01, kind = 8 ), & real ( 0.19738535180254062496D-01, kind = 8 ), & real ( 0.86193088287161479900D-02, kind = 8 ), & real ( 0.40224788162540127227D-02, kind = 8 ), & real ( 0.19718658458164884826D-02, kind = 8 ), & real ( 0.10045868340133538505D-02, kind = 8 ), & real ( 0.15726917263304498649D-03, kind = 8 ), & real ( 0.10352666912350263437D-04, kind = 8 ), & real ( 0.91229759190956745069E-06, kind = 8 ), & real ( 0.25628287737952698742E-09, kind = 8 ) /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & real ( 0.0019531250D+00, kind = 8 ), & real ( 0.0078125000D+00, kind = 8 ), & real ( 0.0312500000D+00, kind = 8 ), & real ( 0.1250000000D+00, kind = 8 ), & real ( 0.5000000000D+00, kind = 8 ), & real ( 1.0000000000D+00, kind = 8 ), & real ( 1.2500000000D+00, kind = 8 ), & real ( 1.5000000000D+00, kind = 8 ), & real ( 1.8750000000D+00, kind = 8 ), & real ( 2.0000000000D+00, kind = 8 ), & real ( 2.1250000000D+00, kind = 8 ), & real ( 3.0000000000D+00, kind = 8 ), & real ( 4.0000000000D+00, kind = 8 ), & real ( 5.0000000000D+00, kind = 8 ), & real ( 6.0000000000D+00, kind = 8 ), & real ( 7.0000000000D+00, kind = 8 ), & real ( 10.0000000000D+00, kind = 8 ), & real ( 15.0000000000D+00, kind = 8 ), & real ( 20.0000000000D+00, kind = 8 ), & real ( 40.0000000000D+00, kind = 8 ) /) 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 ) !******************************************************************************* ! !! 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 ) !******************************************************************************* ! !! ABRAM1_VALUES returns some values of the Abramowitz1 function. ! ! Discussion: ! ! The function is defined by: ! ! ABRAM1(x) = Integral ( 0 <= t < infinity ) t * exp ( -t^2 - x / t ) dt ! ! The data was reported by McLeod. ! ! Modified: ! ! 21 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz and Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & real ( 0.49828219848799921792D+00, kind = 8 ), & real ( 0.49324391773047288556D+00, kind = 8 ), & real ( 0.47431612784691234649D+00, kind = 8 ), & real ( 0.41095983258760410149D+00, kind = 8 ), & real ( 0.25317617388227035867D+00, kind = 8 ), & real ( 0.14656338138597777543D+00, kind = 8 ), & real ( 0.11421547056018366587D+00, kind = 8 ), & real ( 0.90026307383483764795D-01, kind = 8 ), & real ( 0.64088214170742303375D-01, kind = 8 ), & real ( 0.57446614314166191085D-01, kind = 8 ), & real ( 0.51581624564800730959D-01, kind = 8 ), & real ( 0.25263719555776416016D-01, kind = 8 ), & real ( 0.11930803330196594536D-01, kind = 8 ), & real ( 0.59270542280915272465D-02, kind = 8 ), & real ( 0.30609215358017829567D-02, kind = 8 ), & real ( 0.16307382136979552833D-02, kind = 8 ), & real ( 0.28371851916959455295D-03, kind = 8 ), & real ( 0.21122150121323238154D-04, kind = 8 ), & real ( 0.20344578892601627337D-05, kind = 8 ), & real ( 0.71116517236209642290E-09, kind = 8 ) /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & real ( 0.0019531250D+00, kind = 8 ), & real ( 0.0078125000D+00, kind = 8 ), & real ( 0.0312500000D+00, kind = 8 ), & real ( 0.1250000000D+00, kind = 8 ), & real ( 0.5000000000D+00, kind = 8 ), & real ( 1.0000000000D+00, kind = 8 ), & real ( 1.2500000000D+00, kind = 8 ), & real ( 1.5000000000D+00, kind = 8 ), & real ( 1.8750000000D+00, kind = 8 ), & real ( 2.0000000000D+00, kind = 8 ), & real ( 2.1250000000D+00, kind = 8 ), & real ( 3.0000000000D+00, kind = 8 ), & real ( 4.0000000000D+00, kind = 8 ), & real ( 5.0000000000D+00, kind = 8 ), & real ( 6.0000000000D+00, kind = 8 ), & real ( 7.0000000000D+00, kind = 8 ), & real ( 10.0000000000D+00, kind = 8 ), & real ( 15.0000000000D+00, kind = 8 ), & real ( 20.0000000000D+00, kind = 8 ), & real ( 40.0000000000D+00, kind = 8 ) /) 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 ) !******************************************************************************* ! !! 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 = real ( 6.0D+00, kind = 8 ) 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 ) !******************************************************************************* ! !! ABRAM2_VALUES returns some values of the Abramowitz2 function. ! ! Discussion: ! ! The function is defined by: ! ! ABRAM2(x) = Integral ( 0 <= t < infinity ) t^2 * exp ( -t^2 - x / t ) dt ! ! The data was reported by McLeod. ! ! Modified: ! ! 22 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz and Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & real ( 0.44213858162107913430D+00, kind = 8 ), & real ( 0.43923379545684026308D+00, kind = 8 ), & real ( 0.42789857297092602234D+00, kind = 8 ), & real ( 0.38652825661854504406D+00, kind = 8 ), & real ( 0.26538204413231368110D+00, kind = 8 ), & real ( 0.16848734838334595000D+00, kind = 8 ), & real ( 0.13609200032513227112D+00, kind = 8 ), & real ( 0.11070330027727917352D+00, kind = 8 ), & real ( 0.82126019995530382267D-01, kind = 8 ), & real ( 0.74538781999594581763D-01, kind = 8 ), & real ( 0.67732034377612811390D-01, kind = 8 ), & real ( 0.35641808698811851022D-01, kind = 8 ), & real ( 0.17956589956618269083D-01, kind = 8 ), & real ( 0.94058737143575370625D-02, kind = 8 ), & real ( 0.50809356204299213556D-02, kind = 8 ), & real ( 0.28149565414209719359D-02, kind = 8 ), & real ( 0.53808696422559303431D-03, kind = 8 ), & real ( 0.44821756380146327259D-04, kind = 8 ), & real ( 0.46890678427324100410D-05, kind = 8 ), & real ( 0.20161544850996420504D-08, kind = 8 ) /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & real ( 0.0019531250D+00, kind = 8 ), & real ( 0.0078125000D+00, kind = 8 ), & real ( 0.0312500000D+00, kind = 8 ), & real ( 0.1250000000D+00, kind = 8 ), & real ( 0.5000000000D+00, kind = 8 ), & real ( 1.0000000000D+00, kind = 8 ), & real ( 1.2500000000D+00, kind = 8 ), & real ( 1.5000000000D+00, kind = 8 ), & real ( 1.8750000000D+00, kind = 8 ), & real ( 2.0000000000D+00, kind = 8 ), & real ( 2.1250000000D+00, kind = 8 ), & real ( 3.0000000000D+00, kind = 8 ), & real ( 4.0000000000D+00, kind = 8 ), & real ( 5.0000000000D+00, kind = 8 ), & real ( 6.0000000000D+00, kind = 8 ), & real ( 7.0000000000D+00, kind = 8 ), & real ( 10.0000000000D+00, kind = 8 ), & real ( 15.0000000000D+00, kind = 8 ), & real ( 20.0000000000D+00, kind = 8 ), & real ( 40.0000000000D+00, kind = 8 ) /) 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 ) !******************************************************************************* ! !! 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 ) !******************************************************************************* ! !! AIRY_AI_INT_VALUES returns some values of the integral of the Airy function. ! ! Discussion: ! ! The function is defined by: ! ! AIRY_AI_INT(x) = Integral ( 0 <= t <= x ) Ai(t) dt ! ! The data was reported by McLeod. ! ! Modified: ! ! 22 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz and Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & real ( -0.75228838916610124300D+00, kind = 8 ), & real ( -0.57348350185854889466D+00, kind = 8 ), & real ( -0.76569840313421291743D+00, kind = 8 ), & real ( -0.65181015505382467421D+00, kind = 8 ), & real ( -0.55881974894471876922D+00, kind = 8 ), & real ( -0.56902352870716815309D+00, kind = 8 ), & real ( -0.47800749642926168100D+00, kind = 8 ), & real ( -0.46567398346706861416D+00, kind = 8 ), & real ( -0.96783140945618013679D-01, kind = 8 ), & real ( -0.34683049857035607494D-03, kind = 8 ), & real ( 0.34658366917927930790D-03, kind = 8 ), & real ( 0.27657581846051227124D-02, kind = 8 ), & real ( 0.14595330491185717833D+00, kind = 8 ), & real ( 0.23631734191710977960D+00, kind = 8 ), & real ( 0.33289264538612212697D+00, kind = 8 ), & real ( 0.33318759129779422976D+00, kind = 8 ), & real ( 0.33332945170523851439D+00, kind = 8 ), & real ( 0.33333331724248357420D+00, kind = 8 ), & real ( 0.33333333329916901594D+00, kind = 8 ), & real ( 0.33333333333329380187D+00, kind = 8 ) /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & real ( -12.0000000000D+00, kind = 8 ), & real ( -11.0000000000D+00, kind = 8 ), & real ( -10.0000000000D+00, kind = 8 ), & real ( -9.5000000000D+00, kind = 8 ), & real ( -9.0000000000D+00, kind = 8 ), & real ( -6.5000000000D+00, kind = 8 ), & real ( -4.0000000000D+00, kind = 8 ), & real ( -1.0000000000D+00, kind = 8 ), & real ( -0.2500000000D+00, kind = 8 ), & real ( -0.0009765625D+00, kind = 8 ), & real ( 0.0009765625D+00, kind = 8 ), & real ( 0.0078125000D+00, kind = 8 ), & real ( 0.5000000000D+00, kind = 8 ), & real ( 1.0000000000D+00, kind = 8 ), & real ( 4.0000000000D+00, kind = 8 ), & real ( 4.5000000000D+00, kind = 8 ), & real ( 6.0000000000D+00, kind = 8 ), & real ( 8.0000000000D+00, kind = 8 ), & real ( 10.0000000000D+00, kind = 8 ), & real ( 12.0000000000D+00, kind = 8 ) /) 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 ) !******************************************************************************* ! !! 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 = real ( 8.0, kind = 8 ) real ( kind = 8 ), parameter :: four = real ( 4.0D+00, kind = 8 ) 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 ) !******************************************************************************* ! !! AIRY_BI_INT_VALUES returns some values of the integral of the Airy function. ! ! Discussion: ! ! The function is defined by: ! ! AIRY_BI_INT(x) = Integral ( 0 <= t <= x ) Bi(t) dt ! ! The data was reported by McLeod. ! ! Modified: ! ! 23 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz and Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & real ( 0.17660819031554631869D-01, kind = 8 ), & real ( -0.15040424806140020451D-01, kind = 8 ), & real ( 0.14756446293227661920D-01, kind = 8 ), & real ( -0.11847304264848446271D+00, kind = 8 ), & real ( -0.64916741266165856037D-01, kind = 8 ), & real ( 0.97260832464381044540D-01, kind = 8 ), & real ( 0.50760058495287539119D-01, kind = 8 ), & real ( -0.37300500963429492179D+00, kind = 8 ), & real ( -0.13962988442666578531D+00, kind = 8 ), & real ( -0.12001735266723296160D-02, kind = 8 ), & real ( 0.12018836117890354598D-02, kind = 8 ), & real ( 0.36533846550952011043D+00, kind = 8 ), & real ( 0.87276911673800812196D+00, kind = 8 ), & real ( 0.48219475263803429675D+02, kind = 8 ), & real ( 0.44006525804904178439D+06, kind = 8 ), & real ( 0.17608153976228301458D+07, kind = 8 ), & real ( 0.73779211705220007228D+07, kind = 8 ), & real ( 0.14780980310740671617D+09, kind = 8 ), & real ( 0.97037614223613433849D+11, kind = 8 ), & real ( 0.11632737638809878460D+15, kind = 8 ) /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & real ( -12.0000000000D+00, kind = 8 ), & real ( -10.0000000000D+00, kind = 8 ), & real ( -8.0000000000D+00, kind = 8 ), & real ( -7.5000000000D+00, kind = 8 ), & real ( -7.0000000000D+00, kind = 8 ), & real ( -6.5000000000D+00, kind = 8 ), & real ( -4.0000000000D+00, kind = 8 ), & real ( -1.0000000000D+00, kind = 8 ), & real ( -0.2500000000D+00, kind = 8 ), & real ( -0.0019531250D+00, kind = 8 ), & real ( 0.0019531250D+00, kind = 8 ), & real ( 0.5000000000D+00, kind = 8 ), & real ( 1.0000000000D+00, kind = 8 ), & real ( 4.0000000000D+00, kind = 8 ), & real ( 8.0000000000D+00, kind = 8 ), & real ( 8.5000000000D+00, kind = 8 ), & real ( 9.0000000000D+00, kind = 8 ), & real ( 10.0000000000D+00, kind = 8 ), & real ( 12.0000000000D+00, kind = 8 ), & real ( 14.0000000000D+00, kind = 8 ) /) 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 ) !******************************************************************************* ! !! 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 = real ( 4.0D+00, kind = 8 ) integer, parameter :: nterm1 = 28 integer, parameter :: nterm2 = 23 integer, parameter :: nterm3 = 39 integer nterm4 integer nterm5 integer nterm6 real ( kind = 8 ), parameter :: one = real ( 1.0, kind = 8 ) 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 ) !******************************************************************************* ! !! AIRY_GI_VALUES returns some values of the Airy Gi function. ! ! Discussion: ! ! The function is defined by: ! ! AIRY_GI(x) = Integral ( 0 <= t < infinity ) sin ( x*t+t^3/3) dt / pi ! ! The data was reported by McLeod. ! ! Modified: ! ! 24 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz and Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & real ( 0.20468308070040542435D+00, kind = 8 ), & real ( 0.18374662832557904078D+00, kind = 8 ), & real ( -0.11667221729601528265D+00, kind = 8 ), & real ( 0.31466934902729557596D+00, kind = 8 ), & real ( -0.37089040722426257729D+00, kind = 8 ), & real ( -0.25293059772424019694D+00, kind = 8 ), & real ( 0.28967410658692701936D+00, kind = 8 ), & real ( -0.34644836492634090590D+00, kind = 8 ), & real ( 0.28076035913873049496D+00, kind = 8 ), & real ( 0.21814994508094865815D+00, kind = 8 ), & real ( 0.20526679000810503329D+00, kind = 8 ), & real ( 0.22123695363784773258D+00, kind = 8 ), & real ( 0.23521843981043793760D+00, kind = 8 ), & real ( 0.82834303363768729338D-01, kind = 8 ), & real ( 0.45757385490989281893D-01, kind = 8 ), & real ( 0.44150012014605159922D-01, kind = 8 ), & real ( 0.39951133719508907541D-01, kind = 8 ), & real ( 0.35467706833949671483D-01, kind = 8 ), & real ( 0.31896005100679587981D-01, kind = 8 ), & real ( 0.26556892713512410405D-01, kind = 8 ) /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & real ( -0.0019531250D+00, kind = 8 ), & real ( -0.1250000000D+00, kind = 8 ), & real ( -1.0000000000D+00, kind = 8 ), & real ( -4.0000000000D+00, kind = 8 ), & real ( -8.0000000000D+00, kind = 8 ), & real ( -8.2500000000D+00, kind = 8 ), & real ( -9.0000000000D+00, kind = 8 ), & real ( -10.0000000000D+00, kind = 8 ), & real ( -11.0000000000D+00, kind = 8 ), & real ( -13.0000000000D+00, kind = 8 ), & real ( 0.0019531250D+00, kind = 8 ), & real ( 0.1250000000D+00, kind = 8 ), & real ( 1.0000000000D+00, kind = 8 ), & real ( 4.0000000000D+00, kind = 8 ), & real ( 7.0000000000D+00, kind = 8 ), & real ( 7.2500000000D+00, kind = 8 ), & real ( 8.0000000000D+00, kind = 8 ), & real ( 9.0000000000D+00, kind = 8 ), & real ( 10.0000000000D+00, kind = 8 ), & real ( 12.0000000000D+00, kind = 8 ) /) 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 ) !******************************************************************************* ! !! 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 = real ( 4.0D+00, kind = 8 ) integer, parameter :: nterm1 = 29 integer, parameter :: nterm2 = 17 integer, parameter :: nterm3 = 22 integer nterm4 integer nterm5 real ( kind = 8 ), parameter :: one = real ( 1.0, kind = 8 ) 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 ) !******************************************************************************* ! !! AIRY_HI_VALUES returns some values of the Airy Hi function. ! ! Discussion: ! ! The function is defined by: ! ! AIRY_HI(x) = Integral ( 0 <= t < infinity ) exp ( x * t - t^3 / 3 ) dt / pi ! ! The data was reported by McLeod. ! ! Modified: ! ! 24 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz and Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Allan McLeod, ! Algorithm 757, MISCFUN: A software package to compute uncommon ! special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 before the ! first call. On each call, the routine increments N_DATA by 1, and ! returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = 8 ) X, the argument of the function. ! ! Output, real ( kind = 8 ) FX, the value of the function. ! implicit none integer, parameter :: n_max = 20 real ( kind = 8 ) fx real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ & real ( 0.40936798278458884024D+00, kind = 8 ), & real ( 0.37495291608048868619D+00, kind = 8 ), & real ( 0.22066960679295989454D+00, kind = 8 ), & real ( 0.77565356679703713590D-01, kind = 8 ), & real ( 0.39638826473124717315D-01, kind = 8 ), & real ( 0.38450072575004151871D-01, kind = 8 ), & real ( 0.35273216868317898556D-01, kind = 8 ), & real ( 0.31768535282502272742D-01, kind = 8 ), & real ( 0.28894408288051391369D-01, kind = 8 ), & real ( 0.24463284011678541180D-01, kind = 8 ), & real ( 0.41053540139998941517D+00, kind = 8 ), & real ( 0.44993502381204990817D+00, kind = 8 ), & real ( 0.97220515514243332184D+00, kind = 8 ), & real ( 0.83764237105104371193D+02, kind = 8 ), & real ( 0.80327744952044756016D+05, kind = 8 ), & real ( 0.15514138847749108298D+06, kind = 8 ), & real ( 0.11995859641733262114D+07, kind = 8 ), & real ( 0.21472868855967642259D+08, kind = 8 ), & real ( 0.45564115351632913590D+09, kind = 8 ), & real ( 0.32980722582904761929D+12, kind = 8 ) /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & real ( -0.0019531250D+00, kind = 8 ), & real ( -0.1250000000D+00, kind = 8 ), & real ( -1.0000000000D+00, kind = 8 ), & real ( -4.0000000000D+00, kind = 8 ), & real ( -8.0000000000D+00, kind = 8 ), & real ( -8.2500000000D+00, kind = 8 ), & real ( -9.0000000000D+00, kind = 8 ), & real ( -10.0000000000D+00, kind = 8 ), & real ( -11.0000000000D+00, kind = 8 ), & real ( -13.0000000000D+00, kind = 8 ), & real ( 0.0019531250D+00, kind = 8 ), & real ( 0.1250000000D+00, kind = 8 ), & real ( 1.0000000000D+00, kind = 8 ), & real ( 4.0000000000D+00, kind = 8 ), & real ( 7.0000000000D+00, kind = 8 ), & real ( 7.2500000000D+00, kind = 8 ), & real ( 8.0000000000D+00, kind = 8 ), & real ( 9.0000000000D+00, kind = 8 ), & real ( 10.0000000000D+00, kind = 8 ), & real ( 12.0000000000D+00, kind = 8 ) /) 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 ) !******************************************************************************* ! !! 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 = real ( 1.0, kind = 8 ) real ( kind = 8 ) t real ( kind = 8 ), parameter :: twobpi = real ( 0.63661977236758134308d0, & kind = 8 ) real ( kind = 8 ) x real ( kind = 8 ), parameter :: xlow = real ( 7.4505806d-9, kind = 8 ) real ( kind = 8 ), parameter :: xupper = real ( 4.5036d15, kind = 8 ) 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 ) !******************************************************************************* ! !! 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 = (/ & real ( 0.19531241721588483191D-02, kind = 8 ), & real ( -0.39062433772980711281D-02, kind = 8 ), & real ( 0.78124470192576499535D-02, kind = 8 ), & real ( 0.15624576181996527280D-01, kind = 8 ), & real ( -0.31246610349485401551D-01, kind = 8 ), & real ( 0.62472911335014397321D-01, kind = 8 ), & real ( 0.12478419717389654039D+00, kind = 8 ), & real ( -0.24830175098230686908D+00, kind = 8 ), & real ( 0.48722235829452235711D+00, kind = 8 ), & real ( 0.91596559417721901505D+00, kind = 8 ), & real ( 0.12749694484943800618D+01, kind = 8 ), & real ( -0.15760154034463234224D+01, kind = 8 ), & real ( 0.24258878412859089996D+01, kind = 8 ), & real ( 0.33911633326292997361D+01, kind = 8 ), & real ( 0.44176450919422186583D+01, kind = 8 ), & real ( -0.47556713749547247774D+01, kind = 8 ), & real ( 0.50961912150934111303D+01, kind = 8 ), & real ( 0.53759175735714876256D+01, kind = 8 ), & real ( -0.61649904785027487422D+01, kind = 8 ), & real ( 0.72437843013083534973D+01, kind = 8 ) /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & real ( 0.0019531250D+00, kind = 8 ), & real ( -0.0039062500D+00, kind = 8 ), & real ( 0.0078125000D+00, kind = 8 ), & real ( 0.0156250000D+00, kind = 8 ), & real ( -0.0312500000D+00, kind = 8 ), & real ( 0.0625000000D+00, kind = 8 ), & real ( 0.1250000000D+00, kind = 8 ), & real ( -0.2500000000D+00, kind = 8 ), & real ( 0.5000000000D+00, kind = 8 ), & real ( 1.0000000000D+00, kind = 8 ), & real ( 1.5000000000D+00, kind = 8 ), & real ( -2.0000000000D+00, kind = 8 ), & real ( 4.0000000000D+00, kind = 8 ), & real ( 8.0000000000D+00, kind = 8 ), & real ( 16.0000000000D+00, kind = 8 ), & real ( -20.0000000000D+00, kind = 8 ), & real ( 25.0000000000D+00, kind = 8 ), & real ( 30.0000000000D+00, kind = 8 ), & real ( -50.0000000000D+00, kind = 8 ), & real ( 100.0000000000D+00, kind = 8 ) /) 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 ) !******************************************************************************* ! !! 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 = real ( 18.0D+00, kind = 8 ) 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 = & real ( 0.91893853320467274178d0, kind = 8 ) integer, parameter :: nterm1 = 25 integer, parameter :: nterm2 = 27 real ( kind = 8 ) t real ( kind = 8 ) temp real ( kind = 8 ), parameter :: thirt6 = real ( 36.0D+00, kind = 8 ) real ( kind = 8 ), parameter :: three = 3.0D+00 real ( kind = 8 ) x real ( kind = 8 ), parameter :: xhigh = real ( 713.758339d0, kind = 8 ) real ( kind = 8 ), parameter :: xlow = real ( 0.5161914d-7, kind = 8 ) 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 ) !******************************************************************************* ! !! 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 = (/ & real ( 0.19531256208818052282D-02, kind = 8 ), & real ( -0.39062549670565734544D-02, kind = 8 ), & real ( 0.62520348032546565850D-01, kind = 8 ), & real ( 0.12516285581366971819D+00, kind = 8 ), & real ( -0.51051480879740303760D+00, kind = 8 ), & real ( 0.10865210970235898158D+01, kind = 8 ), & real ( 0.27750019054282535299D+01, kind = 8 ), & real ( -0.13775208868039716639D+02, kind = 8 ), & real ( 0.46424372058106108576D+03, kind = 8 ), & real ( 0.64111867658021584522D+07, kind = 8 ), & real ( -0.10414860803175857953D+08, kind = 8 ), & real ( 0.44758598913855743089D+08, kind = 8 ), & real ( -0.11852985311558287888D+09, kind = 8 ), & real ( 0.31430078220715992752D+09, kind = 8 ), & real ( -0.83440212900794309620D+09, kind = 8 ), & real ( 0.22175367579074298261D+10, kind = 8 ), & real ( 0.58991731842803636487D+10, kind = 8 ), & real ( -0.41857073244691522147D+11, kind = 8 ), & real ( 0.79553885818472357663D+12, kind = 8 ), & real ( 0.15089715082719201025D+17, kind = 8 ) /) integer n_data real ( kind = 8 ) x real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ & real ( 0.0019531250D+00, kind = 8 ), & real ( -0.0039062500D+00, kind = 8 ), & real ( 0.0625000000D+00, kind = 8 ), & real ( 0.1250000000D+00, kind = 8 ), & real ( -0.5000000000D+00, kind = 8 ), & real ( 1.0000000000D+00, kind = 8 ), & real ( 2.0000000000D+00, kind = 8 ), & real ( -4.0000000000D+00, kind = 8 ), & real ( 8.0000000000D+00, kind = 8 ), & real ( 18.0000000000D+00, kind = 8 ), & real ( -18.5000000000D+00, kind = 8 ), & real ( 20.0000000000D+00, kind = 8 ), & real ( -21.0000000000D+00, kind = 8 ), & real ( 22.0000000000D+00, kind = 8 ), & real ( -23.0000000000D+00, kind = 8 ), & real ( 24.0000000000D+00, kind = 8 ), & real ( 25.0000000000D+00, kind = 8 ), & real ( -27.0000000000D+00, kind = 8 ), & real ( 30.0000000000D+00, kind = 8 ), & real ( 40.0000000000D+00, kind = 8 ) /) 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 ) !******************************************************************************* ! !! 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 ) !******************************************************************************* ! !! 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 = (/ & real ( 0.97656242238978822427D-03, kind = 8 ), & real ( 0.39062450329491108875D-02, kind = 8 ), & real ( -0.62479657927917933620D-01, kind = 8 ), & real ( 0.12483733492120479139D+00, kind = 8 ), & real ( -0.48968050664604505505D+00, kind = 8 ), & real ( 0.91973041008976023931D+00, kind = 8 ), & real ( -0.14257702931970265690D+01, kind = 8 ), & real ( 0.10247341594606064818D+01, kind = 8 ), & real ( -0.12107468348304501655D+01, kind = 8 ), & real ( 0.11008652032736190799D+01, kind = 8 ), & real ( -0.10060334829904124192D+01, kind = 8 ), & real ( 0.81330572662485953519D+00, kind = 8 ), & real ( -0.10583788214211277585D+01, kind = 8 ), & real ( 0.87101492116545875169D+00, kind = 8 ), & real ( -0.88424908882