-
Notifications
You must be signed in to change notification settings - Fork 11
/
hankel.jl
207 lines (180 loc) · 7.45 KB
/
hankel.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
# Hankel functions
#
# besseljy(nu, x)
# besselh(nu, x), hankelh1(nu, x), hankelh2(nu, x)
#
# A numerical routine to compute both Bessel functions of the first J_{ν}(x) and second kind Y_{ν}(x)
# for real orders and arguments of positive or negative value. Please see notes in src/besselj.jl and src/bessely.jl
# for details on implementation as most of the routine is similar. A key difference is when the methods to compute bessely
# don't also give besselj. We then rely on a continued fraction approach to compute J_{ν}(x)/J_{ν-1}(x) to get J_{ν}(x)
# from Y_{ν}(x) and Y_{ν-1}(x)[1]. The continued fraction approach is more quickly converging when nu and x are small in magnitude.
# When x and nu are large and nu = x + ϵ, we fall back to computing J_{ν}(x) and Y_{ν}(x) separately as this was found to be more efficient.
#
# [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."
# Computer physics communications 76.3 (1993): 381-388.
#
#####
##### Generic routine for `besseljy`
#####
"""
besseljy(nu, x::T) where T <: Float64
Return both Bessel functions of the first ``J_{nu}(x)`` and
second ``Y_{nu}(x)`` kind. This method will be faster than calling (besselj(nu, x), bessely(nu, x)) separately,
unless nu is slightly larger than x.
Results may be slightly different than individual functions in some domains due to using different algorithms.
"""
function besseljy(nu::Real, x::T) where T
isinteger(nu) && return besseljy(Int(nu), x)
abs_nu = abs(nu)
abs_x = abs(x)
Jnu, Ynu = besseljy_positive_args(abs_nu, abs_x)
if nu >= zero(T)
if x >= zero(T)
return Jnu, Ynu
else
return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
#return Ynu * cispi(-nu) + 2im * besselj_positive_args(abs_nu, abs_x) * cospi(abs_nu)
end
else
spi, cpi = sincospi(abs_nu)
out = Jnu * cpi - Ynu * spi
if x >= zero(T)
return out, Ynu * cpi + Jnu * spi
else
return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
#return cpi * (Ynu * cispi(nu) + 2im * Jnu * cpi) + Jnu * spi * cispi(abs_nu)
end
end
end
function besseljy(nu::Integer, x::T) where T
abs_nu = abs(nu)
abs_x = abs(x)
sg = iseven(abs_nu) ? 1 : -1
Jnu, Ynu = besseljy_positive_args(abs_nu, abs_x)
if nu >= zero(T)
if x >= zero(T)
return Jnu, Ynu
else
return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
#return Jnu * sg, Ynu * sg + 2im * sg * besselj_positive_args(abs_nu, abs_x)
end
else
if x >= zero(T)
return Jnu * sg, Ynu * sg
else
return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
#spi, cpi = sincospi(abs_nu)
#return (cpi*Jnu - spi*Ynu) * sg, Ynu + 2im * besselj_positive_args(abs_nu, abs_x)
end
end
end
#####
##### `besseljy` for positive arguments and orders
#####
function besseljy_positive_args(nu::Real, x::T) where T
nu == 0 && return (besselj0(x), bessely0(x))
nu == 1 && return (besselj1(x), bessely1(x))
iszero(x) && return (besselj_positive_args(nu, x), bessely_positive_args(nu, x))
# x < ~nu branch see src/U_polynomials.jl
besseljy_debye_cutoff(nu, x) && return besseljy_debye(nu, x)
# large argument branch see src/asymptotics.jl
besseljy_large_argument_cutoff(nu, x) && return besseljy_large_argument(nu, x)
# x > ~nu branch see src/U_polynomials.jl on computing Hankel function
if hankel_debye_cutoff(nu, x)
H = hankel_debye(nu, x)
return real(H), imag(H)
end
# use forward recurrence if nu is an integer up until the continued fraction becomes inefficient
if isinteger(nu) && nu < 150
Y0 = bessely0(x)
Y1 = bessely1(x)
Ynm1, Yn = besselj_up_recurrence(x, Y1, Y0, 1, nu-1)
ratio_Jv_Jvm1 = besselj_ratio_jnu_jnum1(nu, x)
Jn = 2 / (π*x * (Ynm1 - Yn / ratio_Jv_Jvm1))
return Jn, Yn
end
# use power series for small x and for when nu > x
if bessely_series_cutoff(nu, x) && besselj_series_cutoff(nu, x)
Yn, Jn = bessely_power_series(nu, x)
return Jn, Yn
end
# for x ∈ (6, 19) we use Chebyshev approximation and forward recurrence
if besseljy_chebyshev_cutoff(nu, x)
Yn, Ynp1 = bessely_chebyshev(nu, x)
ratio_Jvp1_Jv = besselj_ratio_jnu_jnum1(nu+1, x)
Jnp1 = 2 / (π*x * (Yn - Ynp1 / ratio_Jvp1_Jv))
return Jnp1 / ratio_Jvp1_Jv, Yn
end
# at this point x > 19.0 (for Float64) and fairly close to nu
# shift nu down and use the debye expansion for Hankel function (valid x > nu) then use forward recurrence
nu_shift = ceil(nu) - floor(Int, -1.5 + x + Base.Math._approx_cbrt(-411*x)) + 2
v2 = maximum((nu - nu_shift, modf(nu)[1] + 1))
Hnu = hankel_debye(v2, x)
Hnum1 = hankel_debye(v2 - 1, x)
# forward recurrence is stable for Hankel when x >= nu
if x >= nu
H = besselj_up_recurrence(x, Hnu, Hnum1, v2, nu)[1]
return real(H), imag(H)
else
# At this point besselj can not be calculated with forward recurrence
# We could calculate it from bessely using the continued fraction approach like the following
# Yn, Ynp1 = besselj_up_recurrence(x, imag(Hnu), imag(Hnum1), v2, nu)
# ratio_Jvp1_Jv = besselj_ratio_jnu_jnum1(nu+1, x)
# Jnp1 = 2 / (π*x * (Yn - Ynp1 / ratio_Jvp1_Jv))
# return Jnp1 / ratio_Jvp1_Jv, Yn
# However, the continued fraction approach is slowly converging for large arguments
# We will fall back to computing besselj separately instead
return besselj_positive_args(nu, x), besselj_up_recurrence(x, imag(Hnu), imag(Hnum1), v2, nu)[1]
end
end
#####
##### Continued fraction for J_{ν}(x)/J_{ν-1}(x)
#####
# implements continued fraction to compute ratio of J_{ν}(x)/J_{ν-1}(x)
# using equation 22 and 24 of [1]
# in general faster converging for small magnitudes of x and nu and nu >> x
function besselj_ratio_jnu_jnum1(n, x::T) where T
MaxIter = 5000
xinv = inv(x)
xinv2 = 2 * xinv
d = x / (n + n)
a = d
h = a
b = muladd(2, n, 2) * xinv
for _ in 1:MaxIter
d = inv(b - d)
a *= muladd(b, d, -1)
h = h + a
b = b + xinv2
abs(a / h) <= eps(T) && break
end
return h
end
#####
##### Hankel functions
#####
"""
besselh(nu, [k=1,] x)
Bessel function of the third kind of order `nu` (the Hankel function). `k` is either 1 or 2,
selecting [`hankelh1`](@ref) or [`hankelh2`](@ref), respectively.
"""
function besselh(nu::Real, k::Integer, x)
Jn, Yn = besseljy(nu, x)
if k == 1
return complex(Jn, Yn)
elseif k == 2
return complex(Jn, -Yn)
else
throw(ArgumentError("k must be 1 or 2"))
end
end
"""
hankelh1(nu, x)
Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x)``.
"""
hankelh1(nu, x) = besselh(nu, 1, x)
"""
hankelh2(nu, x)
Bessel function of the third kind of order `nu`, ``H^{(2)}_\\nu(x)``.
"""
hankelh2(nu, x) = besselh(nu, 2, x)