Stream: helpdesk (published)

Topic: Is there a package with complex bessel functions?

Robbie Rosati (Sep 29 2023 at 23:04):

Neither SpecialFunctions.jl or Bessels.jl support complex order parameters. For my specific problem I need something like hankelh1(-im,1.0) to work. (oops, could an admin move this to the published helpdesk stream? Thanks)

Notification Bot (Sep 30 2023 at 06:40):

This topic was moved here from #helpdesk > Is there a package with complex bessel functions? by Mason Protter.

Robbie Rosati (Oct 04 2023 at 19:07):

Since I only need a few values of the order, I ended up just tabulating these with Mathematica. I'd love a native Julia solution, if I get enough free time to research implementations maybe I'll submit a PR to Bessels.jl . These complex-order Bessel functions seem fairly obscure, even gsl seems also to only support real orders. Apparently mpmath in Python does have them, so you could also use PyCall

Anyway, for posterity, here is the Mathematica code I used:

tab = Table[HankelH1[Sqrt[9/4 - 3], Exp@x], {x, -6, 10, 0.01}];
code = "[" <>
   Map[ToString[Re@#] <> " + im*" <> ToString[Im@#] <> ",\n" &, tab ] <>
    "]";
Export["hankel_table.jl", code, "Text",
 CharacterEncoding -> "Unicode"]

Robbie Rosati (Oct 06 2023 at 23:10):

A not particularly fast implementation (but faster than calling Python) is also possible through identities with HypergeometricFunctions.jl (I think I'll use this one in practice, to avoid needing to tabulate and interpolate). I'm sure this is not as numerically accurate as it could be.

import SpecialFunctions: gamma
import HypergeometricFunctions: pFq
besselj(nu,x) = (x/2)^nu / gamma(nu+1) * pFq((),(1+nu,),-(x/2)^2)
bessely(nu,x) = (besselj(nu,x)*cos(nu*π) - besselj(-nu,x))/sin(nu*π)
hankelh1(nu,x) = besselj(nu,x) + im*bessely(nu,x)
hankelh2(nu,x) = besselj(nu,x) - im*bessely(nu,x)

ederag (Oct 07 2023 at 07:54):

Thanks for sharing your solutions !
For the bessely line, you might be interested in sincospi ?

Robbie Rosati (Oct 09 2023 at 20:20):

Oh thanks! Yes that looks very relevant, I didn't know about that function. It's in Base too

Last updated: Nov 06 2024 at 04:40 UTC