program test_besri
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='besri.dat')
write (1,*) '# n, x, bsin(n,x), b(n) in besri(x,norder+1,b)'
write (1,*) '# arg x from 0 to 20, 201 points'
write (1,*) '# order n from 0 to 20'

norder = 20
nstep = 200
xmax=20.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 besri(x,norder+1,bessvec0)
    bessvec2(i,:)=bessvec0
    do n=0,norder  !loop on rank, for dbesjn only
        bessvec1(i,n)=bsin(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 besri.dat with gnuplot :"
write(*,*) "pl 'besri.dat' u 2:($1==5 ? $3 : 1/0) w l"
write(*,*) "      will plot I5 according to bsin"
write(*,*) "pl 'besri.dat' u 2:($1==5 ? $4 : 1/0) w l"
write(*,*) "      will plot I5 according to besri"
write(*,*) "pl 'besri.dat' u 2:($1==5 ? ($4-$3)/$3 : 1/0) w l"
write(*,*) "      will plot relative shift between besri and bsin"

end program