program test_bsjn
  use fvn_fnlib
  implicit none
  integer :: n,npoints,i
  real(4) :: xmin,xmax,xstep,x,y
  real(4),dimension(200,6) :: bes
  !    interface
  !    real(4) function besjn(n,x,factor,big)
  !        integer :: n
  !        real(4) :: x
  !        integer,optional :: factor
  !        real(4),optional :: big
  !    end function besjn
  !end interface
  !bsyn
  !y=besjn(2,0.5)
  !stop
  open(2,file='bsjn.dat')
  
  xmin=-20.
  xmax=20.
  npoints=200
  xstep=(xmax-xmin)/dble(npoints)
  do i=1,npoints
        x=xmin+i*xstep
        do n=0,5
          bes(i,n)=bsjn(n,x)
        end do
        write(2,'(7e22.14)') x,bes(i,0:5)
  end do
  close(2)
  
  end program