1 | /**********************************************************
|
---|
2 | *
|
---|
3 | * Class CubicSplineShort (Deprecated)
|
---|
4 | *
|
---|
5 | * See Class CubicSplineFast for replacement class
|
---|
6 | * This class has been retained purely for compatibility purposes
|
---|
7 | * and will not be updated
|
---|
8 | *
|
---|
9 | * Class for performing an interpolation using a cubic spline
|
---|
10 | * setTabulatedArrays and interpolate adapted, with modification to
|
---|
11 | * an object-oriented approach, from Numerical Recipes in C (http://www.nr.com/)
|
---|
12 | * Stripped down version of CubicSpline - all data checks have been removed for faster running
|
---|
13 | *
|
---|
14 | *
|
---|
15 | * WRITTEN BY: Dr Michael Thomas Flanagan
|
---|
16 | *
|
---|
17 | * DATE: 26 December 2009 (Stripped down version of CubicSpline: May 2002 - 31 October 2009)
|
---|
18 | * DEPRECATED: 31 Decemebr 2009
|
---|
19 | *
|
---|
20 | * DOCUMENTATION:
|
---|
21 | * See Michael Thomas Flanagan's Java library on-line web page:
|
---|
22 | * http://www.ee.ucl.ac.uk/~mflanaga/java/CubicSplineShort.html
|
---|
23 | * http://www.ee.ucl.ac.uk/~mflanaga/java/
|
---|
24 | *
|
---|
25 | * Copyright (c) 2002 - 2009 Michael Thomas Flanagan
|
---|
26 | *
|
---|
27 | * PERMISSION TO COPY:
|
---|
28 | *
|
---|
29 | * Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee,
|
---|
30 | * provided that an acknowledgement to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies
|
---|
31 | * and associated documentation or publications.
|
---|
32 | *
|
---|
33 | * Redistributions of the source code of this source code, or parts of the source codes, must retain the above copyright notice,
|
---|
34 | * this list of conditions and the following disclaimer and requires written permission from the Michael Thomas Flanagan:
|
---|
35 | *
|
---|
36 | * Redistribution in binary form of all or parts of this class must reproduce the above copyright notice, this list of conditions and
|
---|
37 | * the following disclaimer in the documentation and/or other materials provided with the distribution and requires written permission
|
---|
38 | * from the Michael Thomas Flanagan:
|
---|
39 | *
|
---|
40 | * Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose.
|
---|
41 | * Dr Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software
|
---|
42 | * or its derivatives.
|
---|
43 | *
|
---|
44 | ***************************************************************************************/
|
---|
45 |
|
---|
46 |
|
---|
47 | package agents.anac.y2015.agentBuyogV2.flanagan.interpolation;
|
---|
48 |
|
---|
49 | public class CubicSplineShort{
|
---|
50 |
|
---|
51 | private int nPoints = 0; // no. of tabulated points
|
---|
52 | private double[] y = null; // y=f(x) tabulated function
|
---|
53 | private double[] x = null; // x in tabulated function f(x)
|
---|
54 | private double[] d2ydx2 = null; // second derivatives of y
|
---|
55 | private boolean derivCalculated = false; // = true when the derivatives have been calculated
|
---|
56 |
|
---|
57 | // Constructors
|
---|
58 | // Constructor with data arrays initialised to arrays x and y
|
---|
59 | public CubicSplineShort(double[] x, double[] y){
|
---|
60 | this.nPoints=x.length;
|
---|
61 | this.x = new double[nPoints];
|
---|
62 | this.y = new double[nPoints];
|
---|
63 | this.d2ydx2 = new double[nPoints];
|
---|
64 | for(int i=0; i<this.nPoints; i++){
|
---|
65 | this.x[i]=x[i];
|
---|
66 | this.y[i]=y[i];
|
---|
67 | }
|
---|
68 | this.calcDeriv();
|
---|
69 | }
|
---|
70 |
|
---|
71 | // Constructor with data arrays initialised to zero
|
---|
72 | // Primarily for use by BiCubicSplineShort
|
---|
73 | public CubicSplineShort(int nPoints){
|
---|
74 | this.nPoints=nPoints;
|
---|
75 | this.x = new double[nPoints];
|
---|
76 | this.y = new double[nPoints];
|
---|
77 | this.d2ydx2 = new double[nPoints];
|
---|
78 | }
|
---|
79 |
|
---|
80 | // METHODS
|
---|
81 |
|
---|
82 | // Resets the x y data arrays - primarily for use in BiCubicSplineShort
|
---|
83 | public void resetData(double[] x, double[] y){
|
---|
84 | for(int i=0; i<this.nPoints; i++){
|
---|
85 | this.x[i]=x[i];
|
---|
86 | this.y[i]=y[i];
|
---|
87 | }
|
---|
88 | }
|
---|
89 |
|
---|
90 | // Returns a new CubicSplineShort setting array lengths to n and all array values to zero with natural spline default
|
---|
91 | // Primarily for use in BiCubicSplineShort
|
---|
92 | public static CubicSplineShort zero(int n){
|
---|
93 | if(n<3)throw new IllegalArgumentException("A minimum of three data points is needed");
|
---|
94 | CubicSplineShort aa = new CubicSplineShort(n);
|
---|
95 | return aa;
|
---|
96 | }
|
---|
97 |
|
---|
98 | // Create a one dimensional array of cubic spline objects of length n each of array length m
|
---|
99 | // Primarily for use in BiCubicSplineShort
|
---|
100 | public static CubicSplineShort[] oneDarray(int n, int m){
|
---|
101 | CubicSplineShort[] a =new CubicSplineShort[n];
|
---|
102 | for(int i=0; i<n; i++){
|
---|
103 | a[i]=CubicSplineShort.zero(m);
|
---|
104 | }
|
---|
105 | return a;
|
---|
106 | }
|
---|
107 |
|
---|
108 |
|
---|
109 | // Calculates the second derivatives of the tabulated function
|
---|
110 | // for use by the cubic spline interpolation method (.interpolate)
|
---|
111 | // This method follows the procedure in Numerical Methods C language procedure for calculating second derivatives
|
---|
112 | public void calcDeriv(){
|
---|
113 | double p=0.0D,qn=0.0D,sig=0.0D,un=0.0D;
|
---|
114 | double[] u = new double[nPoints];
|
---|
115 |
|
---|
116 | d2ydx2[0]=u[0]=0.0;
|
---|
117 | for(int i=1;i<=this.nPoints-2;i++){
|
---|
118 | sig=(this.x[i]-this.x[i-1])/(this.x[i+1]-this.x[i-1]);
|
---|
119 | p=sig*this.d2ydx2[i-1]+2.0;
|
---|
120 | this.d2ydx2[i]=(sig-1.0)/p;
|
---|
121 | 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]);
|
---|
122 | u[i]=(6.0*u[i]/(this.x[i+1]-this.x[i-1])-sig*u[i-1])/p;
|
---|
123 | }
|
---|
124 |
|
---|
125 | qn=un=0.0;
|
---|
126 | this.d2ydx2[this.nPoints-1]=(un-qn*u[this.nPoints-2])/(qn*this.d2ydx2[this.nPoints-2]+1.0);
|
---|
127 | for(int k=this.nPoints-2;k>=0;k--){
|
---|
128 | this.d2ydx2[k]=this.d2ydx2[k]*this.d2ydx2[k+1]+u[k];
|
---|
129 | }
|
---|
130 | this.derivCalculated = true;
|
---|
131 | }
|
---|
132 |
|
---|
133 | // INTERPOLATE
|
---|
134 | // Returns an interpolated value of y for a value of x from a tabulated function y=f(x)
|
---|
135 | // after the data has been entered via a constructor.
|
---|
136 | // The derivatives are calculated, bt calcDeriv(), on the first call to this method ands are
|
---|
137 | // then stored for use on all subsequent calls
|
---|
138 | public double interpolate(double xx){
|
---|
139 |
|
---|
140 | double h=0.0D,b=0.0D,a=0.0D, yy=0.0D;
|
---|
141 | int k=0;
|
---|
142 | int klo=0;
|
---|
143 | int khi=this.nPoints-1;
|
---|
144 | while (khi-klo > 1){
|
---|
145 | k=(khi+klo) >> 1;
|
---|
146 | if(this.x[k] > xx){
|
---|
147 | khi=k;
|
---|
148 | }
|
---|
149 | else{
|
---|
150 | klo=k;
|
---|
151 | }
|
---|
152 | }
|
---|
153 | h=this.x[khi]-this.x[klo];
|
---|
154 |
|
---|
155 | if (h == 0.0){
|
---|
156 | throw new IllegalArgumentException("Two values of x are identical: point "+klo+ " ("+this.x[klo]+") and point "+khi+ " ("+this.x[khi]+")" );
|
---|
157 | }
|
---|
158 | else{
|
---|
159 | a=(this.x[khi]-xx)/h;
|
---|
160 | b=(xx-this.x[klo])/h;
|
---|
161 | 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;
|
---|
162 | }
|
---|
163 | return yy;
|
---|
164 | }
|
---|
165 | } |
---|