1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
use std::time::Duration;
#[derive(Debug, Clone, Copy)]
pub struct SpringParams {
pub damping: f64,
pub mass: f64,
pub stiffness: f64,
pub epsilon: f64,
}
#[derive(Debug, Clone, Copy)]
pub struct Spring {
pub from: f64,
pub to: f64,
pub initial_velocity: f64,
pub params: SpringParams,
}
impl SpringParams {
pub fn new(damping_ratio: f64, stiffness: f64, epsilon: f64) -> Self {
let damping_ratio = damping_ratio.max(0.);
let stiffness = stiffness.max(0.);
let epsilon = epsilon.max(0.);
let mass = 1.;
let critical_damping = 2. * (mass * stiffness).sqrt();
let damping = damping_ratio * critical_damping;
Self {
damping,
mass,
stiffness,
epsilon,
}
}
}
impl Spring {
pub fn value_at(&self, t: Duration) -> f64 {
self.oscillate(t.as_secs_f64())
}
// Based on libadwaita (LGPL-2.1-or-later):
// https://gitlab.gnome.org/GNOME/libadwaita/-/blob/1.4.4/src/adw-spring-animation.c,
// which itself is based on (MIT):
// https://github.com/robb/RBBAnimation/blob/master/RBBAnimation/RBBSpringAnimation.m
/// Computes and returns the duration until the spring is at rest.
pub fn duration(&self) -> Duration {
const DELTA: f64 = 0.001;
let beta = self.params.damping / (2. * self.params.mass);
if beta.abs() <= f64::EPSILON || beta < 0. {
return Duration::MAX;
}
let omega0 = (self.params.stiffness / self.params.mass).sqrt();
// As first ansatz for the overdamped solution,
// and general estimation for the oscillating ones
// we take the value of the envelope when it's < epsilon.
let mut x0 = -self.params.epsilon.ln() / beta;
// f64::EPSILON is too small for this specific comparison, so we use
// f32::EPSILON even though it's doubles.
if (beta - omega0).abs() <= f64::from(f32::EPSILON) || beta < omega0 {
return Duration::from_secs_f64(x0);
}
// Since the overdamped solution decays way slower than the envelope
// we need to use the value of the oscillation itself.
// Newton's root finding method is a good candidate in this particular case:
// https://en.wikipedia.org/wiki/Newton%27s_method
let mut y0 = self.oscillate(x0);
let m = (self.oscillate(x0 + DELTA) - y0) / DELTA;
let mut x1 = (self.to - y0 + m * x0) / m;
let mut y1 = self.oscillate(x1);
let mut i = 0;
while (self.to - y1).abs() > self.params.epsilon {
if i > 1000 {
return Duration::ZERO;
}
x0 = x1;
y0 = y1;
let m = (self.oscillate(x0 + DELTA) - y0) / DELTA;
x1 = (self.to - y0 + m * x0) / m;
y1 = self.oscillate(x1);
i += 1;
}
Duration::from_secs_f64(x1)
}
/// Returns the spring position at a given time in seconds.
fn oscillate(&self, t: f64) -> f64 {
let b = self.params.damping;
let m = self.params.mass;
let k = self.params.stiffness;
let v0 = self.initial_velocity;
let beta = b / (2. * m);
let omega0 = (k / m).sqrt();
let x0 = self.from - self.to;
let envelope = (-beta * t).exp();
// Solutions of the form C1*e^(lambda1*x) + C2*e^(lambda2*x)
// for the differential equation m*ẍ+b*ẋ+kx = 0
// f64::EPSILON is too small for this specific comparison, so we use
// f32::EPSILON even though it's doubles.
if (beta - omega0).abs() <= f64::from(f32::EPSILON) {
// Critically damped.
self.to + envelope * (x0 + (beta * x0 + v0) * t)
} else if beta < omega0 {
// Underdamped.
let omega1 = ((omega0 * omega0) - (beta * beta)).sqrt();
self.to
+ envelope
* (x0 * (omega1 * t).cos() + ((beta * x0 + v0) / omega1) * (omega1 * t).sin())
} else {
// Overdamped.
let omega2 = ((beta * beta) - (omega0 * omega0)).sqrt();
self.to
+ envelope
* (x0 * (omega2 * t).cosh() + ((beta * x0 + v0) / omega2) * (omega2 * t).sinh())
}
}
}
|
