Is there a math library in Julia which can compute a modified Bessel function of the second kind with complex order and argument? I am only aware of `SpecialFunctions.jl`

which contains a method to compute a modified Bessel function of the second kind with a complex argument of real order.

# Bessel function of second with complex order

**stevengj**#3

I had a student working on a related problem at one point — my recollection is that we couldn’t find any free code for this, and very little literature. He developed some prototypes, but as I recall he ran into tricky accuracy issues with the second-kind Bessel functions for complex arguments that were close to integers. We also wanted the derivative with respect to the order, however, so maybe the problem was easier if you don’t need that. I’ll bug him to see whether he has any notes he can post.

If you want to use a package that is in development,

`add SpecialFunctions`

and

`add ArbNumerics`

(must use version 0.3.5+, just merged).

```
using ArbNumerics
import SpecialFunctions: besselj, bessely
for F in (:besselj, :bessely)
@eval function $F(order::Complex{T}, argument::Complex{T}) where {T}
result = $F(ArbComplex(order), ArbComplex(argument))
return Complex{T}(result)
end
end
```

then

```
order = 0.0 + 0.5im; argument = 1.0 + 1.0im;
besselj(order, argument)
# 0.6415000256096106 - 0.2518459064184279im
bessely(order, argument)
# 0.12578491928293956 + 0.17534238733104368im
```

Thank you all for your response. I finally ended up constructing my own function using `SpecialFunctions.jl`

. The lines of codes below seem to be performing very well.

```
using SpecialFunctions
# Define the modified Bessel Function of the second kind with complex order & argument
function CBesselK(nu, z)
Knu = exp(lgamma(nu)) * (z/2)^(-nu)
Kmu = exp(lgamma(-nu)) * (z/2)^nu
bound = 20
precision = 1e-10
s = 0
if (abs(z) < 10)
for i = 0:bound
var = i
facto = factorial(i)
r = 0.5 * Knu * ((z/2)^(2*var) * exp(-lgamma(1+var-nu)+lgamma(1-nu))/facto) + 0.5 * Kmu * ((z/2)^(2*var) * exp(-lgamma(1+var+nu)+lgamma(1+nu))/facto)
s+=r
if ((abs(real(r/s))<precision)&&(abs(imag(r/s))<precision))
break
end
end
else
for i = 0:bound
var = i
facto = factorial(i)
r = sqrt(pi/(2*z)) * exp(-z) * exp(lgamma(0.5+var+nu)+lgamma(0.5+var-nu)-lgamma(0.5+nu)-lgamma(0.5-nu)/facto) * (-1/2*z)^var
s+=r
if ((abs(real(r/s))<precision)&&(abs(imag(r/s))<precision))
break
end
end
end
return s
end
```

Comparing the performance with **Mathematica**:

```
# CBesselK
order = exp(1-3im)
arg = (4+6im)*exp(2*im*pi)
CBesselK(order, arg)
0.011090297922744 - 0.005887726316203im
# Mathematica
N[BesselK[Exp[1 - 3 I], (4 + 6 I)*Exp[2*I*Pi]], 15]
0.011090297922749 - 0.005887726316251 I
```

Cheers!