anchor.diff:  
   contains
   
  -real(dp) function f(x)
  -real(dp), intent(in) :: x
  -f = spherical_bessel_jn(n, x)
  +    real(dp) function f(x)
  +        real(dp), intent(in) :: x
  +        f = spherical_bessel_jn(n, x)
  +    end function f
  +
   end function
   
  -end function
   
   ! The SPHJ, SPHY, MSTA1, MSTA2 routines below are taken from SciPy's specfun.f.
   ! Authors: Shanjie Zhang and Jianming Jin
   ! Copyrighted but permission granted to use code in programs.
  -
           SUBROUTINE SPHJ(N,X,NM,SJ,DJ)
   !       =======================================================
   !       Purpose: Compute spherical Bessel functions jn(x) and
  
positive.code: 
  end function
  
  real(dp) function spherical_bessel_jn(n, x) result(r)
  integer, intent(in) :: n
  real(dp), intent(in) :: x
  
negative.code: if (nm /= n) call stop_error("spherical_bessel_jn: sphj didn't converge")
  r = sj(n)
  end function
  
  
  function spherical_bessel_jn_zeros(nmax, nzeros, eps) result(zeros)
  ! Calculates 'nzeros' zeros of spherical_bessel_jn() for all n=0, 1, ..., nmax
  