anchor.diff:      if (n .eq. m) dn = 1
     end function dn
   
  +  function aspherical_bessel_jn(n, x)
  +    integer, intent(in) :: n
  +    double precision, intent(in) :: x
  +    double precision, dimension(n+1) :: aspherical_bessel_jn
  +    integer :: nm
  +    double precision :: sj(0:n), dj(0:n)
  +    call sphj(n, x, nm, sj, dj)
  +    if (nm /= n) print *, "spherical_bessel_jn: sphj didn't converge"
  +    aspherical_bessel_jn(1:n+1) = sj(0:n)
  +  end function aspherical_bessel_jn
  +
  +  function aspherical_bessel_yn(n, x)
  +    integer, intent(in) :: n
  +    double precision, intent(in) :: x
  +    double precision, dimension(n+1) :: aspherical_bessel_yn
  +    integer :: nm
  +    double precision :: sj(0:n), dj(0:n)
  +    call sphy(n, x, nm, sj, dj)
  +    if (nm /= n) print *, "spherical_bessel_jn: sphj didn't converge"
  +    aspherical_bessel_yn(1:n+1) = sj(0:n)
  +  end function aspherical_bessel_yn
  +
     double precision function spherical_bessel_jn(n, x)
       integer, intent(in) :: n
       double precision, intent(in) :: x
  
positive.code:          cdexp((0,1)*thetax*m)*spherical_harmonic_(l,m,thetax,phix)
  
    end function G_x
  
    ! formula in [[file:~/downloads/pub/papers/epstein2012convergence.pdf][epstein]]
  
  
negative.code: 
    end function G_y
  
    double complex function G_harmonic_y(y,k,l,m)
      double precision, intent(in) :: y(:)
      double complex, intent(in) :: k
      integer, intent(in) :: l,m
      double precision :: ry, phiy, thetay
  
      ry = norm(y)
      phiy = dacos(y(3)/ry)
      thetay = datan2(y(2),y(1))
  
      G_harmonic_y = (0,1)*realpart(k) * &
           cdexp((0,1)*(-thetay)*m)*spherical_harmonic_(l,m,thetay,phiy)
  
    end function G_harmonic_y
  
    double complex function G_x(x,k,l,m)
      double precision, intent(in) :: x(:)
      double complex, intent(in) :: k
  