Skip to content

Commit 3f4172f

Browse files
heltonmcheltonmc
committed
add complex besseli1 and besselj1
1 parent adce2ca commit 3f4172f

File tree

4 files changed

+126
-1
lines changed

4 files changed

+126
-1
lines changed

‎src/besseli.jl

Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -243,6 +243,125 @@ function besseli0(z::ComplexF32)
243243
return ComplexF32(r)
244244
end
245245

246+
function besseli1(z::ComplexF64)
247+
c = one(z)
248+
# shift phase to 0 < angle(z) < pi/2
249+
if real(z) < 0.0
250+
z = -z
251+
c = -c
252+
end
253+
if imag(z) < 0.0
254+
z = conj(z)
255+
check = true
256+
else
257+
check = false
258+
end
259+
260+
if abs(z) > 17.5
261+
zinv = 1 / z
262+
e = exp(z)
263+
s = sqrt(2 * z * π)
264+
p = evalpoly(zinv*zinv, (1.0, -0.1171875, -0.144195556640625, -0.6765925884246826, -6.883914268109947, -121.59789187653587, -3302.2722944808525, -127641.2726461746, -6.656367718817688e6, -4.502786003050393e8, -3.8338575207427895e10, -4.0118385991331978e12, -5.060568503314726e14, -7.572616461117957e16, -1.3262572853205555e19))
265+
p2 = zinv*evalpoly(zinv*zinv, (0.375, 0.1025390625, 0.2775764465332031, 1.993531733751297, 27.248827311268542, 603.8440767050702, 19718.37591223663, 890297.8767070678, 5.310411010968523e7, 4.043620325107754e9, 3.827011346598606e11, 4.406481417852279e13, 6.065091351222699e15, 9.83388387659068e17))
266+
r = e / s * (p - p2) - im * (p + p2) / (e * s)
267+
elseif abs(z) < 5.2
268+
r = z*evalpoly(z*z, (0.5, 0.0625, 0.0026041666666666665, 5.425347222222222e-5, 6.781684027777778e-7, 5.651403356481481e-9, 3.363930569334215e-11, 1.5017547184527747e-13, 5.214426105738801e-16, 1.4484516960385557e-18, 3.2919356728148996e-21, 6.234726653058522e-24, 9.991549123491221e-27, 1.3724655389411017e-29, 1.6338875463584545e-32, 1.7019661941233902e-35, 1.564307163716351e-38))
269+
else
270+
if angle(z) <= π / 4.4
271+
zz = z*z
272+
zz2 = zz*zz
273+
r = evalpoly(zz2, (0.5, 0.0026041666666666665, 6.781684027777778e-7, 3.363930569334215e-11, 5.214426105738801e-16, 3.2919356728148996e-21, 9.991549123491221e-27, 1.6338875463584545e-32, 1.564307163716351e-38, 9.342315267006072e-45, 3.658488121477942e-51, 9.781133223498595e-58, 1.845775442236299e-64, 2.5281824488224564e-71, 2.574012221626064e-78, 1.9883297967078112e-85, 1.1862954038963049e-92, 5.553068705606664e-100))
274+
r += zz*evalpoly(zz2, (0.0625, 5.425347222222222e-5, 5.651403356481481e-9, 1.5017547184527747e-13, 1.4484516960385557e-18, 6.234726653058522e-24, 1.3724655389411017e-29, 1.7019661941233902e-35, 1.2780287285264307e-41, 6.146260044082942e-48, 1.9797013644361157e-54, 4.4298610613671176e-61, 7.099136316293458e-68, 8.360391695841456e-75, 7.396586843753058e-82, 5.010911786057992e-89, 2.6432607038687723e-96))
275+
r *= z
276+
else
277+
zz = (-1.3144285719254842 + 20.054304662400146im, -1.2128746095534062 + 18.50489060934118im,
278+
-1.0886579984944276 + 16.60971135387317im, -0.9427507983505322 + 14.383597659587654im,
279+
-0.7876981565498402 + 12.017951489232377im, -0.6255235989509672 + 9.543645881424778im,
280+
-0.4528832294653219 + 6.9096628407010305im, -0.321046242530308 + 4.898219116608337im,
281+
16.905702128729647 + 10.872315095446137im, 4.119703326027316 + 2.6529313043347877im,
282+
-1.2665282079632723 + 19.323486333541947im, -0.3857883819775798 + 5.8859932845642575im,
283+
5.040563157078601 + 19.45771628582095im, 8.609436758553485 + 5.532949040120932im,
284+
0.694001756513151 + 20.08801537140881im, -0.5457083486368534 + 8.325900481870493im,
285+
2.790091405909159 + 4.028075216114001im, 12.479613762034912 + 15.75656181882421im,
286+
-1.0195312436724646 + 15.555040882512422im, 13.06890672928664 + 8.39887636916513im,
287+
5.1276393456482685 + 3.295333712441098im, 2.7253844460913528 + 19.914373693917753im,
288+
14.171480465360736 + 9.107457483790645im, 0.8038847904678548 + 4.833608304740307im
289+
)
290+
w = (0.013404709995345458 + 0.0im, -0.10348199668238728 + 0.03550804014042965im,
291+
-0.11679300367073918 + 0.16207579437860672im, -0.03568356238246015 - 0.16027459313035955im,
292+
-0.060434459186566966 - 0.02767438331782018im, -0.07086407749919574 + 0.06788442528741509im,
293+
-0.08650269086297906 - 0.029416238653716714im, -0.00415412382674645 - 0.00607269307464697im,
294+
-0.07587326685881253 + 0.19888070181762688im, -0.08774917184927954 - 0.014101104643344704im,
295+
0.007193185558007104 + 0.06995793660348812im, -0.02672800181750328 - 0.03807797543563465im,
296+
0.19776072630033434 + 0.05342203578372187im, 0.3489984378687633 - 0.29126891724262455im,
297+
-0.010279561959668174 + 0.09872927973051646im, -0.12150771925043644 - 0.012812175864887374im,
298+
-0.1299759980313551 - 0.07418447933244608im, -0.05894210892668888 + 0.24161474004043546im,
299+
-0.2648591701697078 - 0.07362611929106554im, -0.000770797654486471 + 0.00014803398133387874im,
300+
0.12194676645309753 - 0.25480936217584116im, 0.0718337239451382 + 0.20114630545134515im,
301+
0.5172999110471845 - 0.09099393187858433im, -0.023839427516594167 - 0.056052737377955054im
302+
)
303+
f = (-3.215770403317309 + 4.184450606321457im, 3.196349209039463 - 3.528571814265253im,
304+
-3.0032703696390066 + 0.8968536198145782im, 1.7024961674754864 + 2.2913462636360844im,
305+
2.0860054490489937 - 0.919904256366073im, 0.019539554711925467 - 1.3245943976605632im,
306+
-0.5518276420937166 - 0.23608328648174784im, 0.7317073656843707 + 0.7151746200203796im,
307+
0.3926351501115586 + 0.00413652715802052im, 0.372728884764365 + 0.01864352217896087im,
308+
-3.733141424945208 - 2.8413601177219308im, 0.9786880011963323 - 0.618277150144583im,
309+
0.3971643533953302 + 0.007251821636383407im, 0.38646933950897416 + 0.008351088369392122im,
310+
0.33879033458506147 + 0.08702715214401817im, 1.3930956061616984 + 0.6707847701397638im,
311+
0.38110598793307165 + 0.02730656109552859im, 0.39435250888241274 + 0.005949183211885836im,
312+
1.251416864948188 - 2.932439718874822im, 0.39076282374400617 + 0.005394185668824364im,
313+
0.3778378309823485 + 0.014609112728729376im, 0.39660674100748977 + 0.008337843914085801im,
314+
0.39140459884033496 + 0.004960218953458415im, 0.4210441290407154 + 0.10788322118417641im
315+
).*w
316+
317+
s1 = 0.0 + 0.0im
318+
s2 = 0.0 + 0.0im
319+
@fastmath for ind in eachindex(f)
320+
C = inv(z - zz[ind])
321+
s1 += C*f[ind]
322+
s2 += C*w[ind]
323+
end
324+
325+
r = s1 / (s2 * sqrt(z) * exp(-z))
326+
end
327+
end
328+
check && (r = conj(r))
329+
return r*c
330+
end
331+
332+
function besseli1(z::ComplexF32)
333+
z = ComplexF64(z)
334+
c = one(z)
335+
# shift phase to 0 < angle(z) < pi/2
336+
if real(z) < 0.0
337+
z = -z
338+
c = -c
339+
end
340+
if imag(z) < 0.0
341+
z = conj(z)
342+
check = true
343+
else
344+
check = false
345+
end
346+
347+
if abs(z) > 10.0
348+
zinv = 1 / z
349+
e = exp(z)
350+
s = sqrt(2 * z * π)
351+
p = evalpoly(zinv*zinv, (1.0, -0.1171875, -0.144195556640625, -0.6765925884246826, -6.883914268109947, -121.59789187653587, -3302.2722944808525))
352+
p2 = zinv*evalpoly(zinv*zinv, (0.375, 0.1025390625, 0.2775764465332031, 1.993531733751297, 27.248827311268542, 603.8440767050702, 19718.37591223663))
353+
r = e / s * (p - p2) - im * (p + p2) / (e * s)
354+
else
355+
zz = z*z
356+
zz2 = zz*zz
357+
r = evalpoly(zz2, (0.5, 0.0026041666666666665, 6.781684027777778e-7, 3.363930569334215e-11, 5.214426105738801e-16, 3.2919356728148996e-21, 9.991549123491221e-27, 1.6338875463584545e-32, 1.564307163716351e-38, 9.342315267006072e-45))
358+
r += zz*evalpoly(zz2, (0.0625, 5.425347222222222e-5, 5.651403356481481e-9, 1.5017547184527747e-13, 1.4484516960385557e-18, 6.234726653058522e-24, 1.3724655389411017e-29, 1.7019661941233902e-35, 1.2780287285264307e-41, 6.146260044082942e-48))
359+
r *= z
360+
end
361+
check && (r = conj(r))
362+
return ComplexF32(r*c)
363+
end
364+
246365
# Modified Bessel functions of the first kind of order nu
247366
# besseli(nu, x)
248367
#

‎src/besselj.jl

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -177,7 +177,9 @@ end
177177
##### Implementation for complex arguments
178178
#####
179179

180-
besselj0(z::T) where T <: Union{ComplexF32, ComplexF64} = besseli0(z/im)
180+
besselj0(z::T) where T <: Union{ComplexF32, ComplexF64} = besseli0(z/im)
181+
besselj1(z::T) where T <: Union{ComplexF32, ComplexF64} = besseli1(z/im) * im
182+
181183

182184
# Bessel functions of the first kind of order nu
183185
# besselj(nu, x)

‎test/besseli_test.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,8 @@ for x in [0.0, 0.01, 0.5, 1.0, 2.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14
7171
z = x*exp(im*a)
7272
@test isapprox(besseli0(z), SpecialFunctions.besseli(0, z), rtol=2e-14)
7373
@test isapprox(besseli0(ComplexF32(z)), ComplexF32(SpecialFunctions.besseli(0, ComplexF32(z))), rtol=1e-7)
74+
@test isapprox(besseli1(z), SpecialFunctions.besseli(1, z), rtol=2e-14)
75+
@test isapprox(besseli1(ComplexF32(z)), ComplexF32(SpecialFunctions.besseli(1, ComplexF32(z))), rtol=1e-7)
7476
end
7577

7678
# test for besseli

‎test/besselj_test.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -83,6 +83,8 @@ for x in [0.0, 0.01, 0.5, 1.0, 2.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14
8383
z = x*exp(im*a)
8484
@test isapprox(besselj0(z), SpecialFunctions.besselj(0, z), rtol=2e-14)
8585
@test isapprox(besselj0(ComplexF32(z)), ComplexF32(SpecialFunctions.besselj(0, ComplexF32(z))), rtol=1e-7)
86+
@test isapprox(besselj1(z), SpecialFunctions.besselj(1, z), rtol=2e-14)
87+
@test isapprox(besselj1(ComplexF32(z)), ComplexF32(SpecialFunctions.besselj(1, ComplexF32(z))), rtol=1e-7)
8688
end
8789

8890
## Tests for besselj

0 commit comments

Comments
 (0)