program test_bskn
  use fvn_fnlib
  implicit none
  integer :: n,npoints,i

  real(kind=sp_kind) :: xmin,xmax,xstep,x,y
  real(kind=sp_kind),dimension(200,0:5) :: bes

  open(2,file='bskn.dat')
  
  xmin=0.
  xmax=10.
  npoints=200
  xstep=(xmax-xmin)/dble(npoints)
  do i=1,npoints
        x=xmin+i*xstep
        do n=0,5
          bes(i,n)=bskn(n,x)
        end do
        write(2,'(7e22.14)') x,bes(i,0:5)
  end do
  close(2)
  
  ! with gnuplot -> set yrange [-1:1]
  
  end program