1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 package org.hipparchus.analysis.interpolation;
23
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.util.MathArrays;
27 import org.hipparchus.util.MathUtils;
28
29
30
31
32
33 public class TricubicInterpolator
34 implements TrivariateGridInterpolator {
35
36
37
38
39
40
41
42
43 public TricubicInterpolator() {
44
45 }
46
47
48
49
50 @Override
51 public TricubicInterpolatingFunction interpolate(final double[] xval,
52 final double[] yval,
53 final double[] zval,
54 final double[][][] fval)
55 throws MathIllegalArgumentException {
56 if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) {
57 throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
58 }
59 if (xval.length != fval.length) {
60 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
61 xval.length, fval.length);
62 }
63
64 MathArrays.checkOrder(xval);
65 MathArrays.checkOrder(yval);
66 MathArrays.checkOrder(zval);
67
68 final int xLen = xval.length;
69 final int yLen = yval.length;
70 final int zLen = zval.length;
71
72
73 final double[][][] dFdX = new double[xLen][yLen][zLen];
74 final double[][][] dFdY = new double[xLen][yLen][zLen];
75 final double[][][] dFdZ = new double[xLen][yLen][zLen];
76 final double[][][] d2FdXdY = new double[xLen][yLen][zLen];
77 final double[][][] d2FdXdZ = new double[xLen][yLen][zLen];
78 final double[][][] d2FdYdZ = new double[xLen][yLen][zLen];
79 final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen];
80
81 for (int i = 1; i < xLen - 1; i++) {
82 MathUtils.checkDimension(yval.length, fval[i].length);
83
84 final int nI = i + 1;
85 final int pI = i - 1;
86
87 final double nX = xval[nI];
88 final double pX = xval[pI];
89
90 final double deltaX = nX - pX;
91
92 for (int j = 1; j < yLen - 1; j++) {
93 MathUtils.checkDimension(zval.length, fval[i][j].length);
94
95 final int nJ = j + 1;
96 final int pJ = j - 1;
97
98 final double nY = yval[nJ];
99 final double pY = yval[pJ];
100
101 final double deltaY = nY - pY;
102 final double deltaXY = deltaX * deltaY;
103
104 for (int k = 1; k < zLen - 1; k++) {
105 final int nK = k + 1;
106 final int pK = k - 1;
107
108 final double nZ = zval[nK];
109 final double pZ = zval[pK];
110
111 final double deltaZ = nZ - pZ;
112
113 dFdX[i][j][k] = (fval[nI][j][k] - fval[pI][j][k]) / deltaX;
114 dFdY[i][j][k] = (fval[i][nJ][k] - fval[i][pJ][k]) / deltaY;
115 dFdZ[i][j][k] = (fval[i][j][nK] - fval[i][j][pK]) / deltaZ;
116
117 final double deltaXZ = deltaX * deltaZ;
118 final double deltaYZ = deltaY * deltaZ;
119
120 d2FdXdY[i][j][k] = (fval[nI][nJ][k] - fval[nI][pJ][k] - fval[pI][nJ][k] + fval[pI][pJ][k]) / deltaXY;
121 d2FdXdZ[i][j][k] = (fval[nI][j][nK] - fval[nI][j][pK] - fval[pI][j][nK] + fval[pI][j][pK]) / deltaXZ;
122 d2FdYdZ[i][j][k] = (fval[i][nJ][nK] - fval[i][nJ][pK] - fval[i][pJ][nK] + fval[i][pJ][pK]) / deltaYZ;
123
124 final double deltaXYZ = deltaXY * deltaZ;
125
126 d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] -
127 fval[pI][nJ][nK] + fval[pI][pJ][nK] -
128 fval[nI][nJ][pK] + fval[nI][pJ][pK] +
129 fval[pI][nJ][pK] - fval[pI][pJ][pK]) / deltaXYZ;
130 }
131 }
132 }
133
134
135 return new TricubicInterpolatingFunction(xval, yval, zval, fval,
136 dFdX, dFdY, dFdZ,
137 d2FdXdY, d2FdXdZ, d2FdYdZ,
138 d3FdXdYdZ) {
139
140 @Override
141 public boolean isValidPoint(double x, double y, double z) {
142 if (x < xval[1] ||
143 x > xval[xval.length - 2] ||
144 y < yval[1] ||
145 y > yval[yval.length - 2] ||
146 z < zval[1] ||
147 z > zval[zval.length - 2]) {
148 return false;
149 } else {
150 return true;
151 }
152 }
153 };
154 }
155 }