source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/interpolation/CubicSplineShort.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: 6.7 KB
Line 
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
47package agents.anac.y2015.agentBuyogV2.flanagan.interpolation;
48
49public 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}
Note: See TracBrowser for help on using the repository browser.