1 | /**********************************************************
|
---|
2 | *
|
---|
3 | * Class CubicSpline
|
---|
4 | *
|
---|
5 | * Class for performing an interpolation using a cubic spline
|
---|
6 | * setTabulatedArrays and interpolate adapted, with modification to
|
---|
7 | * an object-oriented approach, from Numerical Recipes in C (http://www.nr.com/)
|
---|
8 | *
|
---|
9 | *
|
---|
10 | * WRITTEN BY: Dr Michael Thomas Flanagan
|
---|
11 | *
|
---|
12 | * DATE: May 2002
|
---|
13 | * UPDATE: 29 April 2005, 17 February 2006, 21 September 2006, 4 December 2007
|
---|
14 | * 24 March 2008 (Thanks to Peter Neuhaus, Florida Institute for Human and Machine Cognition)
|
---|
15 | * 21 September 2008
|
---|
16 | * 14 January 2009 - point deletion and check for 3 points reordered (Thanks to Jan Sacha, Vrije Universiteit Amsterdam)
|
---|
17 | * 31 October 2009, 11-14 January 2010
|
---|
18 | * 31 August 2010 - overriding natural spline error, introduced in earlier update, corrected (Thanks to Dagmara Oszkiewicz, University of Helsinki)
|
---|
19 | * 2 February 2011
|
---|
20 | *
|
---|
21 | * DOCUMENTATION:
|
---|
22 | * See Michael Thomas Flanagan's Java library on-line web page:
|
---|
23 | * http://www.ee.ucl.ac.uk/~mflanaga/java/CubicSpline.html
|
---|
24 | * http://www.ee.ucl.ac.uk/~mflanaga/java/
|
---|
25 | *
|
---|
26 | * Copyright (c) 2002 - 2010 Michael Thomas Flanagan
|
---|
27 | *
|
---|
28 | * PERMISSION TO COPY:
|
---|
29 | *
|
---|
30 | * Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee,
|
---|
31 | * provided that an acknowledgement to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies
|
---|
32 | * and associated documentation or publications.
|
---|
33 | *
|
---|
34 | * Redistributions of the source code of this source code, or parts of the source codes, must retain the above copyright notice,
|
---|
35 | * this list of conditions and the following disclaimer and requires written permission from the Michael Thomas Flanagan:
|
---|
36 | *
|
---|
37 | * Redistribution in binary form of all or parts of this class must reproduce the above copyright notice, this list of conditions and
|
---|
38 | * the following disclaimer in the documentation and/or other materials provided with the distribution and requires written permission
|
---|
39 | * from the Michael Thomas Flanagan:
|
---|
40 | *
|
---|
41 | * Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose.
|
---|
42 | * Dr Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software
|
---|
43 | * or its derivatives.
|
---|
44 | *
|
---|
45 | ***************************************************************************************/
|
---|
46 |
|
---|
47 |
|
---|
48 | package agents.anac.y2015.agentBuyogV2.flanagan.interpolation;
|
---|
49 |
|
---|
50 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
|
---|
51 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
|
---|
52 |
|
---|
53 | public class CubicSpline{
|
---|
54 |
|
---|
55 | private int nPoints = 0; // no. of tabulated points
|
---|
56 | private int nPointsOriginal = 0; // no. of tabulated points after any deletions of identical points
|
---|
57 | private double[] y = null; // y=f(x) tabulated function
|
---|
58 | private double[] x = null; // x in tabulated function f(x)
|
---|
59 | private double yy = Double.NaN; // interpolated value of y
|
---|
60 | private double dydx = Double.NaN; // interpolated value of the first derivative, dy/dx
|
---|
61 | private int[]newAndOldIndices; // record of indices on ordering x into ascending order
|
---|
62 | private double xMin = Double.NaN; // minimum x value
|
---|
63 | private double xMax = Double.NaN; // maximum x value
|
---|
64 | private double range = Double.NaN; // xMax - xMin
|
---|
65 | private double[] d2ydx2 = null; // second derivatives of y
|
---|
66 | private double yp1 = Double.NaN; // second derivative at point one
|
---|
67 | // default value = NaN (natural spline)
|
---|
68 | private double ypn = Double.NaN; // second derivative at point n
|
---|
69 | // default value = NaN (natural spline)
|
---|
70 | private boolean derivCalculated = false; // = true when the derivatives have been calculated
|
---|
71 |
|
---|
72 | private boolean checkPoints = false; // = true when points checked for identical values
|
---|
73 | private static boolean supress = false; // if true: warning messages supressed
|
---|
74 |
|
---|
75 | private static boolean averageIdenticalAbscissae = false; // if true: the the ordinate values for identical abscissae are averaged
|
---|
76 | // if false: the abscissae values are separated by 0.001 of the total abscissae range;
|
---|
77 | private static double potentialRoundingError = 5e-15; // potential rounding error used in checking wheter a value lies within the interpolation bounds
|
---|
78 | private static boolean roundingCheck = true; // = true: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
|
---|
79 |
|
---|
80 |
|
---|
81 | // Constructors
|
---|
82 | // Constructor with data arrays initialised to arrays x and y
|
---|
83 | public CubicSpline(double[] x, double[] y){
|
---|
84 | this.nPoints=x.length;
|
---|
85 | this.nPointsOriginal = this.nPoints;
|
---|
86 | if(this.nPoints!=y.length)throw new IllegalArgumentException("Arrays x and y are of different length "+ this.nPoints + " " + y.length);
|
---|
87 | if(this.nPoints<3)throw new IllegalArgumentException("A minimum of three data points is needed");
|
---|
88 | this.x = new double[nPoints];
|
---|
89 | this.y = new double[nPoints];
|
---|
90 | this.d2ydx2 = new double[nPoints];
|
---|
91 | for(int i=0; i<this.nPoints; i++){
|
---|
92 | this.x[i]=x[i];
|
---|
93 | this.y[i]=y[i];
|
---|
94 | }
|
---|
95 | this.orderPoints();
|
---|
96 | this.checkForIdenticalPoints();
|
---|
97 | this.calcDeriv();
|
---|
98 |
|
---|
99 | }
|
---|
100 |
|
---|
101 | // Constructor with data arrays initialised to zero
|
---|
102 | // Primarily for use by BiCubicSpline
|
---|
103 | public CubicSpline(int nPoints){
|
---|
104 | this.nPoints=nPoints;
|
---|
105 | this.nPointsOriginal = this.nPoints;
|
---|
106 | if(this.nPoints<3)throw new IllegalArgumentException("A minimum of three data points is needed");
|
---|
107 | this.x = new double[nPoints];
|
---|
108 | this.y = new double[nPoints];
|
---|
109 | this.d2ydx2 = new double[nPoints];
|
---|
110 | }
|
---|
111 |
|
---|
112 | // METHODS
|
---|
113 | // Reset rounding error check option
|
---|
114 | // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
|
---|
115 | // This method causes this check to be ignored and an exception to be thrown if any poit lies outside the interpolation bounds
|
---|
116 | public static void noRoundingErrorCheck(){
|
---|
117 | CubicSpline.roundingCheck = false;
|
---|
118 | }
|
---|
119 |
|
---|
120 | // Reset potential rounding error value
|
---|
121 | // Default option: points outside the interpolation bounds by less than the potential rounding error rounded to the bounds limit
|
---|
122 | // The default value for the potential rounding error is 5e-15*times the 10^exponent of the value outside the bounds
|
---|
123 | // This method allows the 5e-15 to be reset
|
---|
124 | public static void potentialRoundingError(double potentialRoundingError){
|
---|
125 | CubicSpline.potentialRoundingError = potentialRoundingError;
|
---|
126 | }
|
---|
127 |
|
---|
128 | // Resets the x y data arrays - primarily for use in BiCubicSpline
|
---|
129 | public void resetData(double[] x, double[] y){
|
---|
130 | this.nPoints = this.nPointsOriginal;
|
---|
131 | if(x.length!=y.length)throw new IllegalArgumentException("Arrays x and y are of different length");
|
---|
132 | if(this.nPoints!=x.length)throw new IllegalArgumentException("Original array length not matched by new array length");
|
---|
133 |
|
---|
134 | for(int i=0; i<this.nPoints; i++){
|
---|
135 | this.x[i]=x[i];
|
---|
136 | this.y[i]=y[i];
|
---|
137 | }
|
---|
138 | this.orderPoints();
|
---|
139 | this.checkForIdenticalPoints();
|
---|
140 | }
|
---|
141 |
|
---|
142 | // Reset the default handing of identical abscissae with different ordinates
|
---|
143 | // from the default option of separating the two relevant abscissae by 0.001 of the range
|
---|
144 | // to averaging the relevant ordinates
|
---|
145 | public static void averageIdenticalAbscissae(){
|
---|
146 | CubicSpline.averageIdenticalAbscissae = true;
|
---|
147 | }
|
---|
148 |
|
---|
149 | // Sort points into an ascending abscissa order
|
---|
150 | public void orderPoints(){
|
---|
151 | double[] dummy = new double[nPoints];
|
---|
152 | this.newAndOldIndices = new int[nPoints];
|
---|
153 | // Sort x into ascending order storing indices changes
|
---|
154 | Fmath.selectionSort(this.x, dummy, this.newAndOldIndices);
|
---|
155 | // Sort x into ascending order and make y match the new order storing both new x and new y
|
---|
156 | Fmath.selectionSort(this.x, this.y, this.x, this.y);
|
---|
157 |
|
---|
158 | // Minimum and maximum values and range
|
---|
159 | this.xMin = Fmath.minimum(this.x);
|
---|
160 | this.xMax = Fmath.maximum(this.x);
|
---|
161 | range = xMax - xMin;
|
---|
162 | }
|
---|
163 |
|
---|
164 | // get the maximum value
|
---|
165 | public double getXmax(){
|
---|
166 | return this.xMax;
|
---|
167 | }
|
---|
168 |
|
---|
169 | // get the minimum value
|
---|
170 | public double getXmin(){
|
---|
171 | return this.xMin;
|
---|
172 | }
|
---|
173 |
|
---|
174 | // get the limits of x
|
---|
175 | public double[] getLimits(){
|
---|
176 | double[] limits = {this.xMin, this.xMax};
|
---|
177 | return limits;
|
---|
178 | }
|
---|
179 |
|
---|
180 | // print to screen the limis of x
|
---|
181 | public void displayLimits(){
|
---|
182 | System.out.println("\nThe limits of the abscissae (x-values) are " + this.xMin + " and " + this.xMax +"\n");
|
---|
183 | }
|
---|
184 |
|
---|
185 |
|
---|
186 | // Checks for and removes all but one of identical points
|
---|
187 | // Checks and appropriately handles identical abscissae with differing ordinates
|
---|
188 | public void checkForIdenticalPoints(){
|
---|
189 | int nP = this.nPoints;
|
---|
190 | boolean test1 = true;
|
---|
191 | int ii = 0;
|
---|
192 | while(test1){
|
---|
193 | boolean test2 = true;
|
---|
194 | int jj = ii+1;
|
---|
195 | while(test2){
|
---|
196 | if(this.x[ii]==this.x[jj]){
|
---|
197 | if(this.y[ii]==this.y[jj]){
|
---|
198 | if(!CubicSpline.supress){
|
---|
199 | System.out.print("CubicSpline: Two identical points, " + this.x[ii] + ", " + this.y[ii]);
|
---|
200 | System.out.println(", in data array at indices " + this.newAndOldIndices[ii] + " and " + this.newAndOldIndices[jj] + ", latter point removed");
|
---|
201 | }
|
---|
202 | double[] xx = new double[this.nPoints-1];
|
---|
203 | double[] yy = new double[this.nPoints-1];
|
---|
204 | int[] naoi = new int[this.nPoints-1];
|
---|
205 | for(int i=0; i<jj; i++){
|
---|
206 | xx[i] = this.x[i];
|
---|
207 | yy[i] = this.y[i];
|
---|
208 | naoi[i] = this.newAndOldIndices[i];
|
---|
209 | }
|
---|
210 | for(int i=jj; i<this.nPoints-1; i++){
|
---|
211 | xx[i] = this.x[i+1];
|
---|
212 | yy[i] = this.y[i+1];
|
---|
213 | naoi[i] = this.newAndOldIndices[i+1];
|
---|
214 | }
|
---|
215 | this.nPoints--;
|
---|
216 | this.x = Conv.copy(xx);
|
---|
217 | this.y = Conv.copy(yy);
|
---|
218 | this.newAndOldIndices = Conv.copy(naoi);
|
---|
219 | }
|
---|
220 | else{
|
---|
221 | if(CubicSpline.averageIdenticalAbscissae==true){
|
---|
222 | if(!CubicSpline.supress){
|
---|
223 | System.out.print("CubicSpline: Two identical points on the absicca (x-axis) with different ordinate (y-axis) values, " + x[ii] + ": " + y[ii] + ", " + y[jj]);
|
---|
224 | System.out.println(", average of the ordinates taken");
|
---|
225 | }
|
---|
226 | this.y[ii] = (this.y[ii] + this.y[jj])/2.0D;
|
---|
227 | double[] xx = new double[this.nPoints-1];
|
---|
228 | double[] yy = new double[this.nPoints-1];
|
---|
229 | int[] naoi = new int[this.nPoints-1];
|
---|
230 | for(int i=0; i<jj; i++){
|
---|
231 | xx[i] = this.x[i];
|
---|
232 | yy[i] = this.y[i];
|
---|
233 | naoi[i] = this.newAndOldIndices[i];
|
---|
234 | }
|
---|
235 | for(int i=jj; i<this.nPoints-1; i++){
|
---|
236 | xx[i] = this.x[i+1];
|
---|
237 | yy[i] = this.y[i+1];
|
---|
238 | naoi[i] = this.newAndOldIndices[i+1];
|
---|
239 | }
|
---|
240 | this.nPoints--;
|
---|
241 | this.x = Conv.copy(xx);
|
---|
242 | this.y = Conv.copy(yy);
|
---|
243 | this.newAndOldIndices = Conv.copy(naoi);
|
---|
244 | }
|
---|
245 | else{
|
---|
246 | double sepn = range*0.0005D;
|
---|
247 | if(!CubicSpline.supress)System.out.print("CubicSpline: Two identical points on the absicca (x-axis) with different ordinate (y-axis) values, " + x[ii] + ": " + y[ii] + ", " + y[jj]);
|
---|
248 | boolean check = false;
|
---|
249 | if(ii==0){
|
---|
250 | if(x[2]-x[1]<=sepn)sepn = (x[2]-x[1])/2.0D;
|
---|
251 | if(this.y[0]>this.y[1]){
|
---|
252 | if(this.y[1]>this.y[2]){
|
---|
253 | check = stay(ii, jj, sepn);
|
---|
254 | }
|
---|
255 | else{
|
---|
256 | check = swap(ii, jj, sepn);
|
---|
257 | }
|
---|
258 | }
|
---|
259 | else{
|
---|
260 | if(this.y[2]<=this.y[1]){
|
---|
261 | check = swap(ii, jj, sepn);
|
---|
262 | }
|
---|
263 | else{
|
---|
264 | check = stay(ii, jj, sepn);
|
---|
265 | }
|
---|
266 | }
|
---|
267 | }
|
---|
268 | if(jj==this.nPoints-1){
|
---|
269 | if(x[nP-2]-x[nP-3]<=sepn)sepn = (x[nP-2]-x[nP-3])/2.0D;
|
---|
270 | if(this.y[ii]<=this.y[jj]){
|
---|
271 | if(this.y[ii-1]<=this.y[ii]){
|
---|
272 | check = stay(ii, jj, sepn);
|
---|
273 | }
|
---|
274 | else{
|
---|
275 | check = swap(ii, jj, sepn);
|
---|
276 | }
|
---|
277 | }
|
---|
278 | else{
|
---|
279 | if(this.y[ii-1]<=this.y[ii]){
|
---|
280 | check = swap(ii, jj, sepn);
|
---|
281 | }
|
---|
282 | else{
|
---|
283 | check = stay(ii, jj, sepn);
|
---|
284 | }
|
---|
285 | }
|
---|
286 | }
|
---|
287 | if(ii!=0 && jj!=this.nPoints-1){
|
---|
288 | if(x[ii]-x[ii-1]<=sepn)sepn = (x[ii]-x[ii-1])/2;
|
---|
289 | if(x[jj+1]-x[jj]<=sepn)sepn = (x[jj+1]-x[jj])/2;
|
---|
290 | if(this.y[ii]>this.y[ii-1]){
|
---|
291 | if(this.y[jj]>this.y[ii]){
|
---|
292 | if(this.y[jj]>this.y[jj+1]){
|
---|
293 | if(this.y[ii-1]<=this.y[jj+1]){
|
---|
294 | check = stay(ii, jj, sepn);
|
---|
295 | }
|
---|
296 | else{
|
---|
297 | check = swap(ii, jj, sepn);
|
---|
298 | }
|
---|
299 | }
|
---|
300 | else{
|
---|
301 | check = stay(ii, jj, sepn);
|
---|
302 | }
|
---|
303 | }
|
---|
304 | else{
|
---|
305 | if(this.y[jj+1]>this.y[jj]){
|
---|
306 | if(this.y[jj+1]>this.y[ii-1] && this.y[jj+1]>this.y[ii-1]){
|
---|
307 | check = stay(ii, jj, sepn);
|
---|
308 | }
|
---|
309 | }
|
---|
310 | else{
|
---|
311 | check = swap(ii, jj, sepn);
|
---|
312 | }
|
---|
313 | }
|
---|
314 | }
|
---|
315 | else{
|
---|
316 | if(this.y[jj]>this.y[ii]){
|
---|
317 | if(this.y[jj+1]>this.y[jj]){
|
---|
318 | check = stay(ii, jj, sepn);
|
---|
319 | }
|
---|
320 | }
|
---|
321 | else{
|
---|
322 | if(this.y[jj+1]>this.y[ii-1]){
|
---|
323 | check = stay(ii, jj, sepn);
|
---|
324 | }
|
---|
325 | else{
|
---|
326 | check = swap(ii, jj, sepn);
|
---|
327 | }
|
---|
328 | }
|
---|
329 | }
|
---|
330 | }
|
---|
331 |
|
---|
332 | if(check==false){
|
---|
333 | check = stay(ii, jj, sepn);
|
---|
334 | }
|
---|
335 | if(!CubicSpline.supress)System.out.println(", the two abscissae have been separated by a distance " + sepn);
|
---|
336 | jj++;
|
---|
337 | }
|
---|
338 | }
|
---|
339 | if((this.nPoints-1)==ii)test2 = false;
|
---|
340 | }
|
---|
341 | else{
|
---|
342 | jj++;
|
---|
343 | }
|
---|
344 | if(jj>=this.nPoints)test2 = false;
|
---|
345 | }
|
---|
346 | ii++;
|
---|
347 | if(ii>=this.nPoints-1)test1 = false;
|
---|
348 | }
|
---|
349 | if(this.nPoints<3)throw new IllegalArgumentException("Removal of duplicate points has reduced the number of points to less than the required minimum of three data points");
|
---|
350 |
|
---|
351 | this.checkPoints = true;
|
---|
352 | }
|
---|
353 |
|
---|
354 | // Swap method for checkForIdenticalPoints procedure
|
---|
355 | private boolean swap(int ii, int jj, double sepn){
|
---|
356 | this.x[ii] += sepn;
|
---|
357 | this.x[jj] -= sepn;
|
---|
358 | double hold = this.x[ii];
|
---|
359 | this.x[ii] = this.x[jj];
|
---|
360 | this.x[jj] = hold;
|
---|
361 | hold = this.y[ii];
|
---|
362 | this.y[ii] = this.y[jj];
|
---|
363 | this.y[jj] = hold;
|
---|
364 | return true;
|
---|
365 | }
|
---|
366 |
|
---|
367 | // Stay method for checkForIdenticalPoints procedure
|
---|
368 | private boolean stay(int ii, int jj, double sepn){
|
---|
369 | this.x[ii] -= sepn;
|
---|
370 | this.x[jj] += sepn;
|
---|
371 | return true;
|
---|
372 | }
|
---|
373 |
|
---|
374 | // Supress warning messages in the identifiaction of duplicate points
|
---|
375 | public static void supress(){
|
---|
376 | CubicSpline.supress = true;
|
---|
377 | }
|
---|
378 |
|
---|
379 | // Unsupress warning messages in the identifiaction of duplicate points
|
---|
380 | public static void unsupress(){
|
---|
381 | CubicSpline.supress = false;
|
---|
382 | }
|
---|
383 |
|
---|
384 | // Returns a new CubicSpline setting array lengths to n and all array values to zero with natural spline default
|
---|
385 | // Primarily for use in BiCubicSpline
|
---|
386 | public static CubicSpline zero(int n){
|
---|
387 | if(n<3)throw new IllegalArgumentException("A minimum of three data points is needed");
|
---|
388 | CubicSpline aa = new CubicSpline(n);
|
---|
389 | return aa;
|
---|
390 | }
|
---|
391 |
|
---|
392 | // Create a one dimensional array of cubic spline objects of length n each of array length m
|
---|
393 | // Primarily for use in BiCubicSpline
|
---|
394 | public static CubicSpline[] oneDarray(int n, int m){
|
---|
395 | if(m<3)throw new IllegalArgumentException("A minimum of three data points is needed");
|
---|
396 | CubicSpline[] a =new CubicSpline[n];
|
---|
397 | for(int i=0; i<n; i++){
|
---|
398 | a[i]=CubicSpline.zero(m);
|
---|
399 | }
|
---|
400 | return a;
|
---|
401 | }
|
---|
402 |
|
---|
403 | // Enters the first derivatives of the cubic spline at
|
---|
404 | // the first and last point of the tabulated data
|
---|
405 | // Overrides a natural spline
|
---|
406 | public void setDerivLimits(double yp1, double ypn){
|
---|
407 | this.yp1=yp1;
|
---|
408 | this.ypn=ypn;
|
---|
409 | this.calcDeriv();
|
---|
410 | }
|
---|
411 |
|
---|
412 | // Resets a natural spline
|
---|
413 | // Use above - this kept for backward compatibility
|
---|
414 | public void setDerivLimits(){
|
---|
415 | this.yp1=Double.NaN;
|
---|
416 | this.ypn=Double.NaN;
|
---|
417 | this.calcDeriv();
|
---|
418 | }
|
---|
419 |
|
---|
420 | // Enters the first derivatives of the cubic spline at
|
---|
421 | // the first and last point of the tabulated data
|
---|
422 | // Overrides a natural spline
|
---|
423 | // Use setDerivLimits(double yp1, double ypn) - this kept for backward compatibility
|
---|
424 | public void setDeriv(double yp1, double ypn){
|
---|
425 | this.yp1=yp1;
|
---|
426 | this.ypn=ypn;
|
---|
427 | this.calcDeriv();
|
---|
428 | }
|
---|
429 |
|
---|
430 | // Returns the internal array of second derivatives
|
---|
431 | public double[] getDeriv(){
|
---|
432 | if(!this.derivCalculated)this.calcDeriv();
|
---|
433 | return this.d2ydx2;
|
---|
434 | }
|
---|
435 |
|
---|
436 | // Sets the internal array of second derivatives
|
---|
437 | // Used primarily with BiCubicSpline
|
---|
438 | public void setDeriv(double[] deriv){
|
---|
439 | this.d2ydx2 = deriv;
|
---|
440 | this.derivCalculated = true;
|
---|
441 | }
|
---|
442 |
|
---|
443 | // Calculates the second derivatives of the tabulated function
|
---|
444 | // for use by the cubic spline interpolation method (.interpolate)
|
---|
445 | // This method follows the procedure in Numerical Methods C language procedure for calculating second derivatives
|
---|
446 | public void calcDeriv(){
|
---|
447 | double p=0.0D,qn=0.0D,sig=0.0D,un=0.0D;
|
---|
448 | double[] u = new double[nPoints];
|
---|
449 |
|
---|
450 | if (Double.isNaN(this.yp1)){
|
---|
451 | d2ydx2[0]=u[0]=0.0;
|
---|
452 | }
|
---|
453 | else{
|
---|
454 | this.d2ydx2[0] = -0.5;
|
---|
455 | u[0]=(3.0/(this.x[1]-this.x[0]))*((this.y[1]-this.y[0])/(this.x[1]-this.x[0])-this.yp1);
|
---|
456 | }
|
---|
457 |
|
---|
458 | for(int i=1;i<=this.nPoints-2;i++){
|
---|
459 | sig=(this.x[i]-this.x[i-1])/(this.x[i+1]-this.x[i-1]);
|
---|
460 | p=sig*this.d2ydx2[i-1]+2.0;
|
---|
461 | this.d2ydx2[i]=(sig-1.0)/p;
|
---|
462 | u[i]=(this.y[i+1]-this.y[i])/(this.x[i+1]-this.x[i]) - (this.y[i]-this.y[i-1])/(this.x[i]-this.x[i-1]);
|
---|
463 | u[i]=(6.0*u[i]/(this.x[i+1]-this.x[i-1])-sig*u[i-1])/p;
|
---|
464 | }
|
---|
465 |
|
---|
466 | if (Double.isNaN(this.ypn)){
|
---|
467 | qn=un=0.0;
|
---|
468 | }
|
---|
469 | else{
|
---|
470 | qn=0.5;
|
---|
471 | un=(3.0/(this.x[nPoints-1]-this.x[this.nPoints-2]))*(this.ypn-(this.y[this.nPoints-1]-this.y[this.nPoints-2])/(this.x[this.nPoints-1]-x[this.nPoints-2]));
|
---|
472 | }
|
---|
473 |
|
---|
474 | this.d2ydx2[this.nPoints-1]=(un-qn*u[this.nPoints-2])/(qn*this.d2ydx2[this.nPoints-2]+1.0);
|
---|
475 | for(int k=this.nPoints-2;k>=0;k--){
|
---|
476 | this.d2ydx2[k]=this.d2ydx2[k]*this.d2ydx2[k+1]+u[k];
|
---|
477 | }
|
---|
478 | this.derivCalculated = true;
|
---|
479 | }
|
---|
480 |
|
---|
481 | // INTERPOLATE
|
---|
482 | // Returns an interpolated value of y for a value of x from a tabulated function y=f(x)
|
---|
483 | // after the data has been entered via a constructor.
|
---|
484 | // The derivatives are calculated, bt calcDeriv(), on the first call to this method ands are
|
---|
485 | // then stored for use on all subsequent calls
|
---|
486 | public double interpolate(double xx){
|
---|
487 |
|
---|
488 | // Check for violation of interpolation bounds
|
---|
489 | if (xx<this.x[0]){
|
---|
490 | // if violation is less than potntial rounding error - amend to lie with bounds
|
---|
491 | if(CubicSpline.roundingCheck && Math.abs(x[0]-xx)<=Math.pow(10, Math.floor(Math.log10(Math.abs(this.x[0]))))*CubicSpline.potentialRoundingError){
|
---|
492 | xx = x[0];
|
---|
493 | }
|
---|
494 | else{
|
---|
495 | throw new IllegalArgumentException("x ("+xx+") is outside the range of data points ("+x[0]+" to "+x[this.nPoints-1] + ")");
|
---|
496 | }
|
---|
497 | }
|
---|
498 | if (xx>this.x[this.nPoints-1]){
|
---|
499 | if(CubicSpline.roundingCheck && Math.abs(xx-this.x[this.nPoints-1])<=Math.pow(10, Math.floor(Math.log10(Math.abs(this.x[this.nPoints-1]))))*CubicSpline.potentialRoundingError){
|
---|
500 | xx = this.x[this.nPoints-1];
|
---|
501 | }
|
---|
502 | else{
|
---|
503 | throw new IllegalArgumentException("x ("+xx+") is outside the range of data points ("+x[0]+" to "+x[this.nPoints-1] + ")");
|
---|
504 | }
|
---|
505 | }
|
---|
506 |
|
---|
507 | double h=0.0D,b=0.0D,a=0.0D, yy=0.0D;
|
---|
508 | int k=0;
|
---|
509 | int klo=0;
|
---|
510 | int khi=this.nPoints-1;
|
---|
511 | while (khi-klo > 1){
|
---|
512 | k=(khi+klo) >> 1;
|
---|
513 | if(this.x[k] > xx){
|
---|
514 | khi=k;
|
---|
515 | }
|
---|
516 | else{
|
---|
517 | klo=k;
|
---|
518 | }
|
---|
519 | }
|
---|
520 | h=this.x[khi]-this.x[klo];
|
---|
521 |
|
---|
522 | if (h == 0.0){
|
---|
523 | throw new IllegalArgumentException("Two values of x are identical: point "+klo+ " ("+this.x[klo]+") and point "+khi+ " ("+this.x[khi]+")" );
|
---|
524 | }
|
---|
525 | else{
|
---|
526 | a=(this.x[khi]-xx)/h;
|
---|
527 | b=(xx-this.x[klo])/h;
|
---|
528 | this.yy = a*this.y[klo]+b*this.y[khi]+((a*a*a-a)*this.d2ydx2[klo]+(b*b*b-b)*this.d2ydx2[khi])*(h*h)/6.0;
|
---|
529 | this.dydx = (y[khi] - y[klo])/h - ((3*a*a - 1.0)*this.d2ydx2[klo] - (3*b*b - 1.0)*this.d2ydx2[khi])*h/6.0;
|
---|
530 | }
|
---|
531 | return this.yy;
|
---|
532 | }
|
---|
533 |
|
---|
534 |
|
---|
535 | // Returns an interpolated value of y and of the first derivative dy/dx for a value of x from a tabulated function y=f(x)
|
---|
536 | // after the data has been entered via a constructor.
|
---|
537 | public double[] interpolate_for_y_and_dydx(double xx){
|
---|
538 | this.interpolate(xx);
|
---|
539 | double[] ret = {this.yy, this.dydx};
|
---|
540 | return ret;
|
---|
541 | }
|
---|
542 |
|
---|
543 |
|
---|
544 | // Returns an interpolated value of y for a value of x (xx) from a tabulated function y=f(x)
|
---|
545 | // after the derivatives (deriv) have been calculated independently of calcDeriv().
|
---|
546 | public static double interpolate(double xx, double[] x, double[] y, double[] deriv){
|
---|
547 |
|
---|
548 | if(((x.length != y.length) || (x.length != deriv.length)) || (y.length != deriv.length)){
|
---|
549 | throw new IllegalArgumentException("array lengths are not all equal");
|
---|
550 | }
|
---|
551 | int n = x.length;
|
---|
552 | double h=0.0D, b=0.0D, a=0.0D, yy = 0.0D;
|
---|
553 |
|
---|
554 | int k=0;
|
---|
555 | int klo=0;
|
---|
556 | int khi=n-1;
|
---|
557 | while(khi-klo > 1){
|
---|
558 | k=(khi+klo) >> 1;
|
---|
559 | if(x[k] > xx){
|
---|
560 | khi=k;
|
---|
561 | }
|
---|
562 | else{
|
---|
563 | klo=k;
|
---|
564 | }
|
---|
565 | }
|
---|
566 | h=x[khi]-x[klo];
|
---|
567 |
|
---|
568 | if (h == 0.0){
|
---|
569 | throw new IllegalArgumentException("Two values of x are identical");
|
---|
570 | }
|
---|
571 | else{
|
---|
572 | a=(x[khi]-xx)/h;
|
---|
573 | b=(xx-x[klo])/h;
|
---|
574 | yy=a*y[klo]+b*y[khi]+((a*a*a-a)*deriv[klo]+(b*b*b-b)*deriv[khi])*(h*h)/6.0;
|
---|
575 | }
|
---|
576 | return yy;
|
---|
577 | }
|
---|
578 |
|
---|
579 | }
|
---|