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