1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
#include "cctk.h"
SUBROUTINE intp3d_many ( x, y, z, f, kt, ft, nx, ny, nz, nvars, xt, yt, zt)
!
implicit none
!
!---------------------------------------------------------------------
!
! purpose: interpolation of a function of three variables in an
! equidistant(!!!) table.
!
! method: 8-point Lagrange linear interpolation formula
!
! x input vector of first variable
! y input vector of second variable
! z input vector of third variable
!
! f output vector of interpolated function values
!
! kt vector length of input and output vectors
!
! ft 3d array of tabulated function values
! nx x-dimension of table
! ny y-dimension of table
! nz z-dimension of table
! xt vector of x-coordinates of table
! yt vector of y-coordinates of table
! zt vector of z-coordinates of table
!
!---------------------------------------------------------------------
integer kt,nx,ny,nz,iv,nvars
real*8 :: ft(nx,ny,nz,nvars)
real*8 x(kt),y(kt),z(kt),f(kt,nvars)
real*8 xt(nx),yt(ny),zt(nz)
real*8 d1,d2,d3
!
!
integer,parameter :: ktx = 1
real*8 fh(ktx,8,nvars), delx(ktx), dely(ktx), delz(ktx), &
a1(ktx,nvars), a2(ktx,nvars), a3(ktx,nvars), a4(ktx,nvars), &
a5(ktx,nvars), a6(ktx,nvars), a7(ktx,nvars), a8(ktx,nvars)
real*8 dx,dy,dz,dxi,dyi,dzi,dxyi,dxzi,dyzi,dxyzi
integer n,ix,iy,iz
IF (kt .GT. ktx) call CCTK_WARN(0, '***KTX**')
!
!
!------ determine spacing parameters of (equidistant!!!) table
!
dx = (xt(nx) - xt(1)) / FLOAT(nx-1)
dy = (yt(ny) - yt(1)) / FLOAT(ny-1)
dz = (zt(nz) - zt(1)) / FLOAT(nz-1)
!
dxi = 1. / dx
dyi = 1. / dy
dzi = 1. / dz
!
dxyi = dxi * dyi
dxzi = dxi * dzi
dyzi = dyi * dzi
!
dxyzi = dxi * dyi * dzi
!
!
!------- loop over all points to be interpolated
!
do n = 1, kt
!
!------- determine location in (equidistant!!!) table
!
ix = 2 + INT( (x(n) - xt(1) - 1.e-10) * dxi )
iy = 2 + INT( (y(n) - yt(1) - 1.e-10) * dyi )
iz = 2 + INT( (z(n) - zt(1) - 1.e-10) * dzi )
!
ix = MAX( 2, MIN( ix, nx ) )
iy = MAX( 2, MIN( iy, ny ) )
iz = MAX( 2, MIN( iz, nz ) )
!
! write(*,*) iy-1,iy,iy+1
!
!------- set-up auxiliary arrays for Lagrange interpolation
!
delx(n) = xt(ix) - x(n)
dely(n) = yt(iy) - y(n)
delz(n) = zt(iz) - z(n)
!
do iv = 1, nvars
fh(n,1,iv) = ft(ix , iy , iz, iv )
fh(n,2,iv) = ft(ix-1, iy , iz, iv )
fh(n,3,iv) = ft(ix , iy-1, iz, iv )
fh(n,4,iv) = ft(ix , iy , iz-1, iv)
fh(n,5,iv) = ft(ix-1, iy-1, iz, iv )
fh(n,6,iv) = ft(ix-1, iy , iz-1, iv)
fh(n,7,iv) = ft(ix , iy-1, iz-1, iv)
fh(n,8,iv) = ft(ix-1, iy-1, iz-1, iv)
!
!------ set up coefficients of the interpolation polynomial and
! evaluate function values
!
a1(n,iv) = fh(n,1,iv)
a2(n,iv) = dxi * ( fh(n,2,iv) - fh(n,1,iv) )
a3(n,iv) = dyi * ( fh(n,3,iv) - fh(n,1,iv) )
a4(n,iv) = dzi * ( fh(n,4,iv) - fh(n,1,iv) )
a5(n,iv) = dxyi * ( fh(n,5,iv) - fh(n,2,iv) - fh(n,3,iv) + fh(n,1,iv) )
a6(n,iv) = dxzi * ( fh(n,6,iv) - fh(n,2,iv) - fh(n,4,iv) + fh(n,1,iv) )
a7(n,iv) = dyzi * ( fh(n,7,iv) - fh(n,3,iv) - fh(n,4,iv) + fh(n,1,iv) )
a8(n,iv) = dxyzi * ( fh(n,8,iv) - fh(n,1,iv) + fh(n,2,iv) + fh(n,3,iv) + &
fh(n,4,iv) - fh(n,5,iv) - fh(n,6,iv) - fh(n,7,iv) )
!
f(n,iv) = a1(n,iv) + a2(n,iv) * delx(n) &
+ a3(n,iv) * dely(n) &
+ a4(n,iv) * delz(n) &
+ a5(n,iv) * delx(n) * dely(n) &
+ a6(n,iv) * delx(n) * delz(n) &
+ a7(n,iv) * dely(n) * delz(n) &
+ a8(n,iv) * delx(n) * dely(n) * delz(n)
!
enddo
enddo
!
end SUBROUTINE intp3d_many
|
