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 |
|
---|
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 |
|
---|
46 | public 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 |
|
---|