module m_bessel
!---------------------------------------------------------------------------



  ! --- public routines
  public :: BesselJn

  
  private
  ! --- all is privat by default


  



  
contains


  ! -----------------------------------------------------------------------  
  function BesselJn( order, arg )  result(res)
  ! -----------------------------------------------------------------------  
    use m_precision, only : wp, sp, dp
    use m_constants, only : pi, imu
    !
    implicit none
    integer,     intent(in) :: order
    complex(wp), intent(in) :: arg
    complex(wp)             :: res
    ! --- local
    integer     :: no_underflow, kode, ierr
    real(wp)    :: rorder
    real(wp)                  :: arg_re, arg_im
    real(wp),    dimension(1) :: besselJ_re, besselJ_im
    complex(wp), dimension(1) :: besselJ
    
    ! --- real absolute order
    rorder = real(order,wp)
    !
    ! --- Do scaling of the Bessel function
    kode = 2
    !
    ! --- argument and its real and imag part.
    arg_re  =  real(     arg,wp)
    arg_im  =  real(-imu*arg,wp)
    !
    select case(wp) 
    case(sp)
       ! --- single precission
 !!!!!!!      call cbesj( arg_re, arg_im, abs(rorder), kode, size(BesselJ,1)
 !!!!!!!           &, besselJ_re, besselJ_im, no_underflow, ierr ) 
       !
    case(dp)
       ! --- double precission
       call zbesj( arg_re, arg_im, abs(rorder), kode, size(BesselJ,1)&
            &, besselJ_re, besselJ_im, no_underflow, ierr ) 
       !
    case default
       stop "Error: Unsupported precission"
    end select
    !
    if (ierr         /= 0) write(*,*) "ERROR: Bessel function calculation"
    if (no_underflow /= 0) write(*,*) &
         "WARNING: Underflow occured in the calculation of the Bessel function"
    !
    ! --- the bessel function
    besselJ = besselJ_re + imu* besselJ_im
    if (order<0) besselJ = ((-1._wp)**abs(order)) * besselJ
    !
    ! --- final result 
    if (kode==2) then
       res = BesselJ(1) * exp( abs( arg_im) )
    else
       res = BesselJ(1) 
    end if
    !
    ! ------------------------------------------------------------------------
  end function BesselJn
    !


end module m_bessel
!---------------------------------------------------------------------------