program test_specfunc
use fvn_fnlib
implicit none
integer,parameter :: npoints=200
integer :: i
real(8), dimension(npoints) :: j0
real(8) :: xmin,xmax,xstep,x,y

! bsj0
xmin=-20.
xmax=20.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='bsj0.dat')
do i=1,npoints
x=xmin+i*xstep
y=bsj0(x)
write(2,'(2e22.14)') x,y
end do
close(2)

! bsj1
xmin=-20.
xmax=20.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='bsj1.dat')
do i=1,npoints
x=xmin+i*xstep
y=bsj1(x)
write(2,'(2e22.14)') x,y
end do
close(2)

!bsi0
xmin=-4.
xmax=4.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='bsi0.dat')
do i=1,npoints
x=xmin+i*xstep
y=bsi0(x)
write(2,'(2e22.14)') x,y
end do
close(2)

!bsi1
xmin=-4.
xmax=4.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='bsi1.dat')
do i=1,npoints
x=xmin+i*xstep
y=bsi1(x)
write(2,'(2e22.14)') x,y
end do
close(2)

!bsy0
xmin=0.
xmax=20.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='bsy0.dat')
do i=1,npoints
x=xmin+i*xstep
y=bsy0(x)
write(2,'(2e22.14)') x,y
end do
close(2)

!bsy1
xmin=0.
xmax=20.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='bsy1.dat')
do i=1,npoints
x=xmin+i*xstep
y=bsy1(x)
write(2,'(2e22.14)') x,y
end do
close(2)

!bsk0
xmin=0.
xmax=4.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='bsk0.dat')
do i=1,npoints
x=xmin+i*xstep
y=bsk0(x)
write(2,'(2e22.14)') x,y
end do
close(2)

!bsk1
xmin=0.
xmax=4.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='bsk1.dat')
do i=1,npoints
x=xmin+i*xstep
y=bsk1(x)
write(2,'(2e22.14)') x,y
end do
close(2)

!erf
xmin=-4.
xmax=4.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='erf.dat')
do i=1,npoints
x=xmin+i*xstep
y=erf(x)
write(2,'(2e22.14)') x,y
end do
close(2)

! gamma
xmin=-3.
xmax=7.
xstep=(xmax-xmin)/2000.
open(2,file='gamma.dat')
do i=1,2000
x=xmin+i*xstep
if ((abs(x-nint(x)) >= 1d-6) .or. x>0. )then
    y=gamma(x)
    write(2,'(2e22.14)') x,y
end if
end do
close(2)

! 1/gamma
xmin=-3.
xmax=7.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='gamr.dat')
do i=1,npoints
x=xmin+i*xstep
y=gamr(x)
write(2,'(2e22.14)') x,y
end do
close(2)

! ei
xmin=0.
xmax=1.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='ei.dat')
do i=1,npoints
x=xmin+i*xstep
y=ei(x)
write(2,'(2e22.14)') x,y
end do
close(2)

! e1
xmin=0.
xmax=1.
xstep=(xmax-xmin)/dble(npoints)
open(2,file='e1.dat')
do i=1,npoints
x=xmin+i*xstep
y=e1(x)
write(2,'(2e22.14)') x,y
end do
close(2)



end program test_specfunc