program test_dbesrj
  use fvn_fnlib

  implicit none
  ! Variables locales -----------------------------
  integer :: i,n,nstep,norder
  real(kind=dp_kind), dimension(:), allocatable :: bessvec0
  real(kind=dp_kind), dimension(:,:), allocatable:: bessvec1,bessvec2
  real(kind=dp_kind) :: x,xstep,xmax
  
  
  open (unit=1, file='dbesrj.dat')
  write (1,*) '# n, x, bsjn(n,x), b(n) in dbesrj(x,norder+1,b)'
  write (1,*) '# arg x from 0 to 50, 501 points'
  write (1,*) '# order n from 0 to 50'
  
  norder = 50
  nstep = 500
  xmax=50.d0
  xstep=xmax/nstep
  allocate(bessvec1(0:nstep,0:norder),bessvec2(0:nstep,0:norder))
  allocate(bessvec0(0:norder))
  
  do i=0,nstep !loop on x
      x=i*xstep
      call dbesrj(x,norder+1,bessvec0)
      bessvec2(i,:)=bessvec0
      do n=0,norder  !loop on rank, for dbesjn only
          bessvec1(i,n)=bsjn(n,x)
      enddo
  enddo
  
  do n=0,norder
      do i=0,nstep
          x=i*xstep
          write (1,*) n,x,bessvec1(i,n), bessvec2(i,n)
      enddo
      write(1,*)
  enddo
  
  close(1)
  deallocate(bessvec0,bessvec1,bessvec2)
  
  END