source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/interpolation/TriCubicInterpolation.java

Last change on this file was 127, checked in by Wouter Pasman, 6 years ago

#41 ROLL BACK of rev.126 . So this version is equal to rev. 125

File size: 49.1 KB
Line 
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
38package agents.anac.y2015.agentBuyogV2.flanagan.interpolation;
39
40import java.util.ArrayList;
41
42import agents.anac.y2015.agentBuyogV2.flanagan.math.ArrayMaths;
43import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
44import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
45import agents.anac.y2015.agentBuyogV2.flanagan.math.Matrix;
46
47public 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
Note: See TracBrowser for help on using the repository browser.