From 230ffe8f35c10d8d4e265a3b84a47dac6bb978c8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 14 Feb 2018 13:03:50 +0100 Subject: [PATCH] Bessel functions for Gaussian scatter --- src/math.f90 | 131 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 131 insertions(+) diff --git a/src/math.f90 b/src/math.f90 index ee228834a..e88a5ed03 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -2821,5 +2821,136 @@ real(pReal) pure function math_limit(a, left, right) math_limit = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_limit, 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