1 | /**********************************************************
|
---|
2 | *
|
---|
3 | * TriCubicInterpolation.java
|
---|
4 | *
|
---|
5 | * Class for performing an interpolation on the tabulated
|
---|
6 | * function y = f(x1,x2,x3) using a tricubic interploation procedure
|
---|
7 | **
|
---|
8 | * WRITTEN BY: Dr Michael Thomas Flanagan
|
---|
9 | *
|
---|
10 | * DATE: 12-15 January 2011
|
---|
11 | * UPDATE:
|
---|
12 | *
|
---|
13 | * DOCUMENTATION:
|
---|
14 | * See Michael Thomas Flanagan's Java library on-line web page:
|
---|
15 | * http://www.ee.ucl.ac.uk/~mflanaga/java/TriCubicInterpolation.html
|
---|
16 | * http://www.ee.ucl.ac.uk/~mflanaga/java/
|
---|
17 | *
|
---|
18 | * Copyright (c) 2011 Michael Thomas Flanagan
|
---|
19 | *
|
---|
20 | * PERMISSION TO COPY:
|
---|
21 | *
|
---|
22 | * Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee,
|
---|
23 | * provided that an acknowledgement to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies
|
---|
24 | * and associated documentation or publications.
|
---|
25 | *
|
---|
26 | * Redistributions of the source code of this source code, or parts of the source codes, must retain the above copyright notice, this list of conditions
|
---|
27 | * and the following disclaimer and requires written permission from the Michael Thomas Flanagan:
|
---|
28 | *
|
---|
29 | * Redistribution in binary form of all or parts of this class must reproduce the above copyright notice, this list of conditions and
|
---|
30 | * the following disclaimer in the documentation and/or other materials provided with the distribution and requires written permission from the Michael Thomas Flanagan:
|
---|
31 | *
|
---|
32 | * Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose.
|
---|
33 | * Dr Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software
|
---|
34 | * or its derivatives.
|
---|
35 | *
|
---|
36 | ***************************************************************************************/
|
---|
37 |
|
---|
38 | package agents.anac.y2015.agentBuyogV2.flanagan.interpolation;
|
---|
39 |
|
---|
40 | import java.util.ArrayList;
|
---|
41 |
|
---|
42 | import agents.anac.y2015.agentBuyogV2.flanagan.math.ArrayMaths;
|
---|
43 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
|
---|
44 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
|
---|
45 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Matrix;
|
---|
46 |
|
---|
47 | public class TriCubicInterpolation{
|
---|
48 |
|
---|
49 | // unit cube: front face anticlockwise then backface anticlockwise from bottom left corner
|
---|
50 | int[][] unitCube = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}};
|
---|
51 |
|
---|
52 | private int lPoints = 0; // no. of x1 tabulated points
|
---|
53 | private int mPoints = 0; // no. of x2 tabulated points
|
---|
54 | private int nPoints = 0; // no. of x3 tabulated points
|
---|
55 | private double[] x1 = null; // x1 in tabulated function f(x1,x2,x3)
|
---|
56 | private double[] x2 = null; // x2 in tabulated function f(x1,x2,x3)
|
---|
57 | private double[] x3 = null; // x3 in tabulated function f(x1,x2,x3)
|
---|
58 | private double[][][] y = null; // y=f(x1,x2,x3) tabulated function
|
---|
59 |
|
---|
60 | private double[][][] dydx1 = null; // dy/dx1
|
---|
61 | private double[][][] dydx2 = null; // dy/dx2
|
---|
62 | private double[][][] dydx3 = null; // dy/dx3
|
---|
63 | private double[][][] d2ydx1dx2 = null; // d2y/dx1dx2
|
---|
64 | private double[][][] d2ydx1dx3 = null; // d2y/dx1dx3
|
---|
65 | private double[][][] d2ydx2dx3 = null; // d2y/dx2dx3
|
---|
66 | private double[][][] d3ydx1dx2dx3 = null; // d3y/dx1dx2dx3
|
---|
67 | private boolean derivCalculated = false; // = true when the derivatives have been calculated or entered
|
---|
68 | private TriCubicSpline tcs = null; // TriCubic spline used in calculating the derivatives
|
---|
69 | private double incrX1 = 0; // x1 increment used in calculating the derivatives
|
---|
70 | private double incrX2 = 0; // x2 increment used in calculating the derivatives
|
---|
71 | private double incrX3 = 0; // x3 increment used in calculating the derivatives
|
---|
72 |
|
---|
73 | private double xx1 = Double.NaN; // value of x1 at which an interpolated y value is required
|
---|
74 | private double xx2 = Double.NaN; // value of x2 at which an interpolated y value is required
|
---|
75 | private double xx3 = Double.NaN; // value of x3 at which an interpolated y value is required
|
---|
76 |
|
---|
77 | private ArrayList<Object> coeff = new ArrayList<Object>(); // grid cube coefficients
|
---|
78 |
|
---|
79 | // Weights used in calculating the grid cube coefficients
|
---|
80 | private double[][] weights = new double[64][64];
|
---|
81 |
|
---|
82 | private int[] x1indices = null; // x1 data indices before ordering
|
---|
83 | private int[] x2indices = null; // x2 data indices before ordering
|
---|
84 | private int[] x3indices = null; // x3 data indices before ordering
|
---|
85 |
|
---|
86 | private double[] xMin = new double[3]; // minimum values of x1, x2 and x3
|
---|
87 | private double[] xMax = new double[3]; // maximum values of x1, x2 and x3
|
---|
88 |
|
---|
89 | private double interpolatedValue = Double.NaN; // interpolated value of y
|
---|
90 | private double interpolatedDydx1 = Double.NaN; // interpolated value of dydx1
|
---|
91 | private double interpolatedDydx2 = Double.NaN; // interpolated value of dydx2
|
---|
92 | private double interpolatedDydx3 = Double.NaN; // interpolated value of dydx3
|
---|
93 | private double interpolatedD2ydx1dx2 = Double.NaN; // interpolated value of d2ydx1dx2
|
---|
94 | private double interpolatedD2ydx1dx3 = Double.NaN; // interpolated value of d2ydx1dx3
|
---|
95 | private double interpolatedD2ydx2dx3 = Double.NaN; // interpolated value of d2ydx2dx3
|
---|
96 | private double interpolatedD3ydx1dx2dx3 = Double.NaN; // interpolated value of d3ydx1dx2d3
|
---|
97 |
|
---|
98 | private boolean numerDiffFlag = true; // = true: if numerical differentiation performed h1 and h2 calculated using delta
|
---|
99 | // = false: if numerical differentiation performed h1 and h2 calculated only provided data points
|
---|
100 |
|
---|
101 | private static double delta = 1e-3; // fractional step factor used in calculating the derivatives
|
---|
102 | private static double potentialRoundingError = 5e-15; // potential rounding error used in checking wheter a value lies within the interpolation bounds (static value)
|
---|
103 | private static boolean roundingCheck = false; // = true: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit (static value)
|
---|
104 |
|
---|
105 |
|
---|
106 | // Constructor without derivatives
|
---|
107 | // numerDiffOption = 0 -> numerical differencing using only supplied data points
|
---|
108 | // numerDiffOption = 1 -> numerical differencing using interpolation
|
---|
109 | public TriCubicInterpolation(double[] x1, double[] x2, double[] x3, double[][][] y, int numerDiffOption){
|
---|
110 | // set numerical differencing option
|
---|
111 | if(numerDiffOption==0){
|
---|
112 | this.numerDiffFlag = false;
|
---|
113 | }
|
---|
114 | else{
|
---|
115 | if(numerDiffOption==1){
|
---|
116 | this.numerDiffFlag = true;
|
---|
117 | }
|
---|
118 | else{
|
---|
119 | throw new IllegalArgumentException("The numerical differencing option, " + numerDiffOption + ", must be 0 or 1");
|
---|
120 | }
|
---|
121 | }
|
---|
122 |
|
---|
123 | // initialize the data
|
---|
124 | this.initialize(Conv.copy(x1), Conv.copy(x2), Conv.copy(x3), Conv.copy(y));
|
---|
125 |
|
---|
126 | // calculate the derivatives
|
---|
127 | this.calcDeriv();
|
---|
128 |
|
---|
129 | // calculate grid coefficients for all grid cubes
|
---|
130 | this.gridCoefficients();
|
---|
131 | }
|
---|
132 |
|
---|
133 | // Constructor with derivatives
|
---|
134 | public TriCubicInterpolation(double[] x1, double[] x2, double[] x3, double[][][] y, double[][][] dydx1, double[][][] dydx2, double[][][] dydx3, double[][][] d2ydx1dx2, double[][][] d2ydx1dx3, double[][][] d2ydx2dx3, double[][][] d3ydx1dx2dx3){
|
---|
135 |
|
---|
136 | // initialize the data
|
---|
137 | this.initialize(Conv.copy(x1), Conv.copy(x2), Conv.copy(x3), Conv.copy(y), Conv.copy(dydx1), Conv.copy(dydx2), Conv.copy(dydx3), Conv.copy(d2ydx1dx2), Conv.copy(d2ydx1dx3), Conv.copy(d2ydx2dx3), Conv.copy(d3ydx1dx2dx3));
|
---|
138 |
|
---|
139 | // calculate grid coefficients for all grid cubes
|
---|
140 | this.gridCoefficients();
|
---|
141 | }
|
---|
142 |
|
---|
143 | // Initialize the data
|
---|
144 | private void initialize(double[] x1, double[] x2, double[] x3, double[][][] y){
|
---|
145 | this.initialize(x1, x2, x3, y, null, null, null, null, null, null, null, false);
|
---|
146 | }
|
---|
147 |
|
---|
148 | private void initialize(double[] x1, double[] x2, double[] x3, double[][][] y, double[][][] dydx1, double[][][] dydx2, double[][][] dyd3, double[][][] d2ydx1dx2, double[][][] d2ydx1dx3, double[][][] d2ydx2dx3, double[][][] d3ydx1dx2dx3){
|
---|
149 |
|
---|
150 | this.initialize(x1, x2, x3, y, dydx1, dydx2, dydx3, d2ydx1dx2, d2ydx1dx3, d2ydx2dx3, d3ydx1dx2dx3, true);
|
---|
151 | }
|
---|
152 |
|
---|
153 | private void initialize(double[] x1, double[] x2, double[] x3, double[][][] y, double[][][] dydx1, double[][][] dydx2, double[][][] dyd3, double[][][] d2ydx1dx2, double[][][] d2ydx1dx3, double[][][] d2ydx2dx3, double[][][] d3ydx1dx2dx3, boolean flag){
|
---|
154 | int lPoints=x1.length;
|
---|
155 | int mPoints=x2.length;
|
---|
156 | int nPoints=x3.length;
|
---|
157 | if(lPoints!=y.length)throw new IllegalArgumentException("Array x1 and y-row are of different length " + lPoints + " " + y.length);
|
---|
158 | if(mPoints!=y[0].length)throw new IllegalArgumentException("Array x2 and y-column are of different length "+ mPoints + " " + y[0].length);
|
---|
159 | if(nPoints!=y[0][0].length)throw new IllegalArgumentException("Array x3 and y-column are of different length "+ nPoints + " " + y[0][0].length);
|
---|
160 | if(lPoints<2 || mPoints<2 || nPoints<2 )throw new IllegalArgumentException("The data matrix must have a minimum size of 2 X 2 X 2");
|
---|
161 |
|
---|
162 | // Calculate weighting matrix
|
---|
163 | this.calcWeights();
|
---|
164 |
|
---|
165 | // order data
|
---|
166 | ArrayMaths am = new ArrayMaths(x1);
|
---|
167 | am = am.sort();
|
---|
168 | this.x1indices = am.originalIndices();
|
---|
169 | x1 = am.array();
|
---|
170 | double[][][] hold = new double[lPoints][mPoints][nPoints];
|
---|
171 | double[][][] hold1 = null;
|
---|
172 | double[][][] hold2 = null;
|
---|
173 | double[][][] hold12 = null;
|
---|
174 |
|
---|
175 | for(int i=0; i<lPoints; i++){
|
---|
176 | for(int j=0; j<mPoints; j++){
|
---|
177 | for(int k=0; k<nPoints; k++){
|
---|
178 | hold[i][j][k] = y[this.x1indices[i]][j][k];
|
---|
179 | }
|
---|
180 | }
|
---|
181 | }
|
---|
182 | for(int i=0; i<lPoints; i++){
|
---|
183 | for(int j=0; j<mPoints; j++){
|
---|
184 | for(int k=0; k<nPoints; k++){
|
---|
185 | y[i][j][k] = hold[i][j][k];
|
---|
186 | }
|
---|
187 | }
|
---|
188 | }
|
---|
189 |
|
---|
190 | if(flag){
|
---|
191 | hold1 = new double[lPoints][mPoints][nPoints];
|
---|
192 | hold2 = new double[lPoints][mPoints][nPoints];
|
---|
193 | hold12 = new double[lPoints][mPoints][nPoints];
|
---|
194 | for(int i=0; i<lPoints; i++){
|
---|
195 | for(int j=0; j<mPoints; j++){
|
---|
196 | for(int k=0; k<nPoints; k++){
|
---|
197 | hold1[i][j][k] = dydx1[this.x1indices[i]][j][k];
|
---|
198 | hold2[i][j][k] = dydx2[this.x1indices[i]][j][k];
|
---|
199 | hold12[i][j][k] = d2ydx1dx2[this.x1indices[i]][j][k];
|
---|
200 | }
|
---|
201 | }
|
---|
202 | }
|
---|
203 | for(int i=0; i<lPoints; i++){
|
---|
204 | for(int j=0; j<mPoints; j++){
|
---|
205 | for(int k=0; k<nPoints; k++){
|
---|
206 | dydx1[i][j][k] = hold1[i][j][k];
|
---|
207 | dydx2[i][j][k] = hold2[i][j][k];
|
---|
208 | d2ydx1dx2[i][j][k] = hold12[i][j][k];
|
---|
209 | }
|
---|
210 | }
|
---|
211 | }
|
---|
212 | }
|
---|
213 |
|
---|
214 | am = new ArrayMaths(x2);
|
---|
215 | am = am.sort();
|
---|
216 | this.x2indices = am.originalIndices();
|
---|
217 | x2 = am.array();
|
---|
218 |
|
---|
219 | for(int i=0; i<lPoints; i++){
|
---|
220 | for(int j=0; j<mPoints; j++){
|
---|
221 | for(int k=0; k<nPoints; k++){
|
---|
222 | hold[i][j][k] = y[i][this.x2indices[j]][k];
|
---|
223 | }
|
---|
224 | }
|
---|
225 | }
|
---|
226 | for(int i=0; i<lPoints; i++){
|
---|
227 | for(int j=0; j<mPoints; j++){
|
---|
228 | for(int k=0; k<nPoints; k++){
|
---|
229 | y[i][j][k] = hold[i][j][k];
|
---|
230 | }
|
---|
231 | }
|
---|
232 | }
|
---|
233 |
|
---|
234 | if(flag){
|
---|
235 | for(int i=0; i<lPoints; i++){
|
---|
236 | for(int j=0; j<mPoints; j++){
|
---|
237 | for(int k=0; k<nPoints; k++){
|
---|
238 | hold1[i][j][k] = dydx1[i][this.x2indices[j]][k];
|
---|
239 | hold2[i][j][k] = dydx2[i][this.x2indices[j]][k];
|
---|
240 | hold12[i][j][k] = d2ydx1dx2[i][this.x2indices[j]][k];
|
---|
241 | }
|
---|
242 | }
|
---|
243 | }
|
---|
244 | for(int i=0; i<lPoints; i++){
|
---|
245 | for(int j=0; j<mPoints; j++){
|
---|
246 | for(int k=0; k<nPoints; k++){
|
---|
247 | dydx1[i][j][k] = hold1[i][j][k];
|
---|
248 | dydx2[i][j][k] = hold2[i][j][k];
|
---|
249 | d2ydx1dx2[i][j][k] = hold12[i][j][k];
|
---|
250 | }
|
---|
251 | }
|
---|
252 | }
|
---|
253 | }
|
---|
254 |
|
---|
255 | am = new ArrayMaths(x3);
|
---|
256 | am = am.sort();
|
---|
257 | this.x3indices = am.originalIndices();
|
---|
258 | x3 = am.array();
|
---|
259 |
|
---|
260 | for(int i=0; i<lPoints; i++){
|
---|
261 | for(int j=0; j<mPoints; j++){
|
---|
262 | for(int k=0; k<nPoints; k++){
|
---|
263 | hold[i][j][k] = y[i][j][this.x3indices[k]];
|
---|
264 | }
|
---|
265 | }
|
---|
266 | }
|
---|
267 | for(int i=0; i<lPoints; i++){
|
---|
268 | for(int j=0; j<mPoints; j++){
|
---|
269 | for(int k=0; k<nPoints; k++){
|
---|
270 | y[i][j][k] = hold[i][j][k];
|
---|
271 | }
|
---|
272 | }
|
---|
273 | }
|
---|
274 |
|
---|
275 | if(flag){
|
---|
276 | for(int i=0; i<lPoints; i++){
|
---|
277 | for(int j=0; j<mPoints; j++){
|
---|
278 | for(int k=0; k<nPoints; k++){
|
---|
279 | hold1[i][j][k] = dydx1[i][j][this.x3indices[k]];
|
---|
280 | hold2[i][j][k] = dydx2[i][j][this.x3indices[k]];
|
---|
281 | hold12[i][j][k] = d2ydx1dx2[i][j][this.x3indices[k]];
|
---|
282 | }
|
---|
283 | }
|
---|
284 | }
|
---|
285 | for(int i=0; i<lPoints; i++){
|
---|
286 | for(int j=0; j<mPoints; j++){
|
---|
287 | for(int k=0; k<nPoints; k++){
|
---|
288 | dydx1[i][j][k] = hold1[i][j][k];
|
---|
289 | dydx2[i][j][k] = hold2[i][j][k];
|
---|
290 | d2ydx1dx2[i][j][k] = hold12[i][j][k];
|
---|
291 | }
|
---|
292 | }
|
---|
293 | }
|
---|
294 | }
|
---|
295 |
|
---|
296 | // check for identical x1 values
|
---|
297 | for(int i=1; i<lPoints; i++){
|
---|
298 | if(x1[i]==x1[i-1]){
|
---|
299 | System.out.println("x1["+this.x1indices[i]+"] and x1["+this.x1indices[i+1]+"] are identical, " + x1[i]);
|
---|
300 | double sep = (Fmath.maximum(x1) - Fmath.minimum(x1))/0.5e-3;
|
---|
301 | x1[i-1] -= sep;
|
---|
302 | x1[i]+= sep;
|
---|
303 | System.out.println("They have been separated by" + 2*sep);
|
---|
304 | }
|
---|
305 | }
|
---|
306 |
|
---|
307 | // check for identical x2 values
|
---|
308 | for(int i=1; i<mPoints; i++){
|
---|
309 | if(x2[i]==x2[i-1]){
|
---|
310 | System.out.println("x2["+this.x2indices[i]+"] and x2["+this.x2indices[i+1]+"] are identical, " + x2[i]);
|
---|
311 | double sep = (Fmath.maximum(x2) - Fmath.minimum(x2))/0.5e-3;
|
---|
312 | x2[i-1] -= sep;
|
---|
313 | x2[i]+= sep;
|
---|
314 | System.out.println("They have been separated by" + 2*sep);
|
---|
315 | }
|
---|
316 | }
|
---|
317 |
|
---|
318 | // check for identical x3 values
|
---|
319 | for(int i=1; i<nPoints; i++){
|
---|
320 | if(x3[i]==x3[i-1]){
|
---|
321 | System.out.println("x3["+this.x3indices[i]+"] and x3["+this.x3indices[i+1]+"] are identical, " + x3[i]);
|
---|
322 | double sep = (Fmath.maximum(x3) - Fmath.minimum(x3))/0.5e-3;
|
---|
323 | x3[i-1] -= sep;
|
---|
324 | x3[i]+= sep;
|
---|
325 | System.out.println("They have been separated by" + 2*sep);
|
---|
326 | }
|
---|
327 | }
|
---|
328 |
|
---|
329 | // assign variables
|
---|
330 | this.lPoints = lPoints;
|
---|
331 | this.mPoints = mPoints;
|
---|
332 | this.nPoints = nPoints;
|
---|
333 | this.x1 = new double[this.lPoints];
|
---|
334 | this.x2 = new double[this.mPoints];
|
---|
335 | this.x3 = new double[this.nPoints];
|
---|
336 | this.y = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
337 | this.dydx1 = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
338 | this.dydx2 = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
339 | this.dydx3 = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
340 | this.d2ydx1dx2 = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
341 | this.d2ydx1dx3 = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
342 | this.d2ydx2dx3 = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
343 | this.d3ydx1dx2dx3 = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
344 |
|
---|
345 | for(int i=0; i<this.lPoints; i++){
|
---|
346 | this.x1[i]=x1[i];
|
---|
347 | }
|
---|
348 | for(int j=0; j<this.mPoints; j++){
|
---|
349 | this.x2[j]=x2[j];
|
---|
350 | }
|
---|
351 | for(int k=0; k<this.nPoints; k++){
|
---|
352 | this.x3[k]=x3[k];
|
---|
353 | }
|
---|
354 | for(int i =0; i<this.lPoints; i++){
|
---|
355 | for(int j=0; j<this.mPoints; j++){
|
---|
356 | for(int k=0; k<this.nPoints; k++){
|
---|
357 | this.y[i][j][k]=y[i][j][k];
|
---|
358 | }
|
---|
359 | }
|
---|
360 | }
|
---|
361 |
|
---|
362 | if(flag){
|
---|
363 | for(int i =0; i<this.lPoints; i++){
|
---|
364 | for(int j=0; j<this.mPoints; j++){
|
---|
365 | for(int k=0; k<this.nPoints; k++){
|
---|
366 | this.dydx1[i][j][k]=dydx1[i][j][k];
|
---|
367 | this.dydx2[i][j][k]=dydx2[i][j][k];
|
---|
368 | this.dydx3[i][j][k]=dydx3[i][j][k];
|
---|
369 | this.d2ydx1dx2[i][j][k]=d2ydx1dx2[i][j][k];
|
---|
370 | this.d2ydx1dx3[i][j][k]=d2ydx1dx3[i][j][k];
|
---|
371 | this.d2ydx2dx3[i][j][k]=d2ydx2dx3[i][j][k];
|
---|
372 | this.d3ydx1dx2dx3[i][j][k]=d3ydx1dx2dx3[i][j][k];
|
---|
373 | }
|
---|
374 | }
|
---|
375 | }
|
---|
376 | this.derivCalculated = true;
|
---|
377 | }
|
---|
378 |
|
---|
379 | // limits
|
---|
380 | this.xMin[0] = Fmath.minimum(this.x1);
|
---|
381 | this.xMax[0] = Fmath.maximum(this.x1);
|
---|
382 | this.xMin[1] = Fmath.minimum(this.x2);
|
---|
383 | this.xMax[1] = Fmath.maximum(this.x2);
|
---|
384 | this.xMin[2] = Fmath.minimum(this.x3);
|
---|
385 | this.xMax[2] = Fmath.maximum(this.x3);
|
---|
386 |
|
---|
387 |
|
---|
388 | if(!flag && this.numerDiffFlag){
|
---|
389 | // numerical difference increments
|
---|
390 | double range1 = this.xMax[0] - this.xMin[0];
|
---|
391 | double range2 = this.xMax[1] - this.xMin[1];
|
---|
392 | double range3 = this.xMax[2] - this.xMin[2];
|
---|
393 | double averageSeparation1 = range1/this.lPoints;
|
---|
394 | double averageSeparation2 = range2/this.mPoints;
|
---|
395 | double averageSeparation3 = range3/this.nPoints;
|
---|
396 |
|
---|
397 | double minSep = this.x1[1] - this.x1[0];
|
---|
398 | double minimumSeparation1 = minSep;
|
---|
399 | for(int i=2; i<this.lPoints; i++){
|
---|
400 | minSep = this.x1[i] - this.x1[i-1];
|
---|
401 | if(minSep<minimumSeparation1)minimumSeparation1 = minSep;
|
---|
402 | }
|
---|
403 | minSep = this.x2[1] - this.x2[0];
|
---|
404 | double minimumSeparation2 = minSep;
|
---|
405 | for(int i=2; i<this.mPoints; i++){
|
---|
406 | minSep = this.x2[i] - this.x2[i-1];
|
---|
407 | if(minSep<minimumSeparation2)minimumSeparation2 = minSep;
|
---|
408 | }
|
---|
409 | minSep = this.x3[1] - this.x3[0];
|
---|
410 | double minimumSeparation3 = minSep;
|
---|
411 | for(int i=2; i<this.nPoints; i++){
|
---|
412 | minSep = this.x3[i] - this.x3[i-1];
|
---|
413 | if(minSep<minimumSeparation3)minimumSeparation3 = minSep;
|
---|
414 | }
|
---|
415 |
|
---|
416 | this.incrX1 = range1*TriCubicInterpolation.delta;
|
---|
417 | double defaultIncr = minimumSeparation1;
|
---|
418 | if(minimumSeparation1<averageSeparation1/10.0)defaultIncr = averageSeparation1/10.0;
|
---|
419 | if(this.incrX1>averageSeparation1)this.incrX1 = defaultIncr;
|
---|
420 | this.incrX2 = range2*TriCubicInterpolation.delta;
|
---|
421 | defaultIncr = minimumSeparation2;
|
---|
422 | if(minimumSeparation2<averageSeparation2/10.0)defaultIncr = averageSeparation2/10.0;
|
---|
423 | if(this.incrX2>averageSeparation2)this.incrX2 = defaultIncr;
|
---|
424 | this.incrX3 = range3*TriCubicInterpolation.delta;
|
---|
425 | defaultIncr = minimumSeparation3;
|
---|
426 | if(minimumSeparation3<averageSeparation3/10.0)defaultIncr = averageSeparation3/10.0;
|
---|
427 | if(this.incrX3>averageSeparation3)this.incrX3 = defaultIncr;
|
---|
428 | }
|
---|
429 | }
|
---|
430 |
|
---|
431 |
|
---|
432 | // Calculate the weighting matrix
|
---|
433 | private void calcWeights(){
|
---|
434 | int kk = 0;
|
---|
435 | // substitute unit cube corners in y
|
---|
436 | for(int m=0; m<8; m++){
|
---|
437 | int n = 0;
|
---|
438 | for(int i=0; i<4; i++){
|
---|
439 | for(int j=0; j<4; j++){
|
---|
440 | for(int k=0; k<4; k++){
|
---|
441 | this.weights[kk][n] = Math.pow(this.unitCube[m][0], i)* Math.pow(this.unitCube[m][1], j)*Math.pow(this.unitCube[m][2],k);
|
---|
442 | n++;
|
---|
443 | }
|
---|
444 | }
|
---|
445 | }
|
---|
446 | kk++;
|
---|
447 | }
|
---|
448 | // substitute unit cube corners in dy/dx1
|
---|
449 | for(int m=0; m<8; m++){
|
---|
450 | int n = 0;
|
---|
451 | for(int i=0; i<4; i++){
|
---|
452 | for(int j=0; j<4; j++){
|
---|
453 | for(int k=0; k<4; k++){
|
---|
454 | if(i==0){
|
---|
455 | this.weights[kk][n] = 0.0;
|
---|
456 | }
|
---|
457 | else{
|
---|
458 | this.weights[kk][n] = i*Math.pow(this.unitCube[m][0], i-1)* Math.pow(this.unitCube[m][1], j)*Math.pow(this.unitCube[m][2], k);
|
---|
459 | }
|
---|
460 | n++;
|
---|
461 | }
|
---|
462 | }
|
---|
463 | }
|
---|
464 | kk++;
|
---|
465 | }
|
---|
466 | // substitute unit cube corners in dy/dx2
|
---|
467 | for(int m=0; m<8; m++){
|
---|
468 | int n = 0;
|
---|
469 | for(int i=0; i<4; i++){
|
---|
470 | for(int j=0; j<4; j++){
|
---|
471 | for(int k=0; k<4; k++){
|
---|
472 | if(j==0){
|
---|
473 | this.weights[kk][n] = 0.0;
|
---|
474 | }
|
---|
475 | else{
|
---|
476 | this.weights[kk][n] = j*Math.pow(this.unitCube[m][0], i)* Math.pow(this.unitCube[m][1], j-1)*Math.pow(this.unitCube[m][2], k);
|
---|
477 | }
|
---|
478 | n++;
|
---|
479 | }
|
---|
480 | }
|
---|
481 | }
|
---|
482 | kk++;
|
---|
483 | }
|
---|
484 | // substitute unit cube corners in dy/dx3
|
---|
485 | for(int m=0; m<8; m++){
|
---|
486 | int n = 0;
|
---|
487 | for(int i=0; i<4; i++){
|
---|
488 | for(int j=0; j<4; j++){
|
---|
489 | for(int k=0; k<4; k++){
|
---|
490 | if(k==0){
|
---|
491 | this.weights[kk][n] = 0.0;
|
---|
492 | }
|
---|
493 | else{
|
---|
494 | this.weights[kk][n] = k*Math.pow(this.unitCube[m][0], i)* Math.pow(this.unitCube[m][1], j)*Math.pow(this.unitCube[m][2], k-1);
|
---|
495 | }
|
---|
496 | n++;
|
---|
497 | }
|
---|
498 | }
|
---|
499 | }
|
---|
500 | kk++;
|
---|
501 | }
|
---|
502 | // substitute unit cube corners in d2y/dx1dx2
|
---|
503 | for(int m=0; m<8; m++){
|
---|
504 | int n = 0;
|
---|
505 | for(int i=0; i<4; i++){
|
---|
506 | for(int j=0; j<4; j++){
|
---|
507 | for(int k=0; k<4; k++){
|
---|
508 | if(i==0 || j==0){
|
---|
509 | this.weights[kk][n] = 0.0;
|
---|
510 | }
|
---|
511 | else{
|
---|
512 | this.weights[kk][n] = i*j*Math.pow(this.unitCube[m][0], i-1)* Math.pow(this.unitCube[m][1], j-1)* Math.pow(this.unitCube[m][2], k);
|
---|
513 | }
|
---|
514 | n++;
|
---|
515 | }
|
---|
516 | }
|
---|
517 | }
|
---|
518 | kk++;
|
---|
519 | }
|
---|
520 | // substitute unit cube corners in d2y/dx1dx3
|
---|
521 | for(int m=0; m<8; m++){
|
---|
522 | int n = 0;
|
---|
523 | for(int i=0; i<4; i++){
|
---|
524 | for(int j=0; j<4; j++){
|
---|
525 | for(int k=0; k<4; k++){
|
---|
526 | if(i==0 || k==0){
|
---|
527 | this.weights[kk][n] = 0.0;
|
---|
528 | }
|
---|
529 | else{
|
---|
530 | this.weights[kk][n] = i*k*Math.pow(this.unitCube[m][0], i-1)* Math.pow(this.unitCube[m][1], j)* Math.pow(this.unitCube[m][2], k-1);
|
---|
531 | }
|
---|
532 | n++;
|
---|
533 | }
|
---|
534 | }
|
---|
535 | }
|
---|
536 | kk++;
|
---|
537 | }
|
---|
538 | // substitute unit cube corners in d2y/dx2dx3
|
---|
539 | for(int m=0; m<8; m++){
|
---|
540 | int n = 0;
|
---|
541 | for(int i=0; i<4; i++){
|
---|
542 | for(int j=0; j<4; j++){
|
---|
543 | for(int k=0; k<4; k++){
|
---|
544 | if(j==0 || k==0){
|
---|
545 | this.weights[kk][n] = 0.0;
|
---|
546 | }
|
---|
547 | else{
|
---|
548 | this.weights[kk][n] = j*k*Math.pow(this.unitCube[m][0], i)* Math.pow(this.unitCube[m][1], j-1)* Math.pow(this.unitCube[m][2], k-1);
|
---|
549 | }
|
---|
550 | n++;
|
---|
551 | }
|
---|
552 | }
|
---|
553 | }
|
---|
554 | kk++;
|
---|
555 | }
|
---|
556 | // substitute unit cube corners in d3y/dx1dx2dx3
|
---|
557 | for(int m=0; m<8; m++){
|
---|
558 | int n = 0;
|
---|
559 | for(int i=0; i<4; i++){
|
---|
560 | for(int j=0; j<4; j++){
|
---|
561 | for(int k=0; k<4; k++){
|
---|
562 | if(i==0 || j==0 || k==0){
|
---|
563 | this.weights[kk][n] = 0.0;
|
---|
564 | }
|
---|
565 | else{
|
---|
566 | this.weights[kk][n] = i*j*k*Math.pow(this.unitCube[m][0], i-1)* Math.pow(this.unitCube[m][1], j-1)* Math.pow(this.unitCube[m][2], k-1);
|
---|
567 | }
|
---|
568 | n++;
|
---|
569 | }
|
---|
570 | }
|
---|
571 | }
|
---|
572 | kk++;
|
---|
573 | }
|
---|
574 |
|
---|
575 | // invert the above calculated matrix
|
---|
576 | Matrix mat = new Matrix(this.weights);
|
---|
577 | mat = mat.inverse();
|
---|
578 | this.weights = mat.getArrayCopy();
|
---|
579 | }
|
---|
580 |
|
---|
581 |
|
---|
582 | // Calculate the derivatives
|
---|
583 | private void calcDeriv(){
|
---|
584 |
|
---|
585 | if(this.numerDiffFlag){
|
---|
586 |
|
---|
587 | // Numerical differentiation using delta and interpolation
|
---|
588 | this.tcs = new TriCubicSpline(this.x1, this.x2, this.x3, this.y);
|
---|
589 |
|
---|
590 | double[] x1jp1 = new double[this.lPoints];
|
---|
591 | double[] x1jm1 = new double[this.lPoints];
|
---|
592 | double[] x2jp1 = new double[this.mPoints];
|
---|
593 | double[] x2jm1 = new double[this.mPoints];
|
---|
594 | double[] x3jp1 = new double[this.nPoints];
|
---|
595 | double[] x3jm1 = new double[this.nPoints];
|
---|
596 |
|
---|
597 | for(int i=0; i<this.lPoints; i++){
|
---|
598 | x1jp1[i] = this.x1[i] + this.incrX1;
|
---|
599 | if(x1jp1[i]>this.x1[this.lPoints-1])x1jp1[i]=this.x1[this.lPoints-1];
|
---|
600 | x1jm1[i] = this.x1[i] - this.incrX1;
|
---|
601 | if(x1jm1[i]<this.x1[0])x1jm1[i]=this.x1[0];
|
---|
602 | }
|
---|
603 | for(int i=0; i<this.mPoints; i++){
|
---|
604 | x2jp1[i] = this.x2[i] + this.incrX2;
|
---|
605 | if(x2jp1[i]>this.x2[this.mPoints-1])x2jp1[i]=this.x2[this.mPoints-1];
|
---|
606 | x2jm1[i] = this.x2[i] - this.incrX2;
|
---|
607 | if(x2jm1[i]<this.x2[0])x2jm1[i]=this.x2[0];
|
---|
608 | }
|
---|
609 | for(int i=0; i<this.nPoints; i++){
|
---|
610 | x3jp1[i] = this.x3[i] + this.incrX3;
|
---|
611 | if(x3jp1[i]>this.x3[this.nPoints-1])x3jp1[i]=this.x3[this.nPoints-1];
|
---|
612 | x3jm1[i] = this.x3[i] - this.incrX3;
|
---|
613 | if(x3jm1[i]<this.x3[0])x3jm1[i]=this.x3[0];
|
---|
614 | }
|
---|
615 |
|
---|
616 | for(int i=0; i<this.lPoints; i++){
|
---|
617 | for(int j=0; j<this.mPoints; j++){
|
---|
618 | for(int k=0; k<this.nPoints; k++){
|
---|
619 | this.dydx1[i][j][k] = (tcs.interpolate(x1jp1[i],x2[j],x3[k]) - tcs.interpolate(x1jm1[i],x2[j],x3[k]))/(x1jp1[i] - x1jm1[i]);
|
---|
620 | this.dydx2[i][j][k] = (tcs.interpolate(x1[i],x2jp1[j],x3[k]) - tcs.interpolate(x1[i],x2jm1[j],x3[k]))/(x2jp1[j] - x2jm1[j]);
|
---|
621 | this.dydx3[i][j][k] = (tcs.interpolate(x1[i],x2[j],x3jp1[k]) - tcs.interpolate(x1[i],x2[j],x3jm1[k]))/(x3jp1[k] - x3jm1[k]);
|
---|
622 | this.d2ydx1dx2[i][j][k] = (tcs.interpolate(x1jp1[i],x2jp1[j],x3[k]) - tcs.interpolate(x1jp1[i],x2jm1[j],x3[k]) - tcs.interpolate(x1jm1[i],x2jp1[j],x3[k]) + tcs.interpolate(x1jm1[i],x2jm1[j],x3[k]))/((x1jp1[i] - x1jm1[i])*(x2jp1[j] - x2jm1[j]));
|
---|
623 | this.d2ydx1dx3[i][j][k] = (tcs.interpolate(x1jp1[i],x2[j],x3jp1[k]) - tcs.interpolate(x1jp1[i],x2[j],x3jm1[k]) - tcs.interpolate(x1jm1[i],x2[j],x3jp1[k]) + tcs.interpolate(x1jm1[i],x2[j],x3jm1[k]))/((x1jp1[i] - x1jm1[i])*(x3jp1[k] - x3jm1[k]));
|
---|
624 | this.d2ydx2dx3[i][j][k] = (tcs.interpolate(x1[i],x2jp1[j],x3jp1[k]) - tcs.interpolate(x1[i],x2jp1[j],x3jm1[k]) - tcs.interpolate(x1[i],x2jm1[j],x3jp1[k]) + tcs.interpolate(x1[i],x2jm1[j],x3jm1[k]))/((x2jp1[j] - x2jm1[j])*(x3jp1[k] - x3jm1[k]));
|
---|
625 | this.d3ydx1dx2dx3[i][j][k] = ((tcs.interpolate(x1jp1[i],x2jp1[j],x3jp1[k]) - tcs.interpolate(x1jp1[i],x2jm1[j],x3jp1[k]) - tcs.interpolate(x1jm1[i],x2jp1[j],x3jp1[k]) + tcs.interpolate(x1jm1[i],x2jm1[j],x3jp1[k])) - (tcs.interpolate(x1jp1[i],x2jp1[j],x3jm1[k]) - tcs.interpolate(x1jp1[i],x2jm1[j],x3jm1[k]) - tcs.interpolate(x1jm1[i],x2jp1[j],x3jm1[k]) + tcs.interpolate(x1jm1[i],x2jm1[j],x3jm1[k])))/((x1jp1[i] - x1jm1[i])*(x2jp1[j] - x2jm1[j])*(x3jp1[k] - x3jm1[k]));
|
---|
626 | }
|
---|
627 | }
|
---|
628 | }
|
---|
629 | }
|
---|
630 | else{
|
---|
631 | // Numerical differentiation using only provided data points
|
---|
632 | int iip = 0;
|
---|
633 | int iim = 0;
|
---|
634 | int jjp = 0;
|
---|
635 | int jjm = 0;
|
---|
636 | int kkp = 0;
|
---|
637 | int kkm = 0;
|
---|
638 | for(int i=0; i<this.lPoints; i++){
|
---|
639 | iip = i+1;
|
---|
640 | if(iip>=this.lPoints)iip = this.lPoints-1;
|
---|
641 | iim = i-1;
|
---|
642 | if(iim<0)iim = 0;
|
---|
643 | for(int j=0; j<this.mPoints; j++){
|
---|
644 | jjp = j+1;
|
---|
645 | if(jjp>=this.mPoints)jjp = this.mPoints-1;
|
---|
646 | jjm = j-1;
|
---|
647 | if(jjm<0)jjm = 0;
|
---|
648 | for(int k=0; k<this.nPoints; k++){
|
---|
649 | kkp = k+1;
|
---|
650 | if(kkp>=this.nPoints)kkp = this.nPoints-1;
|
---|
651 | kkm = k-1;
|
---|
652 | if(kkm<0)kkm = 0;
|
---|
653 | this.dydx1[i][j][k] = (this.y[iip][j][k] - this.y[iim][j][k])/(this.x1[iip] - this.x1[iim]);
|
---|
654 | this.dydx2[i][j][k] = (this.y[i][jjp][k] - this.y[i][jjm][k])/(this.x2[jjp] - this.x2[jjm]);
|
---|
655 | this.dydx3[i][j][k] = (this.y[i][j][kkp] - this.y[i][j][kkm])/(this.x3[kkp] - this.x3[kkm]);
|
---|
656 | this.d2ydx1dx2[i][j][k] = (this.y[iip][jjp][k] - this.y[iip][jjm][k] - this.y[iim][jjp][k] + this.y[iim][jjm][k])/((this.x1[iip] - this.x1[iim])*(this.x2[jjp] - this.x2[jjm]));
|
---|
657 | this.d2ydx1dx3[i][j][k] = (this.y[iip][j][kkp] - this.y[iip][j][kkm] - this.y[iim][j][kkp] + this.y[iim][j][kkm])/((this.x1[iip] - this.x1[iim])*(this.x3[kkp] - this.x3[kkm]));
|
---|
658 | this.d2ydx2dx3[i][j][k] = (this.y[i][jjp][kkp] - this.y[i][jjp][kkm] - this.y[i][jjm][kkp] + this.y[i][jjm][kkm])/((this.x2[jjp] - this.x2[jjm])*(this.x3[kkp] - this.x3[kkm]));
|
---|
659 | this.d2ydx1dx2[i][j][k] = (this.y[iip][jjp][kkp] - this.y[iip][jjm][kkp] - this.y[iim][jjp][kkp] + this.y[iim][jjm][kkp] - this.y[iip][jjp][kkm] + this.y[iip][jjm][kkm] + this.y[iim][jjp][kkm] - this.y[iim][jjm][kkm])/((this.x1[iip] - this.x1[iim])*(this.x2[jjp] - this.x2[jjm])*(this.x3[kkp] - this.x3[kkm]));
|
---|
660 | }
|
---|
661 | }
|
---|
662 | }
|
---|
663 | }
|
---|
664 |
|
---|
665 | this.derivCalculated = true;
|
---|
666 | }
|
---|
667 |
|
---|
668 | // Grid coefficients
|
---|
669 | private void gridCoefficients(){
|
---|
670 |
|
---|
671 | double[] yt = new double[8];
|
---|
672 | double[] dydx1t = new double[8];
|
---|
673 | double[] dydx2t = new double[8];
|
---|
674 | double[] dydx3t = new double[8];
|
---|
675 | double[] d2ydx1dx2t = new double[8];
|
---|
676 | double[] d2ydx1dx3t = new double[8];
|
---|
677 | double[] d2ydx2dx3t = new double[8];
|
---|
678 | double[] d3ydx1dx2dx3t = new double[8];
|
---|
679 | double[] ct = new double[64];
|
---|
680 | double[] xt = new double[64];
|
---|
681 | double d1 = 0.0;
|
---|
682 | double d2 = 0.0;
|
---|
683 | double d3 = 0.0;
|
---|
684 | for(int i=0; i<this.lPoints-1; i++){
|
---|
685 | d1 = this.x1[i+1] - this.x1[i];
|
---|
686 | for(int j=0; j<this.mPoints-1; j++){
|
---|
687 | d2 = this.x2[j+1] - this.x2[j];
|
---|
688 | for(int k=0; k<this.nPoints-1; k++){
|
---|
689 | d3 = this.x3[k+1] - this.x3[k];
|
---|
690 | double[][][] cc = new double[4][4][4];
|
---|
691 | coeff.add(new Double(d1));
|
---|
692 | coeff.add(new Double(this.x1[i]));
|
---|
693 | coeff.add(new Double(d2));
|
---|
694 | coeff.add(new Double(this.x2[j]));
|
---|
695 | coeff.add(new Double(d3));
|
---|
696 | coeff.add(new Double(this.x3[k]));
|
---|
697 |
|
---|
698 | for(int ii=0; ii<8; ii++){
|
---|
699 | yt[ii] = this.y[i+unitCube[ii][0]][j+unitCube[ii][1]][k+unitCube[ii][2]];
|
---|
700 | dydx1t[ii] = this.dydx1[i+unitCube[ii][0]][j+unitCube[ii][1]][k+unitCube[ii][2]];
|
---|
701 | dydx2t[ii] = this.dydx2[i+unitCube[ii][0]][j+unitCube[ii][1]][k+unitCube[ii][2]];
|
---|
702 | dydx3t[ii] = this.dydx3[i+unitCube[ii][0]][j+unitCube[ii][1]][k+unitCube[ii][2]];
|
---|
703 | d2ydx1dx2t[ii] = this.d2ydx1dx2[i+unitCube[ii][0]][j+unitCube[ii][1]][k+unitCube[ii][2]];
|
---|
704 | d2ydx1dx3t[ii] = this.d2ydx1dx3[i+unitCube[ii][0]][j+unitCube[ii][1]][k+unitCube[ii][2]];
|
---|
705 | d2ydx2dx3t[ii] = this.d2ydx2dx3[i+unitCube[ii][0]][j+unitCube[ii][1]][k+unitCube[ii][2]];
|
---|
706 | d3ydx1dx2dx3t[ii] = this.d3ydx1dx2dx3[i+unitCube[ii][0]][j+unitCube[ii][1]][k+unitCube[ii][2]];
|
---|
707 | }
|
---|
708 |
|
---|
709 | for(int k2=0; k2<8; k2++){
|
---|
710 | xt[k2] = yt[k2];
|
---|
711 | xt[k2+8] = dydx1t[k2]*d1;
|
---|
712 | xt[k2+16] = dydx2t[k2]*d2;
|
---|
713 | xt[k2+24] = dydx3t[k2]*d3;
|
---|
714 | xt[k2+32] = d2ydx1dx2t[k2]*d1*d2;
|
---|
715 | xt[k2+40] = d2ydx1dx3t[k2]*d1*d3;
|
---|
716 | xt[k2+48] = d2ydx2dx3t[k2]*d2*d3;
|
---|
717 | xt[k2+56] = d3ydx1dx2dx3t[k2]*d1*d2*d3;
|
---|
718 | }
|
---|
719 |
|
---|
720 | double xh = 0.0;
|
---|
721 | for(int k2=0; k2<64; k2++){
|
---|
722 | for(int kk=0; kk<64; kk++){
|
---|
723 | xh += this.weights[k2][kk]*xt[kk];
|
---|
724 | }
|
---|
725 | ct[k2] = xh;
|
---|
726 | xh = 0.0;
|
---|
727 | }
|
---|
728 | int counter = 0;
|
---|
729 | for(int k2=0; k2<4; k2++){
|
---|
730 | for(int kk=0; kk<4; kk++){
|
---|
731 | for(int kkk=0; kkk<4; kkk++){
|
---|
732 | cc[k2][kk][kkk] = ct[counter++];
|
---|
733 | }
|
---|
734 | }
|
---|
735 | }
|
---|
736 |
|
---|
737 | // Add grid coefficient array to ArrayList
|
---|
738 | coeff.add(cc);
|
---|
739 | }
|
---|
740 | }
|
---|
741 | }
|
---|
742 | }
|
---|
743 |
|
---|
744 | // Returns an interpolated value of y for a value of x
|
---|
745 | // from a tabulated function y=f(x1,x2)
|
---|
746 | public double interpolate(double xx1, double xx2, double xx3){
|
---|
747 | // check that xx1 and xx2 are within the limits
|
---|
748 | if(xx1<x1[0]){
|
---|
749 | if(xx1>=x1[0]-TriCubicInterpolation.potentialRoundingError){
|
---|
750 | xx1=this.x1[0];
|
---|
751 | }
|
---|
752 | else{
|
---|
753 | throw new IllegalArgumentException(xx1 + " is outside the limits, " + x1[0] + " - " + x1[this.lPoints-1]);
|
---|
754 | }
|
---|
755 | }
|
---|
756 | if(xx2<x2[0]){
|
---|
757 | if(xx2>=x2[0]-TriCubicInterpolation.potentialRoundingError){
|
---|
758 | xx2=this.x2[0];
|
---|
759 | }
|
---|
760 | else{
|
---|
761 | throw new IllegalArgumentException(xx2 + " is outside the limits, " + x2[0] + " - " + x2[this.mPoints-1]);
|
---|
762 | }
|
---|
763 | }
|
---|
764 | if(xx3<x3[0]){
|
---|
765 | if(xx3>=x3[0]-TriCubicInterpolation.potentialRoundingError){
|
---|
766 | xx3=this.x3[0];
|
---|
767 | }
|
---|
768 | else{
|
---|
769 | throw new IllegalArgumentException(xx1 + " is outside the limits, " + x3[0] + " - " + x3[this.nPoints-1]);
|
---|
770 | }
|
---|
771 | }
|
---|
772 |
|
---|
773 |
|
---|
774 | if(xx1>this.x1[this.lPoints-1]){
|
---|
775 | if(xx1<=this.x1[this.lPoints-1]+TriCubicInterpolation.potentialRoundingError){
|
---|
776 | xx1=this.x1[this.lPoints-1];
|
---|
777 | }
|
---|
778 | else{
|
---|
779 | throw new IllegalArgumentException(xx1 + " is outside the limits, " + this.x1[0] + " - " + this.x1[this.lPoints-1]);
|
---|
780 | }
|
---|
781 | }
|
---|
782 | if(xx2>this.x2[this.mPoints-1]){
|
---|
783 | if(xx2<=this.x2[this.mPoints-1]+TriCubicInterpolation.potentialRoundingError){
|
---|
784 | xx2=this.x2[this.mPoints-1];
|
---|
785 | }
|
---|
786 | else{
|
---|
787 | throw new IllegalArgumentException(xx2 + " is outside the limits, " + this.x2[0] + " - " + this.x2[this.mPoints-1]);
|
---|
788 | }
|
---|
789 | }
|
---|
790 | if(xx3>this.x3[this.nPoints-1]){
|
---|
791 | if(xx3<=this.x3[this.nPoints-1]+TriCubicInterpolation.potentialRoundingError){
|
---|
792 | xx3=this.x3[this.mPoints-1];
|
---|
793 | }
|
---|
794 | else{
|
---|
795 | throw new IllegalArgumentException(xx3 + " is outside the limits, " + this.x3[0] + " - " + this.x3[this.nPoints-1]);
|
---|
796 | }
|
---|
797 | }
|
---|
798 |
|
---|
799 |
|
---|
800 | // assign variables
|
---|
801 | this.xx1 = xx1;
|
---|
802 | this.xx2 = xx2;
|
---|
803 | this.xx3 = xx3;
|
---|
804 |
|
---|
805 | // Find grid surrounding the interpolation point
|
---|
806 | int gridn =0;
|
---|
807 | double distance1 = ((Double)coeff.get(7*gridn)).doubleValue();
|
---|
808 | double x1lower = ((Double)coeff.get(7*gridn+1)).doubleValue();
|
---|
809 | double distance2 = ((Double)coeff.get(7*gridn+2)).doubleValue();
|
---|
810 | double x2lower = ((Double)coeff.get(7*gridn+3)).doubleValue();
|
---|
811 | double distance3 = ((Double)coeff.get(7*gridn+4)).doubleValue();
|
---|
812 | double x3lower = ((Double)coeff.get(7*gridn+5)).doubleValue();
|
---|
813 | boolean test = true;
|
---|
814 | while(test){
|
---|
815 | boolean test1 = false;
|
---|
816 | boolean test2 = false;
|
---|
817 | boolean test3 = false;
|
---|
818 | if(xx1>=x1lower && xx1<=(x1lower+distance1))test1=true;
|
---|
819 | if(xx2>=x2lower && xx2<=(x2lower+distance2))test2=true;
|
---|
820 | if(xx3>=x3lower && xx3<=(x3lower+distance3))test3=true;
|
---|
821 | if(test1 && test2 && test3){
|
---|
822 | test = false;
|
---|
823 | }
|
---|
824 | else{
|
---|
825 | gridn++;
|
---|
826 | distance1 = ((Double)coeff.get(7*gridn)).doubleValue();
|
---|
827 | x1lower = ((Double)coeff.get(7*gridn+1)).doubleValue();
|
---|
828 | distance2 = ((Double)coeff.get(7*gridn+2)).doubleValue();
|
---|
829 | x2lower = ((Double)coeff.get(7*gridn+3)).doubleValue();
|
---|
830 | distance3 = ((Double)coeff.get(7*gridn+4)).doubleValue();
|
---|
831 | x3lower = ((Double)coeff.get(7*gridn+5)).doubleValue();
|
---|
832 | }
|
---|
833 | }
|
---|
834 | double[][][] gCoeff = (double[][][])coeff.get(7*gridn+6);
|
---|
835 | double x1Normalised = (xx1 - x1lower)/distance1;
|
---|
836 | double x2Normalised = (xx2 - x2lower)/distance2;
|
---|
837 | double x3Normalised = (xx3 - x3lower)/distance3;
|
---|
838 |
|
---|
839 | // interpolation
|
---|
840 | this.interpolatedValue = 0.0; // interpolated value of y
|
---|
841 | for(int i=0; i<4; i++){
|
---|
842 | for(int j=0; j<4; j++){
|
---|
843 | for(int k=0; k<4; k++){
|
---|
844 | this.interpolatedValue += gCoeff[i][j][k]*Math.pow(x1Normalised, i)*Math.pow(x2Normalised, j)*Math.pow(x3Normalised, k);
|
---|
845 | }
|
---|
846 | }
|
---|
847 | }
|
---|
848 | this.interpolatedDydx1 = 0.0; // interpolated value of dy/dx1
|
---|
849 | for(int i=1; i<4; i++){
|
---|
850 | for(int j=0; j<4; j++){
|
---|
851 | for(int k=0; k<4; k++){
|
---|
852 | this.interpolatedDydx1 += i*gCoeff[i][j][k]*Math.pow(x1Normalised, i-1)*Math.pow(x2Normalised, j)*Math.pow(x3Normalised, k);
|
---|
853 | }
|
---|
854 | }
|
---|
855 | }
|
---|
856 | this.interpolatedDydx2 = 0.0; // interpolated value of dydx2
|
---|
857 | for(int i=0; i<4; i++){
|
---|
858 | for(int j=1; j<4; j++){
|
---|
859 | for(int k=0; k<4; k++){
|
---|
860 | this.interpolatedDydx2 += j*gCoeff[i][j][k]*Math.pow(x1Normalised, i)*Math.pow(x2Normalised, j-1)*Math.pow(x3Normalised, k);
|
---|
861 | }
|
---|
862 | }
|
---|
863 | }
|
---|
864 | this.interpolatedDydx3 = 0.0; // interpolated value of dydx3
|
---|
865 | for(int i=0; i<4; i++){
|
---|
866 | for(int j=1; j<4; j++){
|
---|
867 | for(int k=0; k<4; k++){
|
---|
868 | this.interpolatedDydx2 += k*gCoeff[i][j][k]*Math.pow(x1Normalised, i)*Math.pow(x2Normalised, j)*Math.pow(x3Normalised, k-1);
|
---|
869 | }
|
---|
870 | }
|
---|
871 | }
|
---|
872 | this.interpolatedD2ydx1dx2 = 0.0; // interpolated value of d2y/dx1dx2
|
---|
873 | for(int i=1; i<4; i++){
|
---|
874 | for(int j=1; j<4; j++){
|
---|
875 | for(int k=0; k<4; k++){
|
---|
876 | this.interpolatedD2ydx1dx2 += i*j*gCoeff[i][j][k]*Math.pow(x1Normalised, i-1)*Math.pow(x2Normalised, j-1)*Math.pow(x3Normalised, k);
|
---|
877 | }
|
---|
878 | }
|
---|
879 | }
|
---|
880 |
|
---|
881 | return this.interpolatedValue;
|
---|
882 | }
|
---|
883 |
|
---|
884 | // Return last interpolated value and the interpolated gradients
|
---|
885 | public double[] getInterpolatedValues(){
|
---|
886 | double[] ret = new double[11];
|
---|
887 | ret[0] = this.interpolatedValue;
|
---|
888 | ret[1] = this.interpolatedDydx1;
|
---|
889 | ret[2] = this.interpolatedDydx2;
|
---|
890 | ret[3] = this.interpolatedDydx3;
|
---|
891 | ret[4] = this.interpolatedD2ydx1dx2;
|
---|
892 | ret[5] = this.interpolatedD2ydx1dx3;
|
---|
893 | ret[6] = this.interpolatedD2ydx2dx3;
|
---|
894 | ret[7] = this.interpolatedD3ydx1dx2dx3;
|
---|
895 | ret[8] = this.xx1;
|
---|
896 | ret[9] = this.xx2;
|
---|
897 | ret[10] = this.xx3;
|
---|
898 | return ret;
|
---|
899 | }
|
---|
900 |
|
---|
901 | // Return grid point values of dydx1
|
---|
902 | public double[][][] getGridDydx1(){
|
---|
903 | double[][][] ret = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
904 | for(int i=0; i<this.lPoints; i++){
|
---|
905 | for(int j=0; j<this.mPoints; j++){
|
---|
906 | for(int k=0; k<this.nPoints; k++){
|
---|
907 | ret[this.x1indices[i]][this.x2indices[j]][this.x3indices[k]] = this.dydx1[i][j][k];
|
---|
908 | }
|
---|
909 | }
|
---|
910 | }
|
---|
911 | return ret;
|
---|
912 | }
|
---|
913 |
|
---|
914 | // Return grid point values of dydx2
|
---|
915 | public double[][][] getGridDydx2(){
|
---|
916 | double[][][] ret = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
917 | for(int i=0; i<this.lPoints; i++){
|
---|
918 | for(int j=0; j<this.mPoints; j++){
|
---|
919 | for(int k=0; k<this.nPoints; k++){
|
---|
920 | ret[this.x1indices[i]][this.x2indices[j]][this.x3indices[k]] = this.dydx2[i][j][k];
|
---|
921 | }
|
---|
922 | }
|
---|
923 | }
|
---|
924 | return ret;
|
---|
925 | }
|
---|
926 |
|
---|
927 | // Return grid point values of dydx3
|
---|
928 | public double[][][] getGridDydx3(){
|
---|
929 | double[][][] ret = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
930 | for(int i=0; i<this.lPoints; i++){
|
---|
931 | for(int j=0; j<this.mPoints; j++){
|
---|
932 | for(int k=0; k<this.nPoints; k++){
|
---|
933 | ret[this.x1indices[i]][this.x2indices[j]][this.x3indices[k]] = this.dydx3[i][j][k];
|
---|
934 | }
|
---|
935 | }
|
---|
936 | }
|
---|
937 | return ret;
|
---|
938 | }
|
---|
939 |
|
---|
940 | // Return grid point values of d2ydx1dx2
|
---|
941 | public double[][][] getGridD2ydx1dx2(){
|
---|
942 | double[][][] ret = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
943 | for(int i=0; i<this.lPoints; i++){
|
---|
944 | for(int j=0; j<this.mPoints; j++){
|
---|
945 | for(int k=0; k<this.nPoints; k++){
|
---|
946 | ret[this.x1indices[i]][this.x2indices[j]][this.x3indices[k]] = this.d2ydx1dx2[i][j][k];
|
---|
947 | }
|
---|
948 | }
|
---|
949 | }
|
---|
950 | return ret;
|
---|
951 | }
|
---|
952 |
|
---|
953 | // Return grid point values of d2ydx1dx3
|
---|
954 | public double[][][] getGridD2ydx1dx3(){
|
---|
955 | double[][][] ret = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
956 | for(int i=0; i<this.lPoints; i++){
|
---|
957 | for(int j=0; j<this.mPoints; j++){
|
---|
958 | for(int k=0; k<this.nPoints; k++){
|
---|
959 | ret[this.x1indices[i]][this.x2indices[j]][this.x3indices[k]] = this.d2ydx1dx3[i][j][k];
|
---|
960 | }
|
---|
961 | }
|
---|
962 | }
|
---|
963 | return ret;
|
---|
964 | }
|
---|
965 |
|
---|
966 | // Return grid point values of d2ydx2dx3
|
---|
967 | public double[][][] getGridD2ydx2dx3(){
|
---|
968 | double[][][] ret = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
969 | for(int i=0; i<this.lPoints; i++){
|
---|
970 | for(int j=0; j<this.mPoints; j++){
|
---|
971 | for(int k=0; k<this.nPoints; k++){
|
---|
972 | ret[this.x1indices[i]][this.x2indices[j]][this.x3indices[k]] = this.d2ydx2dx3[i][j][k];
|
---|
973 | }
|
---|
974 | }
|
---|
975 | }
|
---|
976 | return ret;
|
---|
977 | }
|
---|
978 |
|
---|
979 | // Return grid point values of d3ydx1dx2dx3
|
---|
980 | public double[][][] getGridD3ydx1dx2dx3(){
|
---|
981 | double[][][] ret = new double[this.lPoints][this.mPoints][this.nPoints];
|
---|
982 | for(int i=0; i<this.lPoints; i++){
|
---|
983 | for(int j=0; j<this.mPoints; j++){
|
---|
984 | for(int k=0; k<this.nPoints; k++){
|
---|
985 | ret[this.x1indices[i]][this.x2indices[j]][this.x3indices[k]] = this.d3ydx1dx2dx3[i][j][k];
|
---|
986 | }
|
---|
987 | }
|
---|
988 | }
|
---|
989 | return ret;
|
---|
990 | }
|
---|
991 |
|
---|
992 | // Reset the numerical differentiation incremental factor delta
|
---|
993 | public static void resetDelta(double delta){
|
---|
994 | TriCubicInterpolation.delta = delta;
|
---|
995 | }
|
---|
996 |
|
---|
997 | // Reset rounding error check option
|
---|
998 | // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
|
---|
999 | // This method causes this check to be ignored and an exception to be thrown if any point lies outside the interpolation bounds
|
---|
1000 | public static void noRoundingErrorCheck(){
|
---|
1001 | TriCubicInterpolation.roundingCheck = false;
|
---|
1002 | TriCubicInterpolation.potentialRoundingError = 0.0;
|
---|
1003 | }
|
---|
1004 |
|
---|
1005 | // Reset potential rounding error value
|
---|
1006 | // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
|
---|
1007 | // The default value for the potential rounding error is 5e-15*times the 10^exponent of the value outside the bounds
|
---|
1008 | // This method allows the 5e-15 to be reset
|
---|
1009 | public static void potentialRoundingError(double potentialRoundingError){
|
---|
1010 | TriCubicInterpolation.potentialRoundingError = potentialRoundingError;
|
---|
1011 | }
|
---|
1012 |
|
---|
1013 | // Get minimum limits
|
---|
1014 | public double[] getXmin(){
|
---|
1015 | return this.xMin;
|
---|
1016 | }
|
---|
1017 |
|
---|
1018 | // Get maximum limits
|
---|
1019 | public double[] getXmax(){
|
---|
1020 | return this.xMax;
|
---|
1021 | }
|
---|
1022 |
|
---|
1023 | // Get limits to x
|
---|
1024 | public double[] getLimits(){
|
---|
1025 | double[] limits = {xMin[0], xMax[0], xMin[1], xMax[1], xMin[2], xMax[2]};
|
---|
1026 | return limits;
|
---|
1027 | }
|
---|
1028 |
|
---|
1029 | // Display limits to x
|
---|
1030 | public void displayLimits(){
|
---|
1031 | System.out.println(" ");
|
---|
1032 | for(int i=0; i<3; i++){
|
---|
1033 | System.out.println("The limits to the x array x" + (i+1) + " are " + xMin[i] + " and " + xMax[i]);
|
---|
1034 | }
|
---|
1035 | System.out.println(" ");
|
---|
1036 | }
|
---|
1037 |
|
---|
1038 | }
|
---|
1039 |
|
---|