program inter2d
use fvn_interpol
implicit none
integer(kind=4),parameter :: nx=21,ny=42
integer(kind=4) :: i,j
real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)
real(kind=8) :: tv
intrinsic dble,sin
f(x,y)=sin(x+2.*y)
do i=1,nx
xdata(i)=dble(i-1)/dble(nx-1)
end do
do i=1,ny
ydata(i)=dble(i-1)/dble(ny-1)
end do
do i=1,nx
do j=1,ny
fdata(i,j)=f(xdata(i),ydata(j))
end do
end do
call init_random_seed()
call random_number(x)
call random_number(y)
q=fvn_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata)
tv=f(x,y)
write(*,'("x y z ",2(f8.5))') x,y
write(*,'("Calculated (real) value :",f8.5)') tv
write(*,'("fvn interpolation : ",f8.5)') q
write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)
end program