program inter2d

  use fvn_interpol

  implicit none

  integer(kind=ip_kind),parameter :: nx=21,ny=42
  integer(kind=ip_kind) :: i,j
  real(kind=dp_kind) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)
  real(kind=dp_kind) :: 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