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