source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/interpolation/BiCubicInterpolation.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: 31.7 KB
Line 
1/**********************************************************
2*
3* BiCubicInterpolation.java
4*
5* Class for performing an interpolation on the tabulated
6* function y = f(x1,x2) using a bicubic interploation procedure
7**
8* WRITTEN BY: Dr Michael Thomas Flanagan
9*
10* DATE: 10-12 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/BiCubicInterpolation.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;
45
46public class BiCubicInterpolation{
47
48 private int nPoints = 0; // no. of x1 tabulated points
49 private int mPoints = 0; // no. of x2 tabulated points
50 private double[] x1 = null; // x1 in tabulated function f(x1,x2)
51 private double[] x2 = null; // x2 in tabulated function f(x1,x2)
52 private double[][] y = null; // y=f(x1,x2) tabulated function
53
54 private double[][] dydx1 = null; // dy/dx1
55 private double[][] dydx2 = null; // dy/dx2
56 private double[][] d2ydx1dx2 = null; // d2y/dx1dx2
57 private boolean derivCalculated = false; // = true when the derivatives have been calculated or entered
58 private BiCubicSpline bcs = null; // bicubic spline used in calculating the derivatives
59 private double incrX1 = 0; // x1 increment used in calculating the derivatives
60 private double incrX2 = 0; // x2 increment used in calculating the derivatives
61
62 private double xx1 = Double.NaN; // value of x1 at which an interpolated y value is required
63 private double xx2 = Double.NaN; // value of x2 at which an interpolated y value is required
64
65 private ArrayList<Object> coeff = new ArrayList<Object>(); // grid square coefficients
66
67 // Weights used in calculating the grid square coefficients
68 private double[][] weights = {{1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
69 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
70 {-3.0,0.0,0.0,3.0,0.0,0.0,0.0,0.0,-2.0,0.0,0.0,-1.0,0.0,0.0,0.0,0.0},
71 {2.0,0.0,0.0,-2.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0},
72 {0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
73 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0},
74 {0.0,0.0,0.0,0.0,-3.0,0.0,0.0,3.0,0.0,0.0,0.0,0.0,-2.0,0.0,0.0,-1.0},
75 {0.0,0.0,0.0,0.0,2.0,0.0,0.0,-2.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,1.0},
76 {-3.0,3.0,0.0,0.0,-2.0,-1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
77 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-3.0,3.0,0.0,0.0,-2.0,-1.0,0.0,0.0},
78 {9.0,-9.0,9.0,-9.0,6.0,3.0,-3.0,-6.0,6.0,-6.0,-3.0,3.0,4.0,2.0,1.0,2.0},
79 {-6.0,6.0,-6.0,6.0,-4.0,-2.0,2.0,4.0,-3.0,3.0,3.0,-3.0,-2.0,-1.0,-1.0,-2.0},
80 {2.0,-2.0,0.0,0.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
81 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,2.0,-2.0,0.0,0.0,1.0,1.0,0.0,0.0},
82 {-6.0,6.0,-6.0,6.0,-3.0,-3.0,3.0,3.0,-4.0,4.0,2.0,-2.0,-2.0,-2.0,-1.0,-1.0},
83 {4.0,-4.0,4.0,-4.0,2.0,2.0,-2.0,-2.0,2.0,-2.0,-2.0,2.0,1.0,1.0,1.0,1.0}};
84
85 private int[] x1indices = null; // x1 data indices before ordering
86 private int[] x2indices = null; // x2 data indices before ordering
87
88 private double[] xMin = new double[2]; // minimum values of x1 and x2
89 private double[] xMax = new double[2]; // maximum values of x1 and x2
90
91 private double interpolatedValue = Double.NaN; // interpolated value of y
92 private double interpolatedDydx1 = Double.NaN; // interpolated value of dydx1
93 private double interpolatedDydx2 = Double.NaN; // interpolated value of dydx2
94 private double interpolatedD2ydx1dx2 = Double.NaN; // interpolated value of d2ydx1dx2
95
96 private boolean numerDiffFlag = true; // = true: if numerical differentiation performed h1 and h2 calculated using delta
97 // = false: if numerical differentiation performed h1 and h2 calculated only provided data points
98
99 private static double delta = 1e-3; // fractional step factor used in calculating the derivatives
100 private static double potentialRoundingError = 5e-15; // potential rounding error used in checking wheter a value lies within the interpolation bounds (static value)
101 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)
102
103 // Constructor without derivatives
104 // numerDiffOption = 0 -> numerical differencing using only supplied data points
105 // numerDiffOption = 1 -> numerical differencing using interpolation
106 public BiCubicInterpolation(double[] x1, double[] x2, double[][] y, int numerDiffOption){
107 // set numerical differencing option
108 if(numerDiffOption==0){
109 this.numerDiffFlag = false;
110 }
111 else{
112 if(numerDiffOption==1){
113 this.numerDiffFlag = true;
114 }
115 else{
116 throw new IllegalArgumentException("The numerical differencing option, " + numerDiffOption + ", must be 0 or 1");
117 }
118 }
119
120 // initialize the data
121 this.initialize(Conv.copy(x1), Conv.copy(x2), Conv.copy(y));
122
123 // calculate the derivatives
124 this.calcDeriv();
125
126 // calculate grid coefficients for all grid squares
127 this.gridCoefficients();
128 }
129
130 // Constructor without derivatives
131 // Numerical differencing by interpolation
132 // Retained for compatability
133 public BiCubicInterpolation(double[] x1, double[] x2, double[][] y){
134 // set numerical differencing option
135 this.numerDiffFlag = true;
136
137 // initialize the data
138 this.initialize(Conv.copy(x1), Conv.copy(x2), Conv.copy(y));
139
140 // calculate the derivatives
141 this.calcDeriv();
142
143 // calculate grid coefficients for all grid squares
144 this.gridCoefficients();
145 }
146
147 // Constructor with derivatives
148 public BiCubicInterpolation(double[] x1, double[] x2, double[][] y, double[][] dydx1, double[][] dydx2, double[][] d2ydx1dx2){
149
150 // initialize the data
151 this.initialize(Conv.copy(x1), Conv.copy(x2), Conv.copy(y), Conv.copy(dydx1), Conv.copy(dydx2), Conv.copy(d2ydx1dx2));
152
153 // calculate grid coefficients for all grid squares
154 this.gridCoefficients();
155 }
156
157 // Initialize the data
158 private void initialize(double[] x1, double[] x2, double[][] y){
159 this.initialize(x1, x2, y, null, null, null, false);
160 }
161
162 private void initialize(double[] x1, double[] x2, double[][] y, double[][] dydx1, double[][] dydx2, double[][] d2ydx1dx2){
163 this.initialize(x1, x2, y, dydx1, dydx2, d2ydx1dx2, true);
164 }
165
166 private void initialize(double[] x1, double[] x2, double[][] y, double[][] dydx1, double[][] dydx2, double[][] d2ydx1dx2, boolean flag){
167 int nPoints=x1.length;
168 int mPoints=x2.length;
169 if(nPoints!=y.length)throw new IllegalArgumentException("Arrays x1 and y-row are of different length " + nPoints + " " + y.length);
170 if(mPoints!=y[0].length)throw new IllegalArgumentException("Arrays x2 and y-column are of different length "+ mPoints + " " + y[0].length);
171 if(nPoints<2 || mPoints<2)throw new IllegalArgumentException("The data matrix must have a minimum size of 2 X 2");
172
173 // order data
174 ArrayMaths am = new ArrayMaths(x1);
175 am = am.sort();
176 this.x1indices = am.originalIndices();
177 x1 = am.array();
178 double[][] hold = new double[nPoints][mPoints];
179 double[][] hold1 = null;
180 double[][] hold2 = null;
181 double[][] hold12 = null;
182
183 for(int i=0; i<nPoints; i++){
184 for(int j=0; j<mPoints; j++){
185 hold[i][j] = y[this.x1indices[i]][j];
186 }
187 }
188 for(int i=0; i<nPoints; i++){
189 for(int j=0; j<mPoints; j++){
190 y[i][j] = hold[i][j];
191 }
192 }
193
194 if(flag){
195 hold1 = new double[nPoints][mPoints];
196 hold2 = new double[nPoints][mPoints];
197 hold12 = new double[nPoints][mPoints];
198 for(int i=0; i<nPoints; i++){
199 for(int j=0; j<mPoints; j++){
200 hold1[i][j] = dydx1[this.x1indices[i]][j];
201 hold2[i][j] = dydx2[this.x1indices[i]][j];
202 hold12[i][j] = d2ydx1dx2[this.x1indices[i]][j];
203 }
204 }
205 for(int i=0; i<nPoints; i++){
206 for(int j=0; j<mPoints; j++){
207 dydx1[i][j] = hold1[i][j];
208 dydx2[i][j] = hold2[i][j];
209 d2ydx1dx2[i][j] = hold12[i][j];
210 }
211 }
212 }
213
214 am = new ArrayMaths(x2);
215 am = am.sort();
216 this.x2indices = am.originalIndices();
217 double[] xh = am.array();
218 for(int i=0; i<mPoints; i++)x2[i] = xh[mPoints-1-i];
219 for(int i=0; i<mPoints; i++){
220 for(int j=0; j<nPoints; j++){
221 hold[j][i] = y[j][this.x2indices[i]];
222 }
223 }
224 for(int i=0; i<nPoints; i++){
225 for(int j=0; j<mPoints; j++){
226 y[i][j] = hold[i][mPoints-1-j];
227 }
228 }
229 if(flag){
230 for(int i=0; i<mPoints; i++){
231 for(int j=0; j<nPoints; j++){
232 hold1[j][i] = dydx1[j][this.x2indices[i]];
233 hold2[j][i] = dydx2[j][this.x2indices[i]];
234 hold12[j][i] = d2ydx1dx2[j][this.x2indices[i]];
235 }
236 }
237 for(int i=0; i<nPoints; i++){
238 for(int j=0; j<mPoints; j++){
239 dydx1[i][j] = hold1[i][mPoints-1-j];
240 dydx2[i][j] = hold2[i][mPoints-1-j];
241 d2ydx1dx2[i][j] = hold12[i][mPoints-1-j];
242 }
243 }
244 }
245
246 // check for identical x1 values
247 for(int i=1; i<nPoints; i++){
248 if(x1[i]==x1[i-1]){
249 System.out.println("x1["+this.x1indices[i]+"] and x1["+this.x1indices[i+1]+"] are identical, " + x1[i]);
250 System.out.println("The y values have been averaged and one point has been deleted");
251 for(int j=i; j<nPoints-1; j++){
252 x1[j]=x1[j+1];
253 this.x1indices[j]=this.x1indices[j+1];
254 }
255 for(int j=0; j<mPoints; j++){
256 y[i-1][j] = (y[i-1][j] + y[i][j])/2.0;
257 for(int k=i; k<nPoints-1; k++)y[k][j]=y[k+1][j];
258 if(flag){
259 dydx1[i-1][j] = (dydx1[i-1][j] + dydx1[i][j])/2.0;
260 dydx2[i-1][j] = (dydx2[i-1][j] + dydx2[i][j])/2.0;
261 d2ydx1dx2[i-1][j] = (d2ydx1dx2[i-1][j] + d2ydx1dx2[i][j])/2.0;
262 for(int k=i; k<nPoints-1; k++){
263 dydx1[k][j]=dydx1[k+1][j];
264 dydx2[k][j]=dydx2[k+1][j];
265 d2ydx1dx2[k][j]=d2ydx1dx2[k+1][j];
266 }
267 }
268 }
269 nPoints--;
270 }
271 }
272 // check for identical x2 values
273 for(int i=1; i<mPoints; i++){
274 if(x2[i]==x2[i-1]){
275 System.out.println("x2["+this.x2indices[i]+"] and x2["+this.x2indices[i]+"] are identical, " + x2[i]);
276 System.out.println("The y values have been averaged and one point has been deleted");
277 for(int j=i; j<mPoints-1; j++){
278 x2[j]=x2[j+1];
279 this.x2indices[j]=this.x2indices[j+1];
280 }
281 for(int j=0; j<nPoints; j++){
282 y[j][i-1] = (y[j][i-1] + y[j][i])/2.0;
283 for(int k=i; k<mPoints-1; k++)y[j][k]=y[j][k+1];
284 if(flag){
285 dydx1[j][i-1] = (dydx1[j][i-1] + dydx1[j][i])/2.0;
286 dydx2[j][i-1] = (dydx2[j][i-1] + dydx2[j][i])/2.0;
287 d2ydx1dx2[j][i-1] = (d2ydx1dx2[j][i-1] + d2ydx1dx2[j][i])/2.0;
288 for(int k=i; k<nPoints-1; k++){
289 dydx1[j][k]=dydx1[j][k+1];
290 dydx2[j][k]=dydx2[j][k+1];
291 d2ydx1dx2[j][k]=d2ydx1dx2[j][k+1];
292 }
293 }
294 }
295 mPoints--;
296 }
297 }
298
299 // assign variables
300 this.nPoints = nPoints;
301 this.mPoints = mPoints;
302 this.x1 = new double[this.nPoints];
303 this.x2 = new double[this.mPoints];
304 this.y = new double[this.nPoints][this.mPoints];
305 this.dydx1 = new double[this.nPoints][this.mPoints];
306 this.dydx2 = new double[this.nPoints][this.mPoints];
307 this.d2ydx1dx2 = new double[this.nPoints][this.mPoints];
308
309 for(int i=0; i<this.nPoints; i++){
310 this.x1[i]=x1[i];
311 }
312 for(int j=0; j<this.mPoints; j++){
313 this.x2[j]=x2[j];
314 }
315 for(int i =0; i<this.nPoints; i++){
316 for(int j=0; j<this.mPoints; j++){
317 this.y[i][j]=y[i][j];
318 }
319 if(flag){
320 for(int j=0; j<this.mPoints; j++){
321 this.dydx1[i][j]=dydx1[i][j];
322 this.dydx2[i][j]=dydx2[i][j];
323 this.d2ydx1dx2[i][j]=d2ydx1dx2[i][j];
324 }
325 }
326 }
327 if(flag)this.derivCalculated = true;
328
329 // limits
330 this.xMin[0] = Fmath.minimum(this.x1);
331 this.xMax[0] = Fmath.maximum(this.x1);
332 this.xMin[1] = Fmath.minimum(this.x2);
333 this.xMax[1] = Fmath.maximum(this.x2);
334
335 if(!flag && this.numerDiffFlag){
336 // numerical difference increments
337 double range1 = this.xMax[0] - this.xMin[0];
338 double range2 = this.xMax[1] - this.xMin[1];
339 double averageSeparation1 = range1/this.nPoints;
340 double averageSeparation2 = range2/this.mPoints;
341 double minSep = this.x1[1] - this.x1[0];
342 double minimumSeparation1 = minSep;
343 for(int i=2; i<this.nPoints; i++){
344 minSep = this.x1[i] - this.x1[i-1];
345 if(minSep<minimumSeparation1)minimumSeparation1 = minSep;
346 }
347 minSep = this.x2[1] - this.x2[0];
348 double minimumSeparation2 = minSep;
349 for(int i=2; i<this.mPoints; i++){
350 minSep = this.x2[i] - this.x2[i-1];
351 if(minSep<minimumSeparation2)minimumSeparation2 = minSep;
352 }
353
354 this.incrX1 = range1*BiCubicInterpolation.delta;
355 double defaultIncr = minimumSeparation1;
356 if(minimumSeparation1<averageSeparation1/10.0)defaultIncr = averageSeparation1/10.0;
357 if(this.incrX1>averageSeparation1)this.incrX1 = defaultIncr;
358 this.incrX2 = range2*BiCubicInterpolation.delta;
359 defaultIncr = minimumSeparation2;
360 if(minimumSeparation2<averageSeparation2/10.0)defaultIncr = averageSeparation2/10.0;
361 if(this.incrX2>averageSeparation2)this.incrX2 = defaultIncr;
362 }
363 }
364
365 // Calculate the derivatives
366 private void calcDeriv(){
367
368 if(this.numerDiffFlag){
369 // Numerical differentiation using delta and interpolation
370 this.bcs = new BiCubicSpline(this.x1, this.x2, this.y);
371 double yjp1k = 0.0;
372 double yjm1k = 0.0;
373 double[] x1jp1 = new double[this.nPoints];
374 double[] x1jm1 = new double[this.nPoints];
375 double[] x2jp1 = new double[this.mPoints];
376 double[] x2jm1 = new double[this.mPoints];
377
378 for(int i=0; i<this.nPoints; i++){
379 x1jp1[i] = this.x1[i] + this.incrX1;
380 if(x1jp1[i]>this.x1[this.nPoints-1])x1jp1[i]=this.x1[this.nPoints-1];
381 x1jm1[i] = this.x1[i] - this.incrX1;
382 if(x1jm1[i]<this.x1[0])x1jm1[i]=this.x1[0];
383 }
384 for(int i=0; i<this.mPoints; i++){
385 x2jp1[i] = this.x2[i] + this.incrX2;
386 if(x2jp1[i]>this.x2[0])x2jp1[i]=this.x2[0];
387 x2jm1[i] = this.x2[i] - this.incrX2;
388 if(x2jm1[i]<this.x2[this.mPoints-1])x2jm1[i]=this.x2[this.mPoints-1];
389
390 }
391 for(int i=0; i<this.nPoints; i++){
392 for(int j=0; j<this.mPoints; j++){
393 this.dydx1[i][j] = (bcs.interpolate(x1jp1[i],x2[j]) - bcs.interpolate(x1jm1[i],x2[j]))/(x1jp1[i] - x1jm1[i]);
394 this.dydx2[i][j] = (bcs.interpolate(x1[i],x2jp1[j]) - bcs.interpolate(x1[i],x2jm1[j]))/(x2jp1[j] - x2jm1[j]);
395 this.d2ydx1dx2[i][j] = (bcs.interpolate(x1jp1[i],x2jp1[j]) - bcs.interpolate(x1jp1[i],x2jm1[j]) - bcs.interpolate(x1jm1[i],x2jp1[j]) + bcs.interpolate(x1jm1[i],x2jm1[j]))/((x1jp1[i] - x1jm1[i])*(x2jp1[j] - x2jm1[j]));
396 }
397 }
398 }
399 else{
400 // Numerical differentiation using only provided data points
401 int iip =0;
402 int iim =0;
403 int jjp =0;
404 int jjm =0;
405 for(int i=0; i<this.nPoints; i++){
406 iip = i+1;
407 if(iip>=this.nPoints)iip = this.nPoints-1;
408 iim = i-1;
409 if(iim<0)iim = 0;
410 for(int j=0; j<this.mPoints; j++){
411 jjp = j+1;
412 if(jjp>=this.mPoints)jjp = this.mPoints-1;
413 jjm = j-1;
414 if(jjm<0)jjm = 0;
415 this.dydx1[i][j] = (this.y[iip][j] - this.y[iim][j])/(this.x1[iip] - this.x1[iim]);
416 this.dydx2[i][j] = (this.y[i][jjp] - this.y[i][jjm])/(this.x2[jjp] - this.x2[jjm]);
417 this.d2ydx1dx2[i][j] = (this.y[iip][jjp] - this.y[iip][jjm] - this.y[iim][jjp] + this.y[iim][jjm])/((this.x1[iip] - this.x1[iim])*(this.x2[jjp] - this.x2[jjm]));
418 }
419 }
420 }
421
422 this.derivCalculated = true;
423 }
424
425 // Grid coefficients
426 private void gridCoefficients(){
427
428 double[] yt = new double[4];
429 double[] dydx1t = new double[4];
430 double[] dydx2t = new double[4];
431 double[] d2ydx1dx2t = new double[4];
432 double[] ct = new double[16];
433 double[] xt = new double[16];
434 double d1 = 0.0;
435 double d2 = 0.0;
436 for(int i=0; i<this.nPoints-1; i++){
437 d1 = this.x1[i+1] - this.x1[i];
438 for(int j=0; j<this.mPoints-1; j++){
439
440 // Calculate grid coefficients for 4-point grid square with point i,j at the top left corner
441 double[][] cc = new double[4][4];
442 d2 = this.x2[j] - this.x2[j+1];
443 coeff.add(new Double(d1));
444 coeff.add(new Double(this.x1[i]));
445 coeff.add(new Double(d2));
446 coeff.add(new Double(this.x2[j+1]));
447 yt[0] = this.y[i][j+1];
448 dydx1t[0] = this.dydx1[i][j+1];
449 dydx2t[0] = this.dydx2[i][j+1];
450 d2ydx1dx2t[0] = this.d2ydx1dx2[i][j+1];
451 yt[1] = this.y[i+1][j+1];
452 dydx1t[1] = this.dydx1[i+1][j+1];
453 dydx2t[1] = this.dydx2[i+1][j+1];
454 d2ydx1dx2t[1] = this.d2ydx1dx2[i+1][j+1];
455 yt[2] = this.y[i+1][j];
456 dydx1t[2] = this.dydx1[i+1][j];
457 dydx2t[2] = this.dydx2[i+1][j];
458 d2ydx1dx2t[2] = this.d2ydx1dx2[i+1][j];
459 yt[3] = this.y[i][j];
460 dydx1t[3] = this.dydx1[i][j];
461 dydx2t[3] = this.dydx2[i][j];
462 d2ydx1dx2t[3] = this.d2ydx1dx2[i][j];
463
464 for(int k=0; k<4; k++){
465 xt[k] = yt[k];
466 xt[k+4] = dydx1t[k]*d1;
467 xt[k+8] = dydx2t[k]*d2;
468 xt[k+12] = d2ydx1dx2t[k]*d1*d2;
469 }
470
471 double xh = 0.0;
472 for(int k=0; k<16; k++){
473 for(int kk=0; kk<16; kk++){
474 xh += this.weights[k][kk]*xt[kk];
475 }
476 ct[k] = xh;
477 xh = 0.0;
478 }
479
480 int counter = 0;
481 for(int k=0; k<4; k++){
482 for(int kk=0; kk<4; kk++){
483 cc[k][kk] = ct[counter++];
484 }
485 }
486
487 // Add grid coefficient 4x4 array to ArrayList
488 coeff.add(cc);
489 }
490 }
491 }
492
493 // Returns an interpolated value of y for a value of x
494 // from a tabulated function y=f(x1,x2)
495 public double interpolate(double xx1, double xx2){
496 // check that xx1 and xx2 are within the limits
497 if(xx1<x1[0]){
498 if(xx1>=x1[0]-BiCubicInterpolation.potentialRoundingError){
499 xx1=this.x1[0];
500 }
501 else{
502 throw new IllegalArgumentException(xx1 + " is outside the limits, " + x1[0] + " - " + x1[this.nPoints-1]);
503 }
504 }
505 if(xx2<this.x2[this.mPoints-1]){
506 if(xx2>this.x2[this.mPoints-1]-BiCubicInterpolation.potentialRoundingError){
507 xx2=this.x2[this.mPoints-1];
508 }
509 else{
510 throw new IllegalArgumentException(xx2 + " is outside the limits, " + this.x2[this.mPoints-1] + " - " + this.x2[0]);
511 }
512 }
513 if(xx1>this.x1[this.nPoints-1]){
514 if(xx1<=this.x1[this.nPoints-1]+BiCubicInterpolation.potentialRoundingError){
515 xx1=this.x1[this.nPoints-1];
516 }
517 else{
518 throw new IllegalArgumentException(xx1 + " is outside the limits, " + this.x1[0] + " - " + this.x1[this.nPoints-1]);
519 }
520 }
521 if(xx2>this.x2[0]){
522 if(xx2<=this.x2[0]+BiCubicInterpolation.potentialRoundingError){
523 xx2=this.x2[0];
524 }
525 else{
526 throw new IllegalArgumentException(xx2 + " is outside the limits, " + this.x2[this.nPoints-1] + " - " + this.x2[0]);
527 }
528 }
529
530 // assign variables
531 this.xx1 = xx1;
532 this.xx2 = xx2;
533
534 // Find grid surrounding the interpolation point
535 int grid1 = 0;
536 int grid2 = 0;
537 int counter = 1;
538 boolean test = true;
539 while(test){
540 if(xx1<this.x1[counter]){
541 grid1 = counter - 1;
542 test = false;
543 }
544 else{
545 counter++;
546 if(counter>=this.nPoints){
547 grid1 = this.nPoints-2;
548 test = false;
549 }
550 }
551 }
552 counter = 0;
553 test = true;
554 while(test){
555 if(xx2>=this.x2[counter+1] && xx2<=this.x2[counter]){
556 grid2 = counter;
557 test = false;
558 }
559 else{
560 counter++;
561 }
562 }
563 int gridn = grid1*(this.mPoints-1) + grid2;
564
565 // grid details
566 double distance1 = ((Double)coeff.get(5*gridn)).doubleValue();
567 double x1lower = ((Double)coeff.get(5*gridn+1)).doubleValue();
568 double distance2 = ((Double)coeff.get(5*gridn+2)).doubleValue();
569 double x2lower = ((Double)coeff.get(5*gridn+3)).doubleValue();
570 double[][] gCoeff = (double[][])coeff.get(5*gridn+4);
571 double x1Normalised = (xx1 - x1lower)/distance1;
572 double x2Normalised = (xx2 - x2lower)/distance2;
573
574 // interpolation
575 this.interpolatedValue = 0.0; // interpolated value of y
576 for(int i=0; i<4; i++){
577 for(int j=0; j<4; j++){
578 this.interpolatedValue += gCoeff[i][j]*Math.pow(x1Normalised, i)*Math.pow(x2Normalised, j);
579 }
580 }
581 this.interpolatedDydx1 = 0.0; // interpolated value of dy/dx1
582 for(int i=1; i<4; i++){
583 for(int j=0; j<4; j++){
584 this.interpolatedDydx1 += i*gCoeff[i][j]*Math.pow(x1Normalised, i-1)*Math.pow(x2Normalised, j);
585 }
586 }
587 this.interpolatedDydx2 = 0.0; // interpolated value of dydx2
588 for(int i=0; i<4; i++){
589 for(int j=1; j<4; j++){
590 this.interpolatedDydx2 += j*gCoeff[i][j]*Math.pow(x1Normalised, i)*Math.pow(x2Normalised, j-1);
591 }
592 }
593 this.interpolatedD2ydx1dx2 = 0.0; // interpolated value of d2y/dx1dx2
594 for(int i=1; i<4; i++){
595 for(int j=1; j<4; j++){
596 this.interpolatedD2ydx1dx2 += i*j*gCoeff[i][j]*Math.pow(x1Normalised, i-1)*Math.pow(x2Normalised, j-1);
597 }
598 }
599
600 return this.interpolatedValue;
601 }
602
603 // Return last interpolated value and the interpolated gradients
604 public double[] getInterpolatedValues(){
605 double[] ret = new double[6];
606 ret[0] = this.interpolatedValue;
607 ret[1] = this.interpolatedDydx1;
608 ret[2] = this.interpolatedDydx2;
609 ret[3] = this.interpolatedD2ydx1dx2;
610 ret[4] = this.xx1;
611 ret[5] = this.xx2;
612 return ret;
613 }
614
615 // Return grid point values of dydx1
616 public double[][] getGridDydx1(){
617 double[][] ret = new double[this.nPoints][this.mPoints];
618 for(int i=0; i<this.nPoints; i++){
619 for(int j=0; j<this.mPoints; j++){
620 ret[this.x1indices[i]][this.x2indices[j]] = this.dydx1[i][j];
621 }
622 }
623 return ret;
624 }
625
626 // Return grid point values of dydx2
627 public double[][] getGridDydx2(){
628 double[][] ret = new double[this.nPoints][this.mPoints];
629 for(int i=0; i<this.nPoints; i++){
630 for(int j=0; j<this.mPoints; j++){
631 ret[this.x1indices[i]][this.x2indices[j]] = this.dydx2[i][j];
632 }
633 }
634 return ret;
635 }
636
637 // Return grid point values of d2ydx1dx2
638 public double[][] getGridD2ydx1dx2(){
639 double[][] ret = new double[this.nPoints][this.mPoints];
640 for(int i=0; i<this.nPoints; i++){
641 for(int j=0; j<this.mPoints; j++){
642 ret[this.x1indices[i]][this.x2indices[j]] = this.d2ydx1dx2[i][j];
643 }
644 }
645 return ret;
646 }
647
648 // Reset the numerical differentiation incremental factor delta
649 public static void resetDelta(double delta){
650 BiCubicInterpolation.delta = delta;
651 }
652
653 // Reset rounding error check option
654 // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
655 // This method causes this check to be ignored and an exception to be thrown if any point lies outside the interpolation bounds
656 public static void noRoundingErrorCheck(){
657 BiCubicInterpolation.roundingCheck = false;
658 BiCubicInterpolation.potentialRoundingError = 0.0;
659 }
660
661 // Reset potential rounding error value
662 // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
663 // The default value for the potential rounding error is 5e-15*times the 10^exponent of the value outside the bounds
664 // This method allows the 5e-15 to be reset
665 public static void potentialRoundingError(double potentialRoundingError){
666 BiCubicInterpolation.potentialRoundingError = potentialRoundingError;
667 }
668
669 // Get minimum limits
670 public double[] getXmin(){
671 return this.xMin;
672 }
673
674 // Get maximum limits
675 public double[] getXmax(){
676 return this.xMax;
677 }
678
679 // Get limits to x
680 public double[] getLimits(){
681 double[] limits = {xMin[0], xMax[0], xMin[1], xMax[1]};
682 return limits;
683 }
684
685 // Display limits to x
686 public void displayLimits(){
687 System.out.println(" ");
688 for(int i=0; i<2; i++){
689 System.out.println("The limits to the x array " + i + " are " + xMin[i] + " and " + xMax[i]);
690 }
691 System.out.println(" ");
692 }
693
694}
695
Note: See TracBrowser for help on using the repository browser.