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