anchor.diff:      res = resReal + i_*resImag
   end function dairybi
   
  +! rotate Hankel functions of negative order:
  +function Hrotate(hkind, order) result(factor)
  +    integer, intent(in) :: hkind
  +    real(dp), intent(in) :: order
  +    complex(dp) :: factor
  +
  +    ! NIST Handbook, 10.4.6; note the signs in the exponents
  +    ! (NIST gives formulas for order -v with positive v, our 'order' here is negative)
  +
  +    if(order < 0.0_dp) then
  +        if(hkind == 1) then
  +            factor = exp(-pi*order*i_)
  +        elseif(hkind == 2) then
  +            factor = exp(pi*order*i_)
  +        endif
  +    else
  +        factor = (1.0_dp, 0.0_dp)
  +    endif
  +
  +end function Hrotate
  +
  +! 'reflect' negative integer order Bessel functions if needed
  +function reflect(order) result(res)
  +    integer, intent(in) :: order
  +    integer :: res
  +
  +    if(order < 0) then
  +        res = (-1)**order
  +    else
  +        res = 1
  +    endif
  +end function
  +
   end module
  
positive.code: 
  
  function hankel1v_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=1
      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:     ! convenience function that calls Bessel functions of integer order
      ! (intrinsic to F2008)
  
      if(abs(order) == 0) then
          res = bessel_j0(x) - i_*bessel_y0(x)
      elseif(abs(order) == 1) then
          res = reflect(order)*(bessel_j1(x) - i_*bessel_y1(x))
      else
          res = reflect(order)*(bessel_jn(order, x) - i_*bessel_yn(order, x))
      endif
  end function hankel2n_real
  
  
  function hankel2v_real(order, x) result(res)
      real(dp), intent(in) :: order, x
      complex(dp) :: res
  
  