modified_bessel_temme_algorithm Subroutine

      subroutine modified_bessel_temme_algorithm(v, x, ck, dk, ek)
         !! Implementation of the Temme's algorithm
         !! [N. M. Temme, J. Comput. Phys. 19 (1975) 324] to calculating
         !! modified Bessel functions of the second kind.
         !! This is a direct modernization of the "mbessk" subroutine
         !! in MOLSCAT:
         !! https://github.com/molscat/molscat/blob/master/source_code/rbessk.f
         !---------------------------------------------------------------------!
         real(dp), intent(in) :: v, x
         real(dp), intent(out) :: ck, dk, ek
         real(dp) :: a, b, c, d, e, f, g, h, p, q, s, sk, y, ak, ak1, ex, pi
         integer :: n, m, na, maxit
         real(dp), parameter :: eps = 1.0e-15_dp, xmin = 1.0_dp
         !---------------------------------------------------------------------!
         pi = acos(-1.0_dp)
         if (v < 0.0_dp .or. x <= 0.0_dp) then
            call write_error("modified_bessel_temme_algorithm: Invalid input"//&
               " values for v or x")
         endif

         na = int(v + 0.5_dp)
         a = v - na
         if (x < xmin) then
         ! Small x: Temme's series for small x
            b = x / 2.0_dp
            d = -log(b)
            e = a * d
            c = a * pi
            if (abs(c) < eps) then
             c = 1.0_dp
            else
             c = c / sin(c)
            endif
            if (abs(e) < eps) then
             s = 1.0_dp
            else
             s = sinh(e) / e
            endif
            e = exp(e)
            ! Compute the gamma function and its derivatives P and Q
            ! Replace RGAMMA(A,P,Q) with a modern equivalent if necessary
            g = e * rgamma(a, p, q)
            e = (e + 1.0_dp / e) / 2.0_dp
            f = c * (p * e + q * s * d)
            e = a * a
            p = 0.5_dp * g * c
            q = 0.5_dp / g
            c = 1.0_dp
            d = b * b
            ak = f
            ak1 = p
            do n = 1, maxit
             f = (f * n + p + q) / (n * n - e)
             c = c * d / n
             p = p / (n - a)
             q = q / (n + a)
             g = c * (p - n * f)
             h = c * f
             ak = ak + h
             ak1 = ak1 + g
             if (abs(h / ak) + abs(g / ak1) < eps) exit
            end do
            f = ak
            g = ak1 / b
            ex = 0.0_dp
         else
         ! Large x: Temme's PQ method for large x
            c = 0.25_dp - a * a
            g = 1.0_dp
            f = 0.0_dp
            e = x * cos(a * pi) / pi / eps
            do n = 1, maxit
             h = (2 * (n + x) * g - (n - 1 + c / n) * f) / (n + 1)
             f = g
             g = h
             if (h * n > e) exit
            end do

            p = f / g
            q = p
            b = x + x
            e = b - 2.0_dp
            do m = n, 1, -1
             p = (m - 1 + c / m) / (e + (m + 1) * (2.0_dp - p))
             q = p * (q + 1.0_dp)
            end do
            f = sqrt(pi / b) / (1.0_dp + q)
            g = f * (a + x + 0.5_dp - p) / x
            ex = x
         endif

         ! Upward recursion
         p = 0.0_dp
         if (na > 0) then
            y = 2.0_dp / x
            do n = 1, na
              h = y * (a + n) * g + f
              f = g
              g = h
              if (abs(f) > 4.0_dp) then
                p = p + 1.0_dp
                f = 0.0625_dp * f
                g = 0.0625_dp * g
              endif
            end do
         endif

         ck = f
         dk = (v / x) * f - g
         sk = sqrt(ck * ck + dk * dk)
         ck = ck / sk
         dk = dk / sk
         ek = log(sk) + p * log(16.0_dp) - ex
      end subroutine modified_bessel_temme_algorithm