!=======================================================================
!  Independent Spherical Bessel Function Library (Double Precision Complex)
!=======================================================================
module bessel_lib
  implicit none
  private
  public :: sph_j, sph_n, sph_h1, sph_h2
  public :: sph_dj, sph_dn, sph_dh1, sph_dh2

  ! Mathematical constants
  real*8, parameter     :: pi = 3.14159265358979323846d0
  complex*16, parameter :: ui = (0.0d0, 1.0d0)

contains

  !---------------------------------------------------------------------
  !  Spherical Bessel Function j_l(z) 
  !  Implemented via Backward Recurrence (Miller's Algorithm)
  !---------------------------------------------------------------------
  recursive function sph_j(l, z) result(res)
    integer, intent(in) :: l
    complex*16, intent(in) :: z
    complex*16 :: res, j0
    integer :: i, lmax
    complex*16, allocatable :: j_arr(:)

    if (abs(z) < 1d-15) then
       if (l == 0) then
          res = 1.0d0
       else
          res = 0.0d0
       end if
       return
    end if

    lmax = max(l, int(abs(z))) + 20
    allocate(j_arr(0:lmax))
    
    j_arr(lmax) = (0.0d0, 0.0d0)
    j_arr(lmax-1) = (1.0d-30, 0.0d0) 
    
    do i = lmax-1, 1, -1
       j_arr(i-1) = dble(2*i + 1)/z * j_arr(i) - j_arr(i+1)
    end do
    
    j0 = sin(z)/z
    res = j_arr(l) * (j0 / j_arr(0))
    deallocate(j_arr)
  end function sph_j

  !---------------------------------------------------------------------
  !  Spherical Neumann Function n_l(z)
  !---------------------------------------------------------------------
  recursive function sph_n(l, z) result(res)
    integer, intent(in) :: l
    complex*16, intent(in) :: z
    complex*16 :: res, n0, n1, n_curr
    integer :: i

    if (abs(z) < 1d-15) then
       res = -fact2_val(2*l-1) / (z**(l+1))
       return
    end if

    n0 = -cos(z)/z
    if (l == 0) then
       res = n0
       return
    end if
    
    n1 = n0/z - sin(z)/z
    if (l == 1) then
       res = n1
       return
    end if
    
    do i = 1, l-1
       n_curr = dble(2*i + 1)/z * n1 - n0
       n0 = n1
       n1 = n_curr
    end do
    res = n1
  end function sph_n

  !---------------------------------------------------------------------
  !  Spherical Hankel Function h_l^(1)(z)
  !---------------------------------------------------------------------
  function sph_h1(l, z) result(res)
    integer, intent(in) :: l
    complex*16, intent(in) :: z
    complex*16 :: res, sum_val, term
    integer :: k

    if (abs(z) < 1d-15) then
       res = -ui * fact2_val(2*l-1) / (z**(l+1))
       return
    end if

    sum_val = 1.0d0
    term = 1.0d0
    do k = 1, l
       term = term * dble((l + k) * (l - k + 1)) / (dble(k) * (-2.0d0 * ui * z))
       sum_val = sum_val + term
    end do
    
    res = sum_val * exp(ui * (z - dble(l+1)*pi/2.0d0)) / z
  end function sph_h1

  !---------------------------------------------------------------------
  !  Spherical Hankel Function h_l^(2)(z)
  !---------------------------------------------------------------------
  function sph_h2(l, z) result(res)
    integer, intent(in) :: l
    complex*16, intent(in) :: z
    complex*16 :: res, sum_val, term
    integer :: k

    if (abs(z) < 1d-15) then
       res = ui * fact2_val(2*l-1) / (z**(l+1))
       return
    end if

    sum_val = 1.0d0
    term = 1.0d0
    do k = 1, l
       term = term * dble((l + k) * (l - k + 1)) / (dble(k) * (2.0d0 * ui * z))
       sum_val = sum_val + term
    end do
    
    res = sum_val * exp(-ui * (z - dble(l+1)*pi/2.0d0)) / z
  end function sph_h2

  !---------------------------------------------------------------------
  !  Derivatives
  !---------------------------------------------------------------------
  function sph_dj(l, z) result(res)
    integer, intent(in) :: l
    complex*16, intent(in) :: z
    complex*16 :: res
    if (l == 0) then
       res = -sph_j(1, z)
    else
       res = sph_j(l-1, z) - dble(l+1)/z * sph_j(l, z)
    end if
  end function sph_dj

  function sph_dn(l, z) result(res)
    integer, intent(in) :: l
    complex*16, intent(in) :: z
    complex*16 :: res
    if (l == 0) then
       res = -sph_n(1, z)
    else
       res = sph_n(l-1, z) - dble(l+1)/z * sph_n(l, z)
    end if
  end function sph_dn

  function sph_dh1(l, z) result(res)
    integer, intent(in) :: l
    complex*16, intent(in) :: z
    complex*16 :: res
    if (l == 0) then
       res = -sph_h1(1, z)
    else
       res = sph_h1(l-1, z) - dble(l+1)/z * sph_h1(l, z)
    end if
  end function sph_dh1

  function sph_dh2(l, z) result(res)
    integer, intent(in) :: l
    complex*16, intent(in) :: z
    complex*16 :: res
    if (l == 0) then
       res = -sph_h2(1, z)
    else
       res = sph_h2(l-1, z) - dble(l+1)/z * sph_h2(l, z)
    end if
  end function sph_dh2

  !---------------------------------------------------------------------
  !  Internal utility
  !---------------------------------------------------------------------
  function fact2_val(n) result(res)
    integer, intent(in) :: n
    real*8 :: res
    integer :: i
    res = 1.0d0
    do i = n, 1, -2
       res = res * dble(i)
    end do
  end function fact2_val

end module bessel_lib