anchor.diff:  use optimize, only: bisect
   use amos
   implicit none
  +
   private
   public bessel_jn_zeros, spherical_bessel_jn, spherical_bessel_jn_zeros, &
          besselj, bessely, hankel1, hankel2, besseli, besselk, &  ! note the missing underscores in Bessel functions
  
positive.code: 
  
  function hankel2v_complex(order, z) result(res)
      implicit none
      real(dp), intent(in) :: order
      complex(dp), intent(in) :: z
      complex(dp) :: res
  
      integer, parameter :: scaling=1, length=1, hankelkind=2
      real(dp) :: resReal(length), resImag(length)
      integer :: underflows, ierr
  
      call zbesh(real(z), aimag(z), order, scaling, hankelkind, length, resReal, resImag, underflows, ierr)
      if(ierr > 0) then
          if(ierr == 1) then
              print *, "IERR=1, INPUT ERROR - NO COMPUTATION"
  
negative.code:     ! here, we do rely in an informed user: if the user wanted Bessel functions
      ! of integer order n, he/she should call bessel_jn() (as of F2008) directly.
      ! Hence, here we just convert the argument to a complex number and call the
      ! AMOS routine
      ! TODO: on the long run, wrap the cephes routines!
  
      resCmplx = besseljv_complex(order, cmplx(x, 0.0_dp, kind=dp))  ! negative order is handled in the routine called
      res = real(resCmplx)  ! Bessel functions of real arument are real
  end function besseljv_real
  
  
  function besseljv_complex(order, z) result(res)
      real(dp), intent(in) :: order
      complex(dp), intent(in) :: z
      complex(dp) :: res
  
      integer, parameter :: scaling=1, length=1
      real(dp) :: resReal(length), resImag(length), callorder
      integer :: underflows, ierr
  
      ! make sure AMOS zbesj() is called with postive order
      callorder = abs(order)
  
      call zbesj(real(z), aimag(z), callorder, scaling, length, resReal, resImag, underflows, ierr)
      if(ierr > 0) then
          if(ierr == 1) then
              print *, "IERR=1, INPUT ERROR - NO COMPUTATION"
  