[ Top ] [ package ]
DESCRIPTION
Set of 3D interpolators,
C and Fortran callable (C name and Fortran name are the same)
interpolate the value at point (px, py, pz) fron a 3 Dimensional real(float) array
Fortran : real, dimension(dimi,dimj,dimk) :: array
C : float array[dimk][dimj][dimi]
(Fortran storage order is expected)
the interpolation routines expect to receive the address of point
array(i[+offseti],1[+offsetj],1)
and will interpolate to get the value at point (px,py,pz)
px, py, pz are expressed in "index" space, i.e.
array(1.5, 1, 1) is halfway between array(1,1,1) and array(2,1,1)
array(1 ,1.5, 1) is halfway between array(1,1,1) and array(1,2,1)
array(1 ,1 ,2.5) is halfway between array(1,1,2) and array(1,1,3)
array(2.5,3.5,4.5) is halfway between array(2,3,4) and array(3,4,5)
constant spacing between array points is assumed for the first 2 dimensions.
along the third dimension, non constant spacing is assumed, the positions along
said third dimension are supplied via an opaque descriptor created with a call
to the Vsearch_setup family of routines (returning a C pointer to said descriptor)
4 point lagrange polynomials are used for cubic interpolations. interpolation is always
cubic along the first 2 dimensions, while along the third dimension linear interpolation
is used in the first 2 (pz <= 2) and the last 2 (pz >= nk-1) intervals, and cubic
interpolation is used in the other intervals (2.0 < pz < nk - 1)
there is a also monotonic interpolation option.
in that case, 4 results are returned
- normal interpolation result
- tri-linear interpolation result
- min and max values for the inner 2x2x2 set of values used for tri-linear interpolation
Name Purpose
Vsearch_setup : simple setup
Vsearch_setup_plus : Vsearch_setup + x and y direction offsets
Tricublin_zyx1_n : simple interpolator, source is a single variable in a 3D array
Tricublin_zyx1_n_m : source is a group of m 3D arrays
Tricublin_zyx3_n : source is a 3 variable ( dimension(3, ....) ) 3D array
Tricublin_mono_zyx_n : monotonic version of Tricublin_zyx1_n
Tricublin_mono_zyx_n_m : monotonic version of Tricublin_zyx1_n_m
Tricublin_zyx1_n_m, Tricublin_zyx3_n, Tricublin_mono_zyx_n_m save computing cycles by
reusing the computed interpolation coefficients for multiple arrays. The interpolation
coefficients are dependent upon px, py, pz.
the horizontal (first 2 dimensions) interpolation is bicubic.
(an extra linear interpolation is performed in the monotonic case)
the vertical (3rd dimension) interpolation is normally cubic, except when the target
lies between the first 2 or the last 2 levels, in which case it is linear.
the monotonic interpolators provides 4 outputs
- the result of the centered tri-cubic or bi-cubic+linear interpolation (4x4x4 or 4x4x2)
- the result of a centered tri-linear interpolation (2x2x2)
- the minimum of the 2x2x2 sub-cube used for tri-linear interpolation
- the maximum of the 2x2x2 sub-cube used for tri-linear interpolation
the interpolation routines assume that the address of the source array is the point
px = 1.0, py = 1.0, pz = 1.0 which is why f(1,1,1) is normally passed to these routines.
if this is not the case, offseti and offsetj are used to compensate
EXAMPLES
use ISO_C_BINDING include "tricublin_f90.inc" type(C_PTR) :: lv, mv ! C pointer opaque descriptor object real*8, dimension(NK) :: levels ! list of positions for the third dimension real, dimension(-2:NI+4,-1:NJ+3,NK), target :: f, f1, f2, f3 ! single valued source arrays real, dimension(3,-2:NI+4,-1:NJ+3,NK), target :: f3 ! triple valued source array type(C_PTR), dimension(3) :: f123 real, dimension(NPOINTS) :: d, l, mi, ma ! normal, linear, min, max real, dimension(3,NPOINTS) :: d3, l3, mi3, ma3 ! normal, linear, min, max real, dimension(3,NPOINTS) :: pxpypz ! (px,py,pz) triplets ...... ! levels(1:NK) = .... lv = vsearch_setup(levels, NK, NI+7, NJ+5) ! NI+7, NJ+5 to account for halos mv = vsearch_setup_plus(levels, NK, NI+7, NJ+5, 3, 2) ! passing f instead of f(1,1,1) ...... ! pxpypz(1:3,1:NPOINTS) = ...., [1.0,1.0,1.0] points to to f(1,1,1) call tricublin_zyx1_n(d ,f(1,1,1) ,pxpypz,lv,NPOINTS) ! single interpolation call tricublin_zyx1_n(d ,f ,pxpypz,mv,NPOINTS) ! with offset compensation (halos) call tricublin_zyx3_n(d3,f3(1,1,1,1),pxpypz,lv,NPOINTS) ! triple interpolation f123(1) = C_LOC(f1(1,1,1)) f123(2) = C_LOC(f2(1,1,1)) f123(3) = C_LOC(f3(1,1,1)) call tricublin_zyx1_n_m(d3,f123,pxpypz,lv,NPOINTS,3) ! 3 single interpolations call tricublin_mono_zyx_n(d,l,mi,ma,f(1,1,1),pxpypz,lv,NPOINTS) ! monotonic interpolation, 1 array call tricublin_mono_zyx_n_m(d3,l3,mi3,ma3,f123,pxpypz,lv,NPOINTS,3) ! monotonic interpolation, 3 arrays
[ Top ] [ Functions ]
Synopsis
ARGUMENTS
void Tricublin_mono_zyx_n(float *d, float *l, float *mi, float *ma, float *f, pxpypz *pxyz, ztab *lv, int n)
[ Top ] [ Functions ]
ARGUMENTS
subroutine tricublin_mono_zyx_n(d,l,mi,ma,f,pxyz,lv,n) bind(C,name='Tricublin_mono_zyx_n')
import :: C_PTR
real, dimension(*), intent(OUT) :: d, l, mi, ma
real, dimension(*), intent(IN) :: f
real, dimension(*), intent(IN) :: pxyz
type(C_PTR), intent(IN), value :: lv
integer, intent(IN), value :: n
[ Top ] [ Functions ]
Synopsis
ARGUMENTS
void Tricublin_mono_zyx_n_m(float *d, float *l, float *mi, float *ma, float **fs, pxpypz *pxyz, ztab *lv, int n, int m)
[ Top ] [ Functions ]
ARGUMENTS
subroutine tricublin_mono_zyx_n_m(d,l,mi,ma,f,pxyz,lv,n,m) bind(C,name='Tricublin_mono_zyx_n_m')
import :: C_PTR
real, dimension(*), intent(OUT) :: d, l, mi, ma
type(C_PTR), dimension(*), intent(IN) :: f
real, dimension(*), intent(IN) :: pxyz
type(C_PTR), intent(IN), value :: lv
integer, intent(IN), value :: n, m
[ Top ] [ Functions ]
Synopsis
ARGUMENTS
void Tricublin_zyx1_n(float *d, float *f1, pxpypz *pxyz, ztab *lv, int n)
[ Top ] [ Functions ]
Synopsis
f real input array of dimension(xmin:xmax,ymin,ymax,nk)
f(1,1,1) MUST be passed to this function
xmax MUST be xmin+ni-1, ymax MUST be ymin+nj-1 (see setup routines)
pxyz real array, dimension(3,n), contains the target interpolation positions
positions are in global index space
pxyz(1,i) is the position along x of point i (i = 1..n)
pxyz(2,i) is the position along y of point i
pxyz(3,i) is the position along z of point i
d real array, dimension(n), d(i) is the interpolated value at point i
lv opaque object obtained from the setup routines
n number of points to interpolate
ARGUMENTS
subroutine tricublin_zyx1_n(d,f,pxyz,lv,n) bind(C,name='Tricublin_zyx1_n')
import :: C_PTR
real, dimension(n), intent(OUT) :: d
real, dimension(*), intent(IN) :: f
real, dimension(3,n), intent(IN) :: pxyz
type(C_PTR), intent(IN), value :: lv
integer, intent(IN), value :: n
[ Top ] [ Functions ]
Synopsis
ARGUMENTS
void Tricublin_zyx1_n_m(float *d, float **fs, pxpypz *pxyz, ztab *lv, int n, int m)
[ Top ] [ Functions ]
ARGUMENTS
subroutine tricublin_zyx1_n_m(d,f,pxyz,lv,n,m) bind(C,name='Tricublin_zyx1_n_m')
import :: C_PTR
real, dimension(*), intent(OUT) :: d
type(C_PTR), dimension(*), intent(IN) :: f
real, dimension(*), intent(IN) :: pxyz
type(C_PTR), intent(IN), value :: lv
integer, intent(IN), value :: n, m
[ Top ] [ Functions ]
Synopsis
ARGUMENTS
void Tricublin_zyx3_n(float *d, float *f3, pxpypz *pxyz, ztab *lv, int n)
[ Top ] [ Functions ]
ARGUMENTS
subroutine tricublin_zyx3_n(d,f,pxyz,lv,n) bind(C,name='Tricublin_zyx3_n')
import :: C_PTR
real, dimension(*), intent(OUT) :: d
real, dimension(*), intent(IN) :: f
real, dimension(*), intent(IN) :: pxyz
type(C_PTR), intent(IN), value :: lv
integer, intent(IN), value :: n
[ Top ] [ Functions ]
Synopsis
ARGUMENTS
ztab *Vsearch_setup(double *targets, int nk, int ni, int nj)
[ Top ] [ Functions ]
Synopsis
the function returns an opaque object of type C_PTR, to be passed to the interpolators
levels : real*8 array of levels, dimension(nk)
nk : number of levels
ni, nj : used to compute indexing into source array(s) for interpolators
ni is used to compute the distance between 2 array elements for which index j differs by 1
ni*nj is used to compute the distance between 2 array elements for which index k differs by 1
when vsearch_setup is used, positions will be assumed to be in local index space
(equivalent to offseti = 0 and offsetj = 0 in vsearch_setup_plus)
see interpolators for more information
ARGUMENTS
function vsearch_setup(levels, nk, ni, nj) result (ztab) bind(C,name='Vsearch_setup')
import :: C_PTR, C_DOUBLE, C_INT
real(C_DOUBLE), dimension(nk), intent(IN) :: levels
integer(C_INT), intent(IN), value :: nk, ni, nj
type(C_PTR) :: ztab
[ Top ] [ Functions ]
Synopsis
ARGUMENTS
ztab *Vsearch_setup_plus(double *targets, int nk, int ni, int nj, int offseti, int offsetj)
[ Top ] [ Functions ]
Synopsis
the function returns an opaque object of type C_PTR, to be passed to the interpolators
levels : real*8 array of levels, dimension(nk)
nk : number of levels
ni, nj : used to compute indexing into source array(s) for interpolators
offseti : used to convert global position along i into local array position
offsetj : used to convert global position along j into local array position
ni is used to compute the distance between 2 array elements for which index j differs by 1
ni*nj is used to compute the distance between 2 array elements for which index k differs by 1
effective px = px + offseti
effective py = py + offsetj
see interpolators for more information
ARGUMENTS
function vsearch_setup_plus(levels, nk, ni, nj, offseti, offsetj) result (ztab) bind(C,name='Vsearch_setup_plus')
import :: C_PTR, C_DOUBLE, C_INT
real(C_DOUBLE), dimension(nk), intent(IN) :: levels
integer(C_INT), intent(IN), value :: nk, ni, nj
integer(C_INT), intent(IN), value :: offseti, offsetj
type(C_PTR) :: ztab