Develop

README.md

@@ -70,6 +70,15 @@ type(bspline_3d) :: s
 s = bspline_3d(x,y,z,fcn,kx,ky,kz,iflag,extrap)
 ```
 
+## Spline order
+
+The various `k` inputs (i.e., `kx`, `ky`, etc.) specify the spline order for each dimension. The order is the polynomial degree + 1. For example:
+
+* `k=2` : Linear
+* `k=3` : Quadratic
+* `k=4` : Cubic
+* etc.
+
 ## Extrapolation
 
 The library optionally supports extrapolation for points outside the range of the coefficients. This is disabled by default (in which case an error code is returned for points outside the bounds). To enable extrapolation, use the optional `extrap` input to the various `db*val` subroutines or the `initialize` methods from the object-oriented interface.
@@ -180,3 +189,10 @@ The latest API documentation can be found [here](https://jacobwilliams.github.io
 ## License
 
 The bspline-fortran source code and related files and documentation are distributed under a permissive free software [license](https://github.com/jacobwilliams/bspline-fortran/blob/master/LICENSE) (BSD-style).
+
+## See also
+
+ * [SPLPAK](https://github.com/jacobwilliams/splpak) Multidimensional least-squares cubic spline fitting
+ * [FINTERP](https://github.com/jacobwilliams/finterp) Multidimensional Linear Interpolation with Modern Fortran
+ * [PCHIP](https://github.com/jacobwilliams/PCHIP) Piecewise Cubic Hermite Interpolation.
+ * [Regridpack](https://github.com/jacobwilliams/regridpack) Linear or cubic interpolation for 1D-4D grids.

fpm.toml

@@ -14,6 +14,7 @@ auto-tests = true
 
 [dev-dependencies]
 pyplot-fortran = { git = "https://github.com/jacobwilliams/pyplot-fortran", tag = "3.3.0" }
+finterp = { git = "https://github.com/jacobwilliams/finterp", tag = "1.3.0" }
 
 [library]
 source-dir = "src"

src/bspline_sub_module.f90

@@ -58,6 +58,9 @@ end function b1fqad_func
     end interface
     public :: b1fqad_func
 
+    integer(ip),parameter,public :: bspline_order_linear    = 2_ip !! spline order `k` parameter
+                                                                   !! (for input to the `db*ink` routines)
+                                                                   !! [order = polynomial degree + 1]
     integer(ip),parameter,public :: bspline_order_quadratic = 3_ip !! spline order `k` parameter
                                                                    !! (for input to the `db*ink` routines)
                                                                    !! [order = polynomial degree + 1]

test/bspline_extrap_test.f90

@@ -22,7 +22,7 @@ program bspline_extrap_test
     real(wp) :: x(nx)
     real(wp) :: xval(nxv),fval(nxv)
     real(wp) :: tx(nx+kx)
-    real(wp) :: fcn_1d(nx)
+    real(wp) :: fcn_1d(nx), bcoef(nx)
     real(wp) :: val,tru,err,errmax
     integer(ip) :: i,j,idx,iflag,inbvx,iloy
     logical :: extrap
@@ -47,7 +47,7 @@ program bspline_extrap_test
     iloy  = 1
 
     ! initialize
-    call db1ink(x,nx,fcn_1d,kx,iknot,tx,fcn_1d,iflag)
+    call db1ink(x,nx,fcn_1d,kx,iknot,tx,bcoef,iflag)
 
     if (iflag/=0) then
         write(*,*) 'Error initializing 1D spline: '//get_status_message(iflag)
@@ -67,7 +67,7 @@ program bspline_extrap_test
 
         errmax = 0.0_wp
         do i=1,nxv
-            call db1val(xval(i),idx,tx,nx,kx,fcn_1d,val,iflag,inbvx,w1_1d,extrap=extrap)
+            call db1val(xval(i),idx,tx,nx,kx,bcoef,val,iflag,inbvx,w1_1d,extrap=extrap)
             fval(i) = val  ! save it for plot
             tru    = f1(xval(i))
             err    = abs(tru-val)

test/bspline_linear_test.f90

@@ -0,0 +1,106 @@
+!*****************************************************************************************
+!>
+!  Linear interpolation test
+
+    program bspline_linear_test
+
+    use bspline_module
+    use bspline_kinds_module, only: wp, ip
+    use pyplot_module
+    use linear_interpolation_module
+
+    implicit none
+
+    integer(ip),parameter :: nx    = 10    !! number of points in x
+    integer(ip),parameter :: nxv   = 111   !! number of points to evaluate interpolant
+    integer(ip),parameter :: iknot = 0    !! automatically select the knots
+    integer,parameter :: idx = 0 !! compute value for the spline interpolation
+
+    real(wp) :: x(nx)
+    real(wp) :: xval(nxv),fval(nxv),fval_linear(nxv),fval3(nxv),fval4(nxv)
+    real(wp) :: fcn_1d(nx)
+    real(wp) :: val,err,errmax
+    integer(ip) :: i,iflag,inbvx,iloy
+    integer :: istat  !! pyplot-fortran status flag
+    type(pyplot) :: plt
+    type(bspline_1d) :: b2, b3, b4
+    type(linear_interp_1d) :: s1
+
+    do i=1,nx
+        x(i) = real(i-1,wp) * 10.0_wp
+        fcn_1d(i) = f1(x(i))
+    end do
+    do i=1,nxv
+        xval(i) = real(i-1,wp) - 10
+    end do
+
+    !have to set these before the first evaluate call:
+    inbvx = 1
+    iloy  = 1
+
+    ! initialize
+    call b2%initialize(x,fcn_1d,bspline_order_linear,iflag,extrap=.true.) ! linear
+    if (iflag/=0) error stop 'Error initializing 1D linear spline: '//get_status_message(iflag)
+    call b3%initialize(x,fcn_1d,bspline_order_quadratic,iflag,extrap=.true.) ! quadratic
+    if (iflag/=0) error stop 'Error initializing 1D quadratic spline: '//get_status_message(iflag)
+    call b4%initialize(x,fcn_1d,bspline_order_cubic,iflag,extrap=.true.) ! cubic
+    if (iflag/=0) error stop 'Error initializing 1D cubic spline: '//get_status_message(iflag)
+    call s1%initialize(x,fcn_1d,iflag)
+    if (iflag/=0) error stop 'Error initializing 1D linear interpolator'
+
+    !initialize the plot:
+    call plt%initialize(grid=.true.,xlabel='x (deg)',ylabel='f(x)',&
+                        title='Linear Test',legend=.true.,figsize=[10,5])
+    call plt%add_plot(x,fcn_1d,&
+                        label='Function $f(x) = \sin(x)$',&
+                        linestyle='ko',markersize=5,linewidth=2,istat=istat)
+
+    errmax = 0.0_wp
+    do i=1,nxv
+        call b2%evaluate(xval(i),idx,val,iflag)
+        fval(i) = val  ! save it for plot
+        if (iflag/=0) error stop 'error evaluating linear spline: '//get_status_message(iflag)
+
+        ! also compute linear interpolation:
+        call s1%evaluate(xval(i),val)
+        fval_linear(i) = val ! linear vector for plot
+
+        err    = abs(fval(i) - fval_linear(i))
+        errmax = max(err,errmax)
+
+        ! also others:
+        call b3%evaluate(xval(i),idx,fval3(i),iflag)
+        if (iflag/=0) error stop 'error evaluating quadratic spline: '//get_status_message(iflag)
+        call b4%evaluate(xval(i),idx,fval4(i),iflag)
+        if (iflag/=0) error stop 'error evaluating cubic spline: '//get_status_message(iflag)
+
+    end do
+
+    write(*,*) ''
+    write(*,*) 'Max difference (spline - linear):', errmax
+    write(*,*) ''
+
+    call plt%add_plot(xval,fval_linear,&
+            label='Linear Interpolated',&
+            linestyle='r.',linewidth=1,istat=istat)
+    call plt%add_plot(xval,fval,&
+            label='Linear ($k=2$) Spline Interpolated',&
+            linestyle='r-',linewidth=1,istat=istat)
+    call plt%add_plot(xval,fval3,&
+            label='Quadratic ($k=3$) Spline Interpolated',&
+            linestyle='k.-',linewidth=1,istat=istat)
+    call plt%add_plot(xval,fval4,&
+            label='Cubic ($k=4$) Spline Interpolated',&
+            linestyle='c.-',linewidth=1,istat=istat)
+    call plt%savefig('bspline_linear_test.png',istat=istat)
+
+    contains
+
+        real(wp) function f1(x) !! 1d test function
+        implicit none
+        real(wp),intent(in) :: x
+        real(wp),parameter :: a = acos(-1.0_wp)/18.0_wp
+        f1 = sin(a*x)
+        end function f1
+
+    end program bspline_linear_test

test/bspline_test.f90

@@ -27,12 +27,12 @@ program bspline_test
 
     real(wp) :: x(nx),y(ny),z(nz),q(nq),r(nr),s(ns)
     real(wp) :: tx(nx+kx),ty(ny+ky),tz(nz+kz),tq(nq+kq),tr(nr+kr),ts(ns+ks)
-    real(wp) :: fcn_1d(nx)
-    real(wp) :: fcn_2d(nx,ny)
-    real(wp) :: fcn_3d(nx,ny,nz)
-    real(wp) :: fcn_4d(nx,ny,nz,nq)
-    real(wp) :: fcn_5d(nx,ny,nz,nq,nr)
-    real(wp) :: fcn_6d(nx,ny,nz,nq,nr,ns)
+    real(wp) :: fcn_1d(nx)               , bcoef_1d(nx)
+    real(wp) :: fcn_2d(nx,ny)            , bcoef_2d(nx,ny)
+    real(wp) :: fcn_3d(nx,ny,nz)         , bcoef_3d(nx,ny,nz)
+    real(wp) :: fcn_4d(nx,ny,nz,nq)      , bcoef_4d(nx,ny,nz,nq)
+    real(wp) :: fcn_5d(nx,ny,nz,nq,nr)   , bcoef_5d(nx,ny,nz,nq,nr)
+    real(wp) :: fcn_6d(nx,ny,nz,nq,nr,ns), bcoef_6d(nx,ny,nz,nq,nr,ns)
 
     real(wp),dimension(3*kx)                     :: w1_1d
     real(wp),dimension(ky)                       :: w1_2d
@@ -229,12 +229,12 @@ program bspline_test
     ilos  = 1
 
     ! initialize
-    call db1ink(x,nx,fcn_1d,kx,iknot,tx,fcn_1d,iflag(1))
-    call db2ink(x,nx,y,ny,fcn_2d,kx,ky,iknot,tx,ty,fcn_2d,iflag(2))
-    call db3ink(x,nx,y,ny,z,nz,fcn_3d,kx,ky,kz,iknot,tx,ty,tz,fcn_3d,iflag(3))
-    call db4ink(x,nx,y,ny,z,nz,q,nq,fcn_4d,kx,ky,kz,kq,iknot,tx,ty,tz,tq,fcn_4d,iflag(4))
-    call db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,fcn_5d,kx,ky,kz,kq,kr,iknot,tx,ty,tz,tq,tr,fcn_5d,iflag(5))
-    call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,fcn_6d,kx,ky,kz,kq,kr,ks,iknot,tx,ty,tz,tq,tr,ts,fcn_6d,iflag(6))
+    call db1ink(x,nx,fcn_1d,kx,iknot,tx,bcoef_1d,iflag(1))
+    call db2ink(x,nx,y,ny,fcn_2d,kx,ky,iknot,tx,ty,bcoef_2d,iflag(2))
+    call db3ink(x,nx,y,ny,z,nz,fcn_3d,kx,ky,kz,iknot,tx,ty,tz,bcoef_3d,iflag(3))
+    call db4ink(x,nx,y,ny,z,nz,q,nq,fcn_4d,kx,ky,kz,kq,iknot,tx,ty,tz,tq,bcoef_4d,iflag(4))
+    call db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,fcn_5d,kx,ky,kz,kq,kr,iknot,tx,ty,tz,tq,tr,bcoef_5d,iflag(5))
+    call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,fcn_6d,kx,ky,kz,kq,kr,ks,iknot,tx,ty,tz,tq,tr,ts,bcoef_6d,iflag(6))
 
     if (any(iflag/=0)) then
         do i=1,6
@@ -249,46 +249,46 @@ program bspline_test
      errmax = 0.0_wp
      do i=1,nx
                         call db1val(x(i),idx,&
-                                            tx,nx,kx,fcn_1d,val(1),iflag(1),inbvx,&
+                                            tx,nx,kx,bcoef_1d,val(1),iflag(1),inbvx,&
                                             w1_1d)
                         tru(1)    = f1(x(i))
                         err(1)    = abs(tru(1)-val(1))
                         errmax(1) = max(err(1),errmax(1))
         do j=1,ny
                         call db2val(x(i),y(j),idx,idy,&
-                                            tx,ty,nx,ny,kx,ky,fcn_2d,val(2),iflag(2),&
+                                            tx,ty,nx,ny,kx,ky,bcoef_2d,val(2),iflag(2),&
                                             inbvx,inbvy,iloy,&
                                             w1_2d,w2_2d)
                         tru(2)    = f2(x(i),y(j))
                         err(2)    = abs(tru(2)-val(2))
                         errmax(2) = max(err(2),errmax(2))
            do k=1,nz
                         call db3val(x(i),y(j),z(k),idx,idy,idz,&
-                                            tx,ty,tz,nx,ny,nz,kx,ky,kz,fcn_3d,val(3),iflag(3),&
+                                            tx,ty,tz,nx,ny,nz,kx,ky,kz,bcoef_3d,val(3),iflag(3),&
                                             inbvx,inbvy,inbvz,iloy,iloz,&
                                             w1_3d,w2_3d,w3_3d)
                         tru(3)    = f3(x(i),y(j),z(k))
                         err(3)    = abs(tru(3)-val(3))
                         errmax(3) = max(err(3),errmax(3))
               do l=1,nq
                         call db4val(x(i),y(j),z(k),q(l),idx,idy,idz,idq,&
-                                            tx,ty,tz,tq,nx,ny,nz,nq,kx,ky,kz,kq,fcn_4d,val(4),iflag(4),&
+                                            tx,ty,tz,tq,nx,ny,nz,nq,kx,ky,kz,kq,bcoef_4d,val(4),iflag(4),&
                                             inbvx,inbvy,inbvz,inbvq,iloy,iloz,iloq,&
                                             w1_4d,w2_4d,w3_4d,w4_4d)
                         tru(4)    = f4(x(i),y(j),z(k),q(l))
                         err(4)    = abs(tru(4)-val(4))
                         errmax(4) = max(err(4),errmax(4))
                 do m=1,nr
                         call db5val(x(i),y(j),z(k),q(l),r(m),idx,idy,idz,idq,idr,&
-                                            tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,fcn_5d,val(5),iflag(5),&
+                                            tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,bcoef_5d,val(5),iflag(5),&
                                             inbvx,inbvy,inbvz,inbvq,inbvr,iloy,iloz,iloq,ilor,&
                                             w1_5d,w2_5d,w3_5d,w4_5d,w5_5d)
                         tru(5)    = f5(x(i),y(j),z(k),q(l),r(m))
                         err(5)    = abs(tru(5)-val(5))
                         errmax(5) = max(err(5),errmax(5))
                     do n=1,ns
                         call db6val(x(i),y(j),z(k),q(l),r(m),s(n),idx,idy,idz,idq,idr,ids,&
-                                            tx,ty,tz,tq,tr,ts,nx,ny,nz,nq,nr,ns,kx,ky,kz,kq,kr,ks,fcn_6d,val(6),iflag(6),&
+                                            tx,ty,tz,tq,tr,ts,nx,ny,nz,nq,nr,ns,kx,ky,kz,kq,kr,ks,bcoef_6d,val(6),iflag(6),&
                                             inbvx,inbvy,inbvz,inbvq,inbvr,inbvs,iloy,iloz,iloq,ilor,ilos,&
                                             w1_6d,w2_6d,w3_6d,w4_6d,w5_6d,w6_6d)
                         tru(6)    = f6(x(i),y(j),z(k),q(l),r(m),s(n))

test/bspline_test_2.f90

@@ -30,12 +30,12 @@ program bspline_test_2
     real(wp) :: x(nx),y(ny),z(nz),q(nq),r(nr),s(ns)
     real(wp) :: xval(nx+2),yval(ny+2),zval(nz+2),qval(nq+2),rval(nr+2),sval(ns+2)
     real(wp) :: tx(nx+kx),ty(ny+ky),tz(nz+kz),tq(nq+kq),tr(nr+kr),ts(ns+ks)
-    real(wp) :: fcn_1d(nx)
-    real(wp) :: fcn_2d(nx,ny)
-    real(wp) :: fcn_3d(nx,ny,nz)
-    real(wp) :: fcn_4d(nx,ny,nz,nq)
-    real(wp) :: fcn_5d(nx,ny,nz,nq,nr)
-    real(wp) :: fcn_6d(nx,ny,nz,nq,nr,ns)
+    real(wp) :: fcn_1d(nx)               , bcoef_1d(nx)
+    real(wp) :: fcn_2d(nx,ny)            , bcoef_2d(nx,ny)
+    real(wp) :: fcn_3d(nx,ny,nz)         , bcoef_3d(nx,ny,nz)
+    real(wp) :: fcn_4d(nx,ny,nz,nq)      , bcoef_4d(nx,ny,nz,nq)
+    real(wp) :: fcn_5d(nx,ny,nz,nq,nr)   , bcoef_5d(nx,ny,nz,nq,nr)
+    real(wp) :: fcn_6d(nx,ny,nz,nq,nr,ns), bcoef_6d(nx,ny,nz,nq,nr,ns)
 
     real(wp),dimension(3*kx)                     :: w1_1d
     real(wp),dimension(ky)                       :: w1_2d
@@ -146,12 +146,12 @@ program bspline_test_2
 
     ! initialize
     iflag = 0
-    call db1ink(x,nx,fcn_1d,kx,iknot,tx,fcn_1d,iflag(1))
-    call db2ink(x,nx,y,ny,fcn_2d,kx,ky,iknot,tx,ty,fcn_2d,iflag(2))
-    call db3ink(x,nx,y,ny,z,nz,fcn_3d,kx,ky,kz,iknot,tx,ty,tz,fcn_3d,iflag(3))
-    call db4ink(x,nx,y,ny,z,nz,q,nq,fcn_4d,kx,ky,kz,kq,iknot,tx,ty,tz,tq,fcn_4d,iflag(4))
-    call db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,fcn_5d,kx,ky,kz,kq,kr,iknot,tx,ty,tz,tq,tr,fcn_5d,iflag(5))
-    call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,fcn_6d,kx,ky,kz,kq,kr,ks,iknot,tx,ty,tz,tq,tr,ts,fcn_6d,iflag(6))
+    call db1ink(x,nx,fcn_1d,kx,iknot,tx,bcoef_1d,iflag(1))
+    call db2ink(x,nx,y,ny,fcn_2d,kx,ky,iknot,tx,ty,bcoef_2d,iflag(2))
+    call db3ink(x,nx,y,ny,z,nz,fcn_3d,kx,ky,kz,iknot,tx,ty,tz,bcoef_3d,iflag(3))
+    call db4ink(x,nx,y,ny,z,nz,q,nq,fcn_4d,kx,ky,kz,kq,iknot,tx,ty,tz,tq,bcoef_4d,iflag(4))
+    call db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,fcn_5d,kx,ky,kz,kq,kr,iknot,tx,ty,tz,tq,tr,bcoef_5d,iflag(5))
+    call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,fcn_6d,kx,ky,kz,kq,kr,ks,iknot,tx,ty,tz,tq,tr,ts,bcoef_6d,iflag(6))
 
     if (any(iflag/=0)) then
         do i=1,6
@@ -168,30 +168,30 @@ program bspline_test_2
      err = -99999.9_wp
      do i=1,nx+2
                         call db1val(xval(i),idx,&
-                                            tx,nx,kx,fcn_1d,val(1),iflag(1),inbvx,w1_1d,extrap=.true.)
+                                            tx,nx,kx,bcoef_1d,val(1),iflag(1),inbvx,w1_1d,extrap=.true.)
                         tru(1)    = f1(xval(i))
                         err(1)    = abs(tru(1)-val(1))
                         errmax(1) = max(err(1),errmax(1))
                         if (iflag(1)/=0) write(*,*) 'Error evaluating 1D spline: '//get_status_message(iflag(1))
         do j=1,ny+2
                         call db2val(xval(i),yval(j),idx,idy,&
-                                            tx,ty,nx,ny,kx,ky,fcn_2d,val(2),iflag(2),&
+                                            tx,ty,nx,ny,kx,ky,bcoef_2d,val(2),iflag(2),&
                                             inbvx,inbvy,iloy,w1_2d,w2_2d,extrap=.true.)
                         tru(2)    = f2(xval(i),yval(j))
                         err(2)    = abs(tru(2)-val(2))
                         errmax(2) = max(err(2),errmax(2))
                         if (iflag(2)/=0) write(*,*) 'Error evaluating 2D spline: '//get_status_message(iflag(2))
            do k=1,nz+2
                         call db3val(xval(i),yval(j),zval(k),idx,idy,idz,&
-                                            tx,ty,tz,nx,ny,nz,kx,ky,kz,fcn_3d,val(3),iflag(3),&
+                                            tx,ty,tz,nx,ny,nz,kx,ky,kz,bcoef_3d,val(3),iflag(3),&
                                             inbvx,inbvy,inbvz,iloy,iloz,w1_3d,w2_3d,w3_3d,extrap=.true.)
                         tru(3)    = f3(xval(i),yval(j),zval(k))
                         err(3)    = abs(tru(3)-val(3))
                         errmax(3) = max(err(3),errmax(3))
                         if (iflag(3)/=0) write(*,*) 'Error evaluating 3D spline: '//get_status_message(iflag(3))
               do l=1,nq+2
                         call db4val(xval(i),yval(j),zval(k),qval(l),idx,idy,idz,idq,&
-                                            tx,ty,tz,tq,nx,ny,nz,nq,kx,ky,kz,kq,fcn_4d,val(4),iflag(4),&
+                                            tx,ty,tz,tq,nx,ny,nz,nq,kx,ky,kz,kq,bcoef_4d,val(4),iflag(4),&
                                             inbvx,inbvy,inbvz,inbvq,iloy,iloz,iloq,&
                                             w1_4d,w2_4d,w3_4d,w4_4d,extrap=.true.)
                         tru(4)    = f4(xval(i),yval(j),zval(k),qval(l))
@@ -200,7 +200,7 @@ program bspline_test_2
                         if (iflag(4)/=0) write(*,*) 'Error evaluating 4D spline: '//get_status_message(iflag(4))
                 do m=1,nr+2
                         call db5val(xval(i),yval(j),zval(k),qval(l),rval(m),idx,idy,idz,idq,idr,&
-                                            tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,fcn_5d,val(5),iflag(5),&
+                                            tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,bcoef_5d,val(5),iflag(5),&
                                             inbvx,inbvy,inbvz,inbvq,inbvr,iloy,iloz,iloq,ilor,&
                                             w1_5d,w2_5d,w3_5d,w4_5d,w5_5d,extrap=.true.)
                         tru(5)    = f5(xval(i),yval(j),zval(k),qval(l),rval(m))
@@ -209,7 +209,7 @@ program bspline_test_2
                         if (iflag(5)/=0) write(*,*) 'Error evaluating 5D spline: '//get_status_message(iflag(5))
                     do n=1,ns+2
                         call db6val(xval(i),yval(j),zval(k),qval(l),rval(m),sval(n),idx,idy,idz,idq,idr,ids,&
-                                            tx,ty,tz,tq,tr,ts,nx,ny,nz,nq,nr,ns,kx,ky,kz,kq,kr,ks,fcn_6d,val(6),iflag(6),&
+                                            tx,ty,tz,tq,tr,ts,nx,ny,nz,nq,nr,ns,kx,ky,kz,kq,kr,ks,bcoef_6d,val(6),iflag(6),&
                                             inbvx,inbvy,inbvz,inbvq,inbvr,inbvs,iloy,iloz,iloq,ilor,ilos,&
                                             w1_6d,w2_6d,w3_6d,w4_6d,w5_6d,w6_6d,extrap=.true.)
                         tru(6)    = f6(xval(i),yval(j),zval(k),qval(l),rval(m),sval(n))