anchor.diff:           cdexp((0,1)*thetax*m)*spherical_harmonic_(l,m,thetax,phix)
   
     end function G_x
  +
  +  double complex function G_harmonic_x(x,k,l,m)
  +    double precision, intent(in) :: x(:)
  +    double complex, intent(in) :: k
  +    integer, intent(in) :: l,m
  +    double precision :: rx, phix, thetax
  +
  +    rx = norm(x)
  +    phix = dacos(x(3)/rx)
  +    thetax = datan2(x(2),x(1))
  +
  +    G_harmonic_x = cdexp((0,1)*thetax*m)*spherical_harmonic_(l,m,thetax,phix)
  +
  +  end function G_harmonic_x
   
     ! formula in [[file:~/downloads/pub/papers/epstein2012convergence.pdf][epstein]]
   
  
positive.code: 
    end function G_y
  
    double complex function G_x(x,k,l,m)
      double precision, intent(in) :: x(:)
      double complex, intent(in) :: k
  
negative.code:     if (nm /= n) print *, "spherical_bessel_jn: sphj didn't converge"
      spherical_bessel_yn = sj(n)
    end function spherical_bessel_yn
  
    function aspherical_hankel_n(n,x)
      double precision, intent(in) :: x
      double complex, dimension(n+1) :: aspherical_hankel_n
      integer, intent(in) :: n
      aspherical_hankel_n(1:n+1) = aspherical_bessel_jn(n,x) + (0,1)*aspherical_bessel_yn(n,x)
    end function aspherical_hankel_n
  
    double complex function spherical_hankel_n(n,x)
      double precision, intent(in) :: x
  