36 lines
1 KiB
Text
36 lines
1 KiB
Text
const N: usize = 32;
|
|
const STEP: f64 = 0.05;
|
|
|
|
fn main() {
|
|
let x: Vec<f64> = (0..N).map(|i| i as f64 * STEP).collect();
|
|
let sin = x.iter().map(|x| x.sin()).collect::<Vec<_>>();
|
|
let cos = x.iter().map(|x| x.cos()).collect::<Vec<_>>();
|
|
let tan = x.iter().map(|x| x.tan()).collect::<Vec<_>>();
|
|
|
|
println!(
|
|
"{}\n{}\n{}",
|
|
6. * thiele(&sin, &x, 0.5),
|
|
3. * thiele(&cos, &x, 0.5),
|
|
4. * thiele(&tan, &x, 1.)
|
|
);
|
|
}
|
|
|
|
fn thiele(x: &[f64], y: &[f64], xin: f64) -> f64 {
|
|
let mut p: Vec<Vec<f64>> = (0..N).map(|i| (i..N).map(|_| 0.0).collect()).collect();
|
|
|
|
(0..N).for_each(|i| p[i][0] = y[i]);
|
|
|
|
(0..N - 1).for_each(|i| p[i][1] = (x[i] - x[i + 1]) / (p[i][0] - p[i + 1][0]));
|
|
|
|
(2..N).for_each(|i| {
|
|
(0..N - i).for_each(|j| {
|
|
p[j][i] = (x[j] - x[j + i]) / (p[j][i - 1] - p[j + 1][i - 1]) + p[j + 1][i - 2];
|
|
})
|
|
});
|
|
|
|
let mut a = 0.;
|
|
(2..N).rev().for_each(|i| {
|
|
a = (xin - x[i - 1]) / (p[0][i] - p[0][i - 2] + a);
|
|
});
|
|
y[0] + (xin - x[0]) / (p[0][1] + a)
|
|
}
|