program test_besrj
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='besrj.dat')
write (1,*) '# n, x, bsjn(n,x), b(n) in besrj(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 besrj(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)

write(*,*) "Exemple use of generated file besrj.dat with gnuplot :"
write(*,*) "pl 'besrj.dat' u 2:($1==5 ? $3 : 1/0) w l"
write(*,*) "      will plot J5 according to bsjn"
write(*,*) "pl 'besrj.dat' u 2:($1==5 ? $4 : 1/0) w l"
write(*,*) "      will plot J5 according to besrj"
write(*,*) "pl 'besrj.dat' u 2:($1==5 ? ($4-$3)/$3 : 1/0) w l"
write(*,*) "      will plot relative shift between besrj and bsjn"



end program