|
| 1 | +# Spherical Bessel functions |
| 2 | +# |
| 3 | +# sphericalbesselj(nu, x), sphericalbessely(nu, x) |
| 4 | +# |
| 5 | +# A numerical routine to compute the spherical bessel functions of the first and second kind. |
| 6 | +# For small arguments, the power series series is used for sphericalbesselj. If nu is not a big integer |
| 7 | +# then forward recurrence is used if x >= nu. If x < nu, forward recurrence for sphericalbessely is used |
| 8 | +# and then a continued fraction and wronskian is used to compute sphericalbesselj [1]. For all other values, |
| 9 | +# we directly call besselj and bessely routines. |
| 10 | +# |
| 11 | +# [1] Ratis, Yu L., and P. Fernández de Córdoba. "A code to calculate (high order) Bessel functions based on the continued fractions method." |
| 12 | +# Computer physics communications 76.3 (1993): 381-388. |
| 13 | +# |
| 14 | + |
| 15 | +##### |
| 16 | +##### Generic routine for `sphericalbesselj` |
| 17 | +##### |
| 18 | + |
| 19 | +""" |
| 20 | + sphericalbesselj(nu, x) |
| 21 | +
|
| 22 | +Spherical bessel function of the first kind of order `nu`, ``j_ν(x)``. This is the non-singular |
| 23 | +solution to the radial part of the Helmholz equation in spherical coordinates. |
| 24 | +""" |
| 25 | +sphericalbesselj(nu::Real, x::Real) = _sphericalbesselj(nu, float(x)) |
| 26 | + |
| 27 | +_sphericalbesselj(nu, x::Float16) = Float16(_sphericalbesselj(nu, Float32(x))) |
| 28 | + |
| 29 | +function _sphericalbesselj(nu::Real, x::T) where T <: Union{Float32, Float64} |
| 30 | + x < zero(T) && return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported")) |
| 31 | + if ~isfinite(x) |
| 32 | + isnan(x) && return x |
| 33 | + isinf(x) && return zero(x) |
| 34 | + end |
| 35 | + return nu < zero(T) ? sphericalbesselj_generic(nu, x) : sphericalbesselj_positive_args(nu, x) |
| 36 | +end |
| 37 | + |
| 38 | +sphericalbesselj_generic(nu, x::T) where T = SQPIO2(T) * besselj(nu + one(T)/2, x) / sqrt(x) |
| 39 | + |
| 40 | +##### |
| 41 | +##### Positive arguments for `sphericalbesselj` |
| 42 | +##### |
| 43 | + |
| 44 | +function sphericalbesselj_positive_args(nu::Real, x::T) where T |
| 45 | + isinteger(nu) && return sphericalbesselj_positive_args(Int(nu), x) |
| 46 | + return sphericalbesselj_small_args_cutoff(nu, x) ? sphericalbesselj_small_args(nu, x) : sphericalbesselj_generic(nu, x) |
| 47 | +end |
| 48 | + |
| 49 | +function sphericalbesselj_positive_args(nu::Integer, x::T) where T |
| 50 | + sphericalbesselj_small_args_cutoff(nu, x) && return sphericalbesselj_small_args(nu, x) |
| 51 | + return (x >= nu && nu < 250) || (x < nu && nu < 60) ? sphericalbesselj_recurrence(nu, x) : sphericalbesselj_generic(nu, x) |
| 52 | +end |
| 53 | + |
| 54 | +##### |
| 55 | +##### Power series expansion for small arguments |
| 56 | +##### |
| 57 | + |
| 58 | +sphericalbesselj_small_args_cutoff(nu, x::T) where T = x^2 / (4*nu + 110) < eps(T) |
| 59 | + |
| 60 | +function sphericalbesselj_small_args(nu, x::T) where T |
| 61 | + iszero(x) && return iszero(nu) ? one(T) : zero(T) |
| 62 | + x2 = x^2 / 4 |
| 63 | + coef = evalpoly(x2, (1, -inv(T(3)/2 + nu), -inv(5 + nu), -inv(T(21)/2 + nu), -inv(18 + nu))) |
| 64 | + a = SQPIO2(T) / (gamma(T(3)/2 + nu) * 2^(nu + T(1)/2)) |
| 65 | + return x^nu * a * coef |
| 66 | +end |
| 67 | + |
| 68 | +##### |
| 69 | +##### Integer recurrence and/or wronskian with continued fraction |
| 70 | +##### |
| 71 | + |
| 72 | +# very accurate approach however need to consider some performance issues |
| 73 | +# if recurrence is stable (x>=nu) can generate very fast up to orders around 250 |
| 74 | +# for larger orders it is more efficient to use other expansions |
| 75 | +# if (x<nu) we can use forward recurrence from sphericalbessely_recurrence and |
| 76 | +# then use a continued fraction approach. However, for largish orders (>60) the |
| 77 | +# continued fraction is slower converging and more efficient to use other methods |
| 78 | +function sphericalbesselj_recurrence(nu::Integer, x::T) where T |
| 79 | + if x >= nu |
| 80 | + # forward recurrence if stable |
| 81 | + xinv = inv(x) |
| 82 | + s, c = sincos(x) |
| 83 | + sJ0 = s * xinv |
| 84 | + sJ1 = (sJ0 - c) * xinv |
| 85 | + |
| 86 | + nu_start = one(T) |
| 87 | + while nu_start < nu + 0.5 |
| 88 | + sJ0, sJ1 = sJ1, muladd((2*nu_start + 1)*xinv, sJ1, -sJ0) |
| 89 | + nu_start += 1 |
| 90 | + end |
| 91 | + return sJ0 |
| 92 | + elseif x < nu |
| 93 | + # compute sphericalbessely with forward recurrence and use continued fraction |
| 94 | + sYnm1, sYn = sphericalbessely_forward_recurrence(nu, x) |
| 95 | + H = besselj_ratio_jnu_jnum1(nu + T(3)/2, x) |
| 96 | + return inv(x^2 * (H*sYnm1 - sYn)) |
| 97 | + end |
| 98 | +end |
| 99 | + |
| 100 | +##### |
| 101 | +##### Generic routine for `sphericalbessely` |
| 102 | +##### |
| 103 | + |
| 104 | +""" |
| 105 | + sphericalbessely(nu, x) |
| 106 | +
|
| 107 | +Spherical bessel function of the second kind at order `nu`, ``y_ν(x)``. This is the singular |
| 108 | +solution to the radial part of the Helmholz equation in spherical coordinates. Sometimes |
| 109 | +known as a spherical Neumann function. |
| 110 | +""" |
| 111 | +sphericalbessely(nu::Real, x::Real) = _sphericalbessely(nu, float(x)) |
| 112 | + |
| 113 | +_sphericalbessely(nu, x::Float16) = Float16(_sphericalbessely(nu, Float32(x))) |
| 114 | + |
| 115 | +function _sphericalbessely(nu::Real, x::T) where T <: Union{Float32, Float64} |
| 116 | + x < zero(T) && return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported")) |
| 117 | + if ~isfinite(x) |
| 118 | + isnan(x) && return x |
| 119 | + isinf(x) && return zero(x) |
| 120 | + end |
| 121 | + return nu < zero(T) ? sphericalbessely_generic(nu, x) : sphericalbessely_positive_args(nu, x) |
| 122 | +end |
| 123 | + |
| 124 | +sphericalbessely_generic(nu, x::T) where T = SQPIO2(T) * bessely(nu + one(T)/2, x) / sqrt(x) |
| 125 | + |
| 126 | +##### |
| 127 | +##### Positive arguments for `sphericalbesselj` |
| 128 | +##### |
| 129 | + |
| 130 | +sphericalbessely_positive_args(nu::Real, x) = isinteger(nu) ? sphericalbessely_positive_args(Int(nu), x) : sphericalbessely_generic(nu, x) |
| 131 | + |
| 132 | +function sphericalbessely_positive_args(nu::Integer, x::T) where T |
| 133 | + if besseljy_debye_cutoff(nu, x) |
| 134 | + # for very large orders use expansion nu >> x to avoid overflow in recurrence |
| 135 | + return SQPIO2(T) * besseljy_debye(nu + one(T)/2, x)[2] / sqrt(x) |
| 136 | + elseif nu < 250 |
| 137 | + return sphericalbessely_forward_recurrence(nu, x)[1] |
| 138 | + else |
| 139 | + return sphericalbessely_generic(nu, x) |
| 140 | + end |
| 141 | +end |
| 142 | + |
| 143 | +##### |
| 144 | +##### Integer recurrence |
| 145 | +##### |
| 146 | + |
| 147 | +function sphericalbessely_forward_recurrence(nu::Integer, x::T) where T |
| 148 | + xinv = inv(x) |
| 149 | + s, c = sincos(x) |
| 150 | + sY0 = -c * xinv |
| 151 | + sY1 = xinv * (sY0 - s) |
| 152 | + |
| 153 | + nu_start = one(T) |
| 154 | + while nu_start < nu + 0.5 |
| 155 | + sY0, sY1 = sY1, muladd((2*nu_start + 1)*xinv, sY1, -sY0) |
| 156 | + nu_start += 1 |
| 157 | + end |
| 158 | + # need to check if NaN resulted during loop |
| 159 | + # this could happen if x is very small and nu is large which eventually results in overflow (-> -Inf) |
| 160 | + # NaN inputs were checked in top level function so if sY0 is NaN we should return -infinity |
| 161 | + return isnan(sY0) ? (-T(Inf), -T(Inf)) : (sY0, sY1) |
| 162 | +end |
0 commit comments