rustyqlib/core/
interpolation.rs

1//Algorithm for computing natural cubic splines
2// https://en.wikipedia.org/w/index.php?title=Spline_%28mathematics%29&oldid=288288033#Algorithm_for_computing_natural_cubic_splines
3//https://stackoverflow.com/questions/1204553/are-there-any-good-libraries-for-solving-cubic-splines-in-c
4
5pub struct CubicSpline<'a>{
6    pub x_vec : &'a Vec<f64>,
7    pub y_vec : &'a Vec<f64>,
8    pub spline_set : Vec<SplineSet>
9}
10
11pub struct SplineSet{
12    a:f64,
13    b:f64,
14    c:f64,
15    d:f64,
16    x:f64,
17}
18impl CubicSpline<'_> {
19    pub fn new<'a>(x: &'a Vec<f64>, y: &'a Vec<f64>) -> CubicSpline<'a> {
20        assert!(x.len() == y.len() && x.len() >= 2 && y.len() >= 2, "Must have at least 2 control points.");
21        let n = x.len()-1;
22        let mut a = vec![x[0]];
23        let mut aa = y.clone();
24        a.append(&mut aa);
25        let mut c = vec![0.0; n+1];
26        let mut b = vec![0.0; n];
27        let mut d = vec![0.0; n];
28        let mut dx = Vec::new();
29        for i in 0..n {
30            dx.push(x[i+1]-x[i]);
31        }
32
33        let mut dy = Vec::new();
34        for i in 0..n {
35            dy.push(y[i+1]-y[i]);
36        }
37
38
39        let mut alpha = vec![0.0];
40        for i in 1..n {
41            alpha.push(3.0*(dy[i]/dx[i] - dy[i-1]/dx[i-1]));
42        }
43
44        let mut l = vec![0.0; n+1];
45        let mut mu = vec![0.0; n+1];
46        let mut z = vec![0.0; n+1];
47        l[0] = 1.0;
48        for i in 1..n{
49            l[i] = 2.0*(x[i+1]-x[i-1]) - dx[i-1]*mu[i-1];
50            mu[i] = dx[i]/l[i];
51            z[i] = (alpha[i]-dx[i-1]*z[i-1])/l[i];
52        }
53        l[n] = 1.0;
54        z[n] = 0.0;
55        //c[n] = 0;
56        for j in (0..n).rev(){
57            c[j] = z[j] - mu[j] * c[j+1];
58            b[j] = (a[j+1]-a[j])/dx[j]-dx[j]*(c[j+1]+2.0*c[j])/3.0;
59            d[j] = (c[j+1]-c[j])/(3.0*dx[j]);
60
61        }
62        let mut output = Vec::new();
63        for i in 0..n {
64            let spline_set = SplineSet{
65                a:a[i],
66                b:b[i],
67                c:c[i],
68                d:d[i],
69                x:x[i]
70            };
71            output.push(spline_set);
72        }
73        let spline = CubicSpline{
74            x_vec:x,
75            y_vec:y,
76            spline_set:output
77        };
78        spline
79    }
80    pub fn interpolation(&self,x:f64) -> f64 {
81        let n = self.x_vec.len();
82        for i in 0..n {
83            if x<=self.x_vec[i] {
84                let diff = self.x_vec[i] - x;
85                return self.spline_set[i].a +
86                self.spline_set[i].b * diff +
87                self.spline_set[i].c * diff * diff +
88                self.spline_set[i].d * diff * diff * diff;
89            }
90        }
91    0.0
92    }
93}
94
95
rustyqlib/core/interpolation.rs

rustyqlib/core/
interpolation.rs
