program test_inter3d
use fvn_interpol
implicit none
integer(kind=ip_kind),parameter :: nx=21,ny=42,nz=18
integer(kind=ip_kind) :: i,j,k
real(kind=dp_kind) :: f,fdata(nx,ny,nz),dble,pi,q,sin,x,xdata(nx),y,ydata(ny),z,zdata(nz)
real(kind=dp_kind) :: tv
intrinsic dble,sin
f(x,y,z)=sin(x+2.*y+3.*z)
do i=1,nx
      xdata(i)=2.*(dble(i-1)/dble(nx-1))
end do
do i=1,ny
      ydata(i)=2.*(dble(i-1)/dble(ny-1))
end do
do i=1,nz
      zdata(i)=2.*(dble(i-1)/dble(nz-1))
end do
do i=1,nx
      do j=1,ny
            do k=1,nz
                  fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))
            end do
      end do
end do
call init_random_seed()
call random_number(x)
call random_number(y)
call random_number(z)
q=fvn_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)
tv=f(x,y,z)
write(*,'("x y z ",3(f8.5))') x,y,z
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