Bessel functions for Gaussian scatter
This commit is contained in:
Martin Diehl
parent
408524daac
commit
230ffe8f35
131
src/math.f90
131
src/math.f90
|
|
@ -2822,4 +2822,135 @@ real(pReal) pure function math_limit(a, left, right)
|
|||
|
||||
end function math_limit
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Modified Bessel I function of order 0
|
||||
!> @author John Burkardt
|
||||
!> @details original version available on https://people.sc.fsu.edu/~jburkardt/f_src/toms715/toms715.html
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
real(pReal) function bessel_i0 (x)
|
||||
use, intrinsic :: IEEE_ARITHMETIC
|
||||
|
||||
implicit none
|
||||
real(pReal), intent(in) :: x
|
||||
integer(pInt) :: i
|
||||
real(pReal) :: sump_p, sump_q, xAbs, xx
|
||||
real(pReal), parameter, dimension(15) :: p_small = real( &
|
||||
[-5.2487866627945699800e-18, -1.5982226675653184646e-14, -2.6843448573468483278e-11, &
|
||||
-3.0517226450451067446e-08, -2.5172644670688975051e-05, -1.5453977791786851041e-02, &
|
||||
-7.0935347449210549190e+00, -2.4125195876041896775e+03, -5.9545626019847898221e+05, &
|
||||
-1.0313066708737980747e+08, -1.1912746104985237192e+10, -8.4925101247114157499e+11, &
|
||||
-3.2940087627407749166e+13, -5.5050369673018427753e+14, -2.2335582639474375249e+15], pReal)
|
||||
real(pReal), parameter, dimension(5) :: q_small = real( &
|
||||
[-3.7277560179962773046e+03, 6.5158506418655165707e+06, -6.5626560740833869295e+09, &
|
||||
3.7604188704092954661e+12, -9.7087946179594019126e+14], pReal)
|
||||
real(pReal), parameter, dimension(8) :: p_large = real( &
|
||||
[-3.9843750000000000000e-01, 2.9205384596336793945e+00, -2.4708469169133954315e+00, &
|
||||
4.7914889422856814203e-01, -3.7384991926068969150e-03, -2.6801520353328635310e-03, &
|
||||
9.9168777670983678974e-05, -2.1877128189032726730e-06], pReal)
|
||||
real(pReal), parameter, dimension(7) :: q_large = real( &
|
||||
[-3.1446690275135491500e+01, 8.5539563258012929600e+01, -6.0228002066743340583e+01, &
|
||||
1.3982595353892851542e+01, -1.1151759188741312645e+00, 3.2547697594819615062e-02, &
|
||||
-5.5194330231005480228e-04], pReal)
|
||||
|
||||
|
||||
xAbs = abs(x)
|
||||
|
||||
argRange: if (xAbs < 5.55e-17_pReal) then
|
||||
bessel_i0 = 1.0_pReal
|
||||
else if (xAbs < 15.0_pReal) then argRange
|
||||
xx = xAbs**2.0_pReal
|
||||
sump_p = p_small(1)
|
||||
do i = 2, 15
|
||||
sump_p = sump_p * xx + p_small(i)
|
||||
end do
|
||||
xx = xx - 225.0_pReal
|
||||
sump_q = ((((xx+q_small(1))*xx+q_small(2))*xx+q_small(3))*xx+q_small(4))*xx+q_small(5)
|
||||
bessel_i0 = sump_p / sump_q
|
||||
else if (xAbs <= 713.986_pReal) then argRange
|
||||
xx = 1.0_pReal / xAbs - 2.0_pReal/30.0_pReal
|
||||
sump_p = ((((((p_large(1)*xx+p_large(2))*xx+p_large(3))*xx+p_large(4))*xx+ &
|
||||
p_large(5))*xx+p_large(6))*xx+p_large(7))*xx+p_large(8)
|
||||
sump_q = ((((((xx+q_large(1))*xx+q_large(2))*xx+q_large(3))*xx+ &
|
||||
q_large(4))*xx+q_large(5))*xx+q_large(6))*xx+q_large(7)
|
||||
bessel_i0 = sump_p / sump_q
|
||||
|
||||
avoidOverflow: if (xAbs > 698.986_pReal) then
|
||||
bessel_i0 = ((bessel_i0*exp(xAbs-40.0_pReal)-p_large(1)*exp(xAbs-40.0_pReal))/sqrt(xAbs))*exp(40.0)
|
||||
else avoidOverflow
|
||||
bessel_i0 = ((bessel_i0*exp(xAbs)-p_large(1)*exp(xAbs))/sqrt(xAbs))
|
||||
endif avoidOverflow
|
||||
|
||||
else argRange
|
||||
bessel_i0 = IEEE_value(bessel_i0,IEEE_positive_inf)
|
||||
end if argRange
|
||||
|
||||
end function bessel_i0
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Modified Bessel I function of order 1
|
||||
!> @author John Burkardt
|
||||
!> @details original version available on https://people.sc.fsu.edu/~jburkardt/f_src/toms715/toms715.html
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
real(pReal) function bessel_i1 (x)
|
||||
use, intrinsic :: IEEE_ARITHMETIC
|
||||
|
||||
implicit none
|
||||
real(pReal), intent(in) :: x
|
||||
integer(pInt) :: i
|
||||
real(pReal) :: sump_p, sump_q, xAbs, xx
|
||||
real(pReal), dimension(15), parameter :: p_small = real( &
|
||||
[-1.9705291802535139930e-19, -6.5245515583151902910e-16, -1.1928788903603238754e-12, &
|
||||
-1.4831904935994647675e-09, -1.3466829827635152875e-06, -9.1746443287817501309e-04, &
|
||||
-4.7207090827310162436e-01, -1.8225946631657315931e+02, -5.1894091982308017540e+04, &
|
||||
-1.0588550724769347106e+07, -1.4828267606612366099e+09, -1.3357437682275493024e+11, &
|
||||
-6.9876779648010090070e+12, -1.7732037840791591320e+14, -1.4577180278143463643e+15], pReal)
|
||||
real(pReal), dimension(5), parameter :: q_small = real( &
|
||||
[-4.0076864679904189921e+03, 7.4810580356655069138e+06, -8.0059518998619764991e+09, &
|
||||
4.8544714258273622913e+12, -1.3218168307321442305e+15], pReal)
|
||||
real(pReal), dimension(8), parameter :: p_large = real( &
|
||||
[-6.0437159056137600000e-02, 4.5748122901933459000e-01, -4.2843766903304806403e-01, &
|
||||
9.7356000150886612134e-02, -3.2457723974465568321e-03, -3.6395264712121795296e-04, &
|
||||
1.6258661867440836395e-05, -3.6347578404608223492e-07], pReal)
|
||||
real(pReal), dimension(6), parameter :: q_large = real( &
|
||||
[-3.8806586721556593450e+00, 3.2593714889036996297e+00, -8.5017476463217924408e-01, &
|
||||
7.4212010813186530069e-02, -2.2835624489492512649e-03, 3.7510433111922824643e-05], pReal)
|
||||
real(pReal), parameter :: pbar = 3.98437500e-01
|
||||
|
||||
|
||||
xAbs = abs(x)
|
||||
|
||||
argRange: if (xAbs < 5.55e-17_pReal) then
|
||||
bessel_i1 = 0.5_pReal * xAbs
|
||||
else if (xAbs < 15.0_pReal) then argRange
|
||||
xx = xAbs**2.0_pReal
|
||||
sump_p = p_small(1)
|
||||
do i = 2, 15
|
||||
sump_p = sump_p * xx + p_small(i)
|
||||
end do
|
||||
xx = xx - 225.0_pReal
|
||||
sump_q = ((((xx+q_small(1))*xx+q_small(2))*xx+q_small(3))*xx+q_small(4)) * xx + q_small(5)
|
||||
bessel_i1 = (sump_p / sump_q) * xAbs
|
||||
else if (xAbs <= 713.986_pReal) then argRange
|
||||
xx = 1.0_pReal / xAbs - 2.0_pReal/30.0_pReal
|
||||
sump_p = ((((((p_large(1)*xx+p_large(2))*xx+p_large(3))*xx+p_large(4))*xx+&
|
||||
p_large(5))*xx+p_large(6))*xx+p_large(7))*xx+p_large(8)
|
||||
sump_q = (((((xx+q_large(1))*xx+q_large(2))*xx+q_large(3))*xx+ q_large(4))*xx+q_large(5))*xx+q_large(6)
|
||||
bessel_i1 = sump_p / sump_q
|
||||
|
||||
avoidOverflow: if (xAbs > 698.986_pReal) then
|
||||
bessel_i1 = ((bessel_i1 * exp(xAbs-40.0_pReal) + pbar * exp(xAbs-40.0_pReal)) / sqrt(xAbs)) * exp(40.0_pReal)
|
||||
else avoidOverflow
|
||||
bessel_i1 = ((bessel_i1 * exp(xAbs) + pbar * exp(xAbs)) / sqrt(xAbs))
|
||||
endif avoidOverflow
|
||||
|
||||
else argRange
|
||||
bessel_i1 = IEEE_value(bessel_i1,IEEE_positive_inf)
|
||||
end if argRange
|
||||
|
||||
if (x < 0.0_pReal) bessel_i1 = -bessel_i1
|
||||
|
||||
end function bessel_i1
|
||||
|
||||
end module math
|
||||
|
|
|
|||
Loading…
Reference in New Issue