Hankel functions overflowing for |z| < 10^{-1}

[joeydumont]

For low |z|, it seems that the real part of the Hankel functions, i.e. zbesh, overflow instead of going to zero. The original documentation indicates that tests are performed to limit overflows, but I need to confirm this.

My initial tests show that SciPy also suffers the same issue.

The code below shows the issue:

#include <complex_bessel.h>
#include <iostream>
#include <array>

using namespace std::complex_literals;
using namespace sp_bessel;

static constexpr int size = 20;
static constexpr double exp_min = -15;
static constexpr double exp_max = 1;
static constexpr double dx = (exp_max-exp_min)/size;
static const int nu = 4;

int main (int argc, char* argv[])
{
  std::array<double,size> exponent;
  std::array<double,size> z;

  for (auto i = 0; i < size; i++)
  {
    exponent[i] = exp_min + i*dx;
    z[i] = std::pow(10,exponent[i]);
    std::cout << "Hankel function: z[i]" << z[i] << ", " << hankelH2(nu, z[i]) << std::endl;
    std::cout << "Bessel function: z[i]" << z[i] << ", " << besselJ(nu, z[i])-1i*besselY(nu,z[i]) << std::endl;
    //std::cout << "Bessel function of the second kind: " << besselY(nu, z[i]) << std::endl;
  }

  //double z = 8.693214466689956e-11;


  return 0;
}

Thanks to @NiDingRte for reporting the issue.