( ) (

This section describes the trilinear interpolation algorithm using pseudo-Fortran code.

First, number the eight vertices of a hexahedron.

x1 = x(i ,j ,k ) x2 = x(i+1,j ,k ) x3 = x(i ,j+1,k ) x4 = x(i+1,j+1,k ) x5 = x(i ,j ,k+1) x6 = x(i+1,j ,k+1) x7 = x(i ,j+1,k+1) x8 = x(i+1,j+1,k+1)

Defined as below, notice that xa = x1 when a = −1 and that xa = x2 when a = 1 and similar observations for the rest of the equations.

xa = [(1-a)*x1 + (1+a)*x2] / 2 xb = [(1-a)*x3 + (1+a)*x4] / 2 xc = [(1-a)*x5 + (1+a)*x6] / 2 xd = [(1-a)*x7 + (1+a)*x8] / 2 xe = [(1-b)*xa + (1+b)*xb] / 2 xf = [(1-b)*xc + (1+b)*xd] / 2 x = [(1-g)*xe + (1+g)*xf] / 2

Take three steps of substitutions.

x = {(1-g) * [(1-b)*xa + (1+b)*xb] + (1+g) * [(1-b)*xc + (1+b)*xd]} / 4 x = {(1-g) * [(1-b) * ((1-a)*x1 + (1+a)*x2) + (1+b) * ((1-a)*x3 + (1+a)*x4)] + (1+g) * [(1-b) * ((1-a)*x5 + (1+a)*x6) + (1+b) * ((1-a)*x7 + (1+a)*x8)]} / 8 x = {(1-g)*(1-b)*(1-a)*x1 + (1-g)*(1-b)*(1+a)*x2 + (1-g)*(1+b)*(1-a)*x3 + (1-g)*(1+b)*(1+a)*x4 + (1+g)*(1-b)*(1-a)*x5 + (1+g)*(1-b)*(1+a)*x6 + (1+g)*(1+b)*(1-a)*x7 + (1+g)*(1+b)*(1+a)*x8} / 8

Regrouping

x = (x8 + x7 + x6 + x5 + x4 + x3 + x2 + x1) / 8 + a*(x8 - x7 + x6 - x5 + x4 - x3 + x2 - x1) / 8 + b*(x8 + x7 - x6 - x5 + x4 + x3 - x2 - x1) / 8 + g*(x8 + x7 + x6 + x5 - x4 - x3 - x2 - x1) / 8 + a*b*(x8 - x7 - x6 + x5 + x4 - x3 - x2 + x1) / 8 + a*g*(x8 - x7 + x6 - x5 - x4 + x3 - x2 + x1) / 8 + b*g*(x8 + x7 - x6 - x5 - x4 - x3 + x2 + x1) / 8 + a*b*g*(x8 - x7 - x6 + x5 - x4 + x3 + x2 - x1) / 8

Define

f0 = (x8 + x7 + x6 + x5 + x4 + x3 + x2 + x1) / 8 - x f1 = (x8 - x7 + x6 - x5 + x4 - x3 + x2 - x1) / 8 f2 = (x8 + x7 - x6 - x5 + x4 + x3 - x2 - x1) / 8 f3 = (x8 + x7 + x6 + x5 - x4 - x3 - x2 - x1) / 8 f4 = (x8 - x7 - x6 + x5 + x4 - x3 - x2 + x1) / 8 f5 = (x8 - x7 + x6 - x5 - x4 + x3 - x2 + x1) / 8 f6 = (x8 + x7 - x6 - x5 - x4 - x3 + x2 + x1) / 8 f7 = (x8 - x7 - x6 + x5 - x4 + x3 + x2 - x1) / 8

Then define

f = f0 + f1*a + f2*b + f3*g + f4*a*b + f5*a*g + f6*b*g + f7*a*b*g

Exactly similar relations can be derived for `y` and `z`.

g = g0 + g1*a + g2*b + g3*g + g4*a*b + g5*a*g + g6*b*g + g7*a*b*g h = h0 + h1*a + h2*b + h3*g + h4*a*b + h5*a*g + h6*b*g + h7*a*b*gwhere the

f = 0 g = 0 h = 0

So the three equations can be solved for the unknowns `a`,
`b` and `g`.
A solution can be obtained from Newton's method by letting

df = -f dg = -g dh = -hwhere

df = f1*da + f2*db + f3*dg + f4*b*da + f4*a*db f5*g*da + f5*a*dg f6*g*db + f6*b*dg f7*b*g*da + f7*a*g*db + f7*a*b*dg

Regrouping

df = (f1 + f4*b + f5*g + f7*b*g)*da + (f2 + f4*a + f6*g + f7*a*g)*db + (f3 + f5*a + f6*b + f7*a*b)*dg = -f

Similarly

dg = (g1 + g4*b + g5*g + g7*b*g)*da + (g2 + g4*a + g6*g + g7*a*g)*db + (g3 + g5*a + g6*b + g7*a*b)*dg = -g dh = (h1 + h4*b + h5*g + h7*b*g)*da + (h2 + h4*a + h6*g + h7*a*g)*db + (h3 + h5*a + h6*b + h7*a*b)*dg = -h

Using an initial guess of `a` = `b` = `g` = 0
the three equations can be solved for `da,` `db` and
`dg`.
Then one has

anew = aold + da bnew = bold + db gnew = gold + dg

If the matrix of the coefficients are zero it is because the hexahedron
is degenerate.
The new `a`, `b` and `g` can be used to compute new
coefficients and solve for new deltas and iteration should continue util
the deltas are negligible.
The interpolating program has a hard-wired limit of 20 iterations.
If after 20 iterations or if the absolute values of any of the three
parameters exceed 5.0, then the method fails and the starting stencil is
returned as the "best" one.

After `a`, `b` and `g` have been computed the
interplated value of any variable known at the vertices can be computed
by

v = [(1-g)*(1-b)*(1-a)*v1 + (1-g)*(1-b)*(1+a)*v2 + (1-g)*(1+b)*(1-a)*v3 + (1-g)*(1+b)*(1+a)*v4 + (1+g)*(1-b)*(1-a)*v5 + (1+g)*(1-b)*(1+a)*v6 + (1+g)*(1+b)*(1-a)*v7 + (1+g)*(1+b)*(1+a)*v8] / 8

This equation was copied from the equation for `X` and can easily
be verified by letting `a`, `b` and `g` take on
their extreme values of ±1.
If (`x,y,z`) is inside the hexahedron, then the parameters
`a`, `b` and `g` will be between ±1.
If one of them is outside those limits then the given point is outside
the hexahedron and the above formula becomes an extrapolation, not an
interpolation.