TricubicSplineInterpolatingFunction implements TrivariateFunction { public double value(double x, double y, double z) throws OutOfRangeException { final int i = searchIndex(x, xval); if (i == -1) { throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]); } final int j = searchIndex(y, yval); if (j == -1) { throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]); } final int k = searchIndex(z, zval); if (k == -1) { throw new OutOfRangeException(z, zval[0], zval[zval.length - 1]); } final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]); return splines[i][j][k].value(xN, yN, zN); } TricubicSplineInterpolatingFunction(double[] x,
                                               double[] y,
                                               double[] z,
                                               double[][][] f,
                                               double[][][] dFdX,
                                               double[][][] dFdY,
                                               double[][][] dFdZ,
                                               double[][][] d2FdXdY,
                                               double[][][] d2FdXdZ,
                                               double[][][] d2FdYdZ,
                                               double[][][] d3FdXdYdZ); double value(double x, double y, double z);  }
@Test public void testPlane() { double[] xval = new double[] {3, 4, 5, 6.5}; double[] yval = new double[] {-4, -3, -1, 2, 2.5}; double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5}; TrivariateFunction f = new TrivariateFunction() { public double value(double x, double y, double z) { return 2 * x - 3 * y - 4 * z + 5; } }; double[][][] fval = new double[xval.length][yval.length][zval.length]; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { fval[i][j][k] = f.value(xval[i], yval[j], zval[k]); } } } double[][][] dFdX = new double[xval.length][yval.length][zval.length]; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { dFdX[i][j][k] = 2; } } } double[][][] dFdY = new double[xval.length][yval.length][zval.length]; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { dFdY[i][j][k] = -3; } } } double[][][] dFdZ = new double[xval.length][yval.length][zval.length]; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { dFdZ[i][j][k] = -4; } } } double[][][] d2FdXdY = new double[xval.length][yval.length][zval.length]; double[][][] d2FdXdZ = new double[xval.length][yval.length][zval.length]; double[][][] d2FdYdZ = new double[xval.length][yval.length][zval.length]; double[][][] d3FdXdYdZ = new double[xval.length][yval.length][zval.length]; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { d2FdXdY[i][j][k] = 0; d2FdXdZ[i][j][k] = 0; d2FdYdZ[i][j][k] = 0; d3FdXdYdZ[i][j][k] = 0; } } } TrivariateFunction tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval, fval, dFdX, dFdY, dFdZ, d2FdXdY, d2FdXdZ, d2FdYdZ, d3FdXdYdZ); double x, y, z; double expected, result; x = 4; y = -3; z = 0; expected = f.value(x, y, z); result = tcf.value(x, y, z); Assert.assertEquals("On sample point", expected, result, 1e-15); x = 4.5; y = -1.5; z = -4.25; expected = f.value(x, y, z); result = tcf.value(x, y, z); Assert.assertEquals("Half-way between sample points (middle of the patch)", expected, result, 0.3); x = 3.5; y = -3.5; z = -10; expected = f.value(x, y, z); result = tcf.value(x, y, z); Assert.assertEquals("Half-way between sample points (border of the patch)", expected, result, 0.3); }
@Test public void testWave() { double[] xval = new double[] {3, 4, 5, 6.5}; double[] yval = new double[] {-4, -3, -1, 2, 2.5}; double[] zval = new double[] {-12, -8, -5.5, -3, 0, 4}; final double a = 0.2; final double omega = 0.5; final double kx = 2; final double ky = 1; TrivariateFunction f = new TrivariateFunction() { public double value(double x, double y, double z) { return a * FastMath.cos(omega * z - kx * x - ky * y); } }; double[][][] fval = new double[xval.length][yval.length][zval.length]; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { fval[i][j][k] = f.value(xval[i], yval[j], zval[k]); } } } double[][][] dFdX = new double[xval.length][yval.length][zval.length]; TrivariateFunction dFdX_f = new TrivariateFunction() { public double value(double x, double y, double z) { return a * FastMath.sin(omega * z - kx * x - ky * y) * kx; } }; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { dFdX[i][j][k] = dFdX_f.value(xval[i], yval[j], zval[k]); } } } double[][][] dFdY = new double[xval.length][yval.length][zval.length]; TrivariateFunction dFdY_f = new TrivariateFunction() { public double value(double x, double y, double z) { return a * FastMath.sin(omega * z - kx * x - ky * y) * ky; } }; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { dFdY[i][j][k] = dFdY_f.value(xval[i], yval[j], zval[k]); } } } double[][][] dFdZ = new double[xval.length][yval.length][zval.length]; TrivariateFunction dFdZ_f = new TrivariateFunction() { public double value(double x, double y, double z) { return -a * FastMath.sin(omega * z - kx * x - ky * y) * omega; } }; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { dFdZ[i][j][k] = dFdZ_f.value(xval[i], yval[j], zval[k]); } } } double[][][] d2FdXdY = new double[xval.length][yval.length][zval.length]; TrivariateFunction d2FdXdY_f = new TrivariateFunction() { public double value(double x, double y, double z) { return -a * FastMath.cos(omega * z - kx * x - ky * y) * kx * ky; } }; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { d2FdXdY[i][j][k] = d2FdXdY_f.value(xval[i], yval[j], zval[k]); } } } double[][][] d2FdXdZ = new double[xval.length][yval.length][zval.length]; TrivariateFunction d2FdXdZ_f = new TrivariateFunction() { public double value(double x, double y, double z) { return a * FastMath.cos(omega * z - kx * x - ky * y) * kx * omega; } }; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { d2FdXdZ[i][j][k] = d2FdXdZ_f.value(xval[i], yval[j], zval[k]); } } } double[][][] d2FdYdZ = new double[xval.length][yval.length][zval.length]; TrivariateFunction d2FdYdZ_f = new TrivariateFunction() { public double value(double x, double y, double z) { return a * FastMath.cos(omega * z - kx * x - ky * y) * ky * omega; } }; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { d2FdYdZ[i][j][k] = d2FdYdZ_f.value(xval[i], yval[j], zval[k]); } } } double[][][] d3FdXdYdZ = new double[xval.length][yval.length][zval.length]; TrivariateFunction d3FdXdYdZ_f = new TrivariateFunction() { public double value(double x, double y, double z) { return a * FastMath.sin(omega * z - kx * x - ky * y) * kx * ky * omega; } }; for (int i = 0; i < xval.length; i++) { for (int j = 0; j < yval.length; j++) { for (int k = 0; k < zval.length; k++) { d3FdXdYdZ[i][j][k] = d3FdXdYdZ_f.value(xval[i], yval[j], zval[k]); } } } TrivariateFunction tcf = new TricubicSplineInterpolatingFunction(xval, yval, zval, fval, dFdX, dFdY, dFdZ, d2FdXdY, d2FdXdZ, d2FdYdZ, d3FdXdYdZ); double x, y, z; double expected, result; x = 4; y = -3; z = 0; expected = f.value(x, y, z); result = tcf.value(x, y, z); Assert.assertEquals("On sample point", expected, result, 1e-14); x = 4.5; y = -1.5; z = -4.25; expected = f.value(x, y, z); result = tcf.value(x, y, z); Assert.assertEquals("Half-way between sample points (middle of the patch)", expected, result, 0.1); x = 3.5; y = -3.5; z = -10; expected = f.value(x, y, z); result = tcf.value(x, y, z); Assert.assertEquals("Half-way between sample points (border of the patch)", expected, result, 0.1); }