sva_samples/physics/
chaigne_askenfelt.rs

1// Concern: one hammer/unison-string-set call site's finite-difference state | Non-concern: argument evaluation, builtin dispatch | IO: (params, sr) -> a site; () -> a sample
2
3//! Chaigne & Askenfelt 1994's coupled hammer/stiff-string model. A physics citation, not an
4//! instrument: nothing here, or anywhere this is wired in, may name one.
5
6use crate::error::SampleError;
7use crate::physics::{Solver, Varies};
8
9use crate::physics::ball::Ball;
10use crate::physics::bound::Bound::*;
11use crate::physics::bound::all;
12use crate::physics::hammer::Hammer;
13use crate::physics::stiff_string::{
14    StringGrid, Wire, dispersive_grid, grid_tension, point_weights, read_at, spread, stencil_update,
15};
16use crate::physics::string_energy::{Felt, energy, press};
17use crate::physics::unison_energy::{unison_energy, unison_stable};
18
19#[derive(Clone, Debug, PartialEq)]
20pub struct ChaigneAskenfeltParams {
21    pub f0: f64,
22    pub b: f64,
23    pub strike_pos: f64,
24    pub vel: f64,
25    pub hammer_mass: f64,
26    pub hammer_k: f64,
27    pub hammer_p: f64,
28    pub damp_dc: f64,
29    pub damp_freq: f64,
30    pub unison_count: f64,
31    pub detune: f64,
32    pub bridge_coupling: f64,
33    /// kg; `0` is massless.
34    pub bridge_mass: f64,
35    /// Cents off `detune`'s placing.
36    pub string_cents: [f64; MAX_UNISON],
37    pub string_hammer_k_ratio: [f64; MAX_UNISON],
38    pub damper_pos: f64,
39    /// N s/m.
40    pub damper_r: f64,
41    /// N/m.
42    pub damper_k: f64,
43}
44
45impl ChaigneAskenfeltParams {
46    /// Chaigne & Askenfelt's own C4 reference set, at the fundamental asked for.
47    pub fn at(f0: f64) -> ChaigneAskenfeltParams {
48        ChaigneAskenfeltParams {
49            f0,
50            b: 0.00021,
51            strike_pos: 0.125,
52            vel: 3.2,
53            hammer_mass: 2.9e-3,
54            hammer_k: 2.6646e8,
55            hammer_p: 2.5,
56            damp_dc: 0.6,
57            damp_freq: 1.6e-4,
58            unison_count: 1.0,
59            detune: 2f64.powf(3.0 / 1200.0),
60            bridge_coupling: 1000.0,
61            bridge_mass: 0.0,
62            string_cents: [0.0; MAX_UNISON],
63            string_hammer_k_ratio: [1.0; MAX_UNISON],
64            damper_pos: DAMPER_POS,
65            damper_r: 0.0,
66            damper_k: 0.0,
67        }
68    }
69}
70
71impl ChaigneAskenfeltParams {
72    pub fn valid(&self) -> bool {
73        let [c1, c2, c3] = self.string_cents;
74        let [k1, k2, k3] = self.string_hammer_k_ratio;
75        all(&[
76            (self.f0, Positive),
77            (self.b, NonNegative),
78            (self.strike_pos, OpenUnit),
79            (self.vel, Positive),
80            (self.hammer_mass, Positive),
81            (self.hammer_k, Positive),
82            (self.hammer_p, Positive),
83            (self.damp_dc, NonNegative),
84            (self.damp_freq, NonNegative),
85            (self.unison_count, Within(1.0, MAX_UNISON as f64)),
86            (self.detune, AtLeast(1.0)),
87            (self.bridge_coupling, Positive),
88            (self.bridge_mass, NonNegative),
89            (c1, Finite),
90            (c2, Finite),
91            (c3, Finite),
92            (k1, Positive),
93            (k2, Positive),
94            (k3, Positive),
95            (self.damper_pos, OpenUnit),
96            (self.damper_r, NonNegative),
97            (self.damper_k, NonNegative),
98        ])
99    }
100}
101
102/// The dashpot moves every sample; the spring stores energy, so it only jumps.
103pub const VARYING: &[(&str, Varies)] = &[
104    ("damper_r", Varies::PerSample),
105    ("damper_k", Varies::Piecewise),
106];
107
108/// Fitted to MAPS ENSTDkCl forte key-off slopes (Emiya, Badeau & David 2010) beside a caller's
109/// `damper_r = 0.1 (262/f0)^2` ramped in over 0.03 s; unpressed, the felt is lifted.
110const DAMPER_POS: f64 = 0.15;
111const FELT_LENGTH_M: f64 = 0.04;
112
113const WIRE_DENSITY_KG_M3: f64 = 7850.0;
114/// A generic wire radius, tuned against this module's bridge-force and contact-time tests.
115const WIRE_RADIUS_M: f64 = 0.6e-3;
116/// A generic string tension, tuned alongside [`WIRE_RADIUS_M`].
117const STRING_TENSION_N: f64 = 1500.0;
118/// `L = c/(2 f0)` on a generic wire.
119fn build_grid(params: &ChaigneAskenfeltParams, f0: f64, sr: f64) -> Option<StringGrid> {
120    let rho = std::f64::consts::PI * WIRE_RADIUS_M * WIRE_RADIUS_M * WIRE_DENSITY_KG_M3;
121    let c = (STRING_TENSION_N / rho).sqrt();
122    let length = c / (2.0 * f0);
123    dispersive_grid(
124        Wire { rho, c, length },
125        f0,
126        params.b,
127        params.damp_dc,
128        params.damp_freq,
129        sr,
130    )
131}
132
133/// Strings placed symmetrically in log frequency around `f0` at `+-sqrt(detune)`, then each
134/// moved by its own `cents`.
135fn unison_frequencies(f0: f64, detune: f64, unison_count: usize, cents: &[f64]) -> Vec<f64> {
136    let spread = detune.sqrt();
137    let placed = match unison_count {
138        1 => vec![f0],
139        2 => vec![f0 / spread, f0 * spread],
140        _ => vec![f0 / spread, f0, f0 * spread],
141    };
142    placed
143        .iter()
144        .zip(cents)
145        .map(|(&f, &c)| f * 2f64.powf(c / 1200.0))
146        .collect()
147}
148
149/// Bounds `valid()`: the force buffer is a stack array.
150const MAX_UNISON: usize = 3;
151
152#[derive(Clone)]
153pub struct ChaigneAskenfeltSite {
154    pub(crate) strings: Vec<StringGrid>,
155    /// Per string, `rho (lambda dx/dt)^2`.
156    pub(crate) tensions: Vec<f64>,
157    strike: Vec<Vec<f64>>,
158    hammer: Hammer,
159    detached: Vec<bool>,
160    pub(crate) bridge_now: f64,
161    pub(crate) bridge_prev: f64,
162    bridge_coupling: f64,
163    pub(crate) bridge_mass: f64,
164    pub(crate) dt: f64,
165    pub(crate) felt: Vec<Vec<Felt>>,
166    /// Per string, nodes under the felt, each share, `dt/(rho dx)`.
167    under_felt: Vec<(Vec<usize>, f64, f64)>,
168    /// `(damper_r, damper_k)` as opened, and as the felt is designed for.
169    own: (f64, f64),
170    designed_for: (f64, f64),
171    pub(crate) steps: u64,
172}
173
174/// Every node under the felt, or the nearest, sharing it evenly.
175fn under_felt(
176    grid: &StringGrid,
177    params: &ChaigneAskenfeltParams,
178    dt: f64,
179) -> (Vec<usize>, f64, f64) {
180    let n = grid.n;
181    let at = params.damper_pos * n as f64;
182    let half = FELT_LENGTH_M / 2.0 / grid.dx;
183    let mut nodes: Vec<usize> = (1..n).filter(|&j| (j as f64 - at).abs() <= half).collect();
184    if nodes.is_empty() {
185        nodes.push((at.round() as usize).clamp(1, n - 1));
186    }
187    let psi = 1.0 / nodes.len() as f64;
188    (nodes, psi, dt / (grid.rho * grid.dx))
189}
190
191fn designed(under: &[(Vec<usize>, f64, f64)], dt: f64, (r, k): (f64, f64)) -> Vec<Vec<Felt>> {
192    under
193        .iter()
194        .map(|(nodes, psi, unit)| {
195            nodes
196                .iter()
197                .map(|&j| (j, k * psi * dt * unit, r * psi * unit / 2.0))
198                .collect()
199        })
200        .collect()
201}
202
203impl ChaigneAskenfeltSite {
204    pub fn new(params: &ChaigneAskenfeltParams, sr: f64) -> Result<Self, SampleError> {
205        let unison_count = params.unison_count.round().clamp(1.0, MAX_UNISON as f64) as usize;
206        let freqs =
207            unison_frequencies(params.f0, params.detune, unison_count, &params.string_cents);
208        let strings = freqs
209            .iter()
210            .map(|&f0| build_grid(params, f0, sr))
211            .collect::<Option<Vec<StringGrid>>>()
212            .ok_or(SampleError::StringPastRate {
213                model: "chaigne_askenfelt",
214            })?;
215        let dt = 1.0 / sr;
216        let tensions = strings.iter().map(|g| grid_tension(g, dt)).collect();
217        let strike = strings
218            .iter()
219            .map(|g| point_weights(g.n, params.strike_pos))
220            .collect();
221        let under: Vec<_> = strings.iter().map(|g| under_felt(g, params, dt)).collect();
222        let own = (params.damper_r, params.damper_k);
223        let site = ChaigneAskenfeltSite {
224            felt: designed(&under, dt, own),
225            under_felt: under,
226            own,
227            designed_for: own,
228            steps: 0,
229            detached: vec![false; strings.len()],
230            strings,
231            tensions,
232            strike,
233            hammer: Hammer::new(
234                params.hammer_mass,
235                params.hammer_k,
236                params.hammer_p,
237                params.vel,
238                dt,
239            )
240            .with_anvil_ratios(&params.string_hammer_k_ratio[..unison_count]),
241            bridge_now: 0.0,
242            bridge_prev: 0.0,
243            bridge_coupling: params.bridge_coupling,
244            bridge_mass: params.bridge_mass,
245            dt,
246        };
247        match site.strings.len() == 1 || unison_stable(&site) {
248            true => Ok(site),
249            false => Err(SampleError::BridgeUnstable {
250                model: "chaigne_askenfelt",
251            }),
252        }
253    }
254}
255
256impl ChaigneAskenfeltSite {
257    pub(crate) fn advance(&mut self) -> f64 {
258        // One string ends on a rigid pin; a unison ends on its shared bridge.
259        let sample = match self.strings.len() {
260            1 => self.step_single(),
261            _ => self.step_unison(),
262        };
263        self.steps += 1;
264        sample
265    }
266}
267
268impl ChaigneAskenfeltSite {
269    /// No later step is driven.
270    pub fn let_go(&self) -> bool {
271        self.detached.iter().all(|d| *d)
272    }
273
274    /// The discrete energy in joules: the strings', and a unison's bridge's.
275    pub fn energy(&self) -> f64 {
276        match self.strings.as_slice() {
277            [grid] => energy(grid, self.dt, self.springs(0)).0,
278            _ => unison_energy(self).0,
279        }
280    }
281
282    /// The bridge's dashpot `R_B`, `bridge_coupling sqrt(T rho)` of the first string.
283    pub(crate) fn bridge_r(&self) -> f64 {
284        self.bridge_coupling * (self.tensions[0] * self.strings[0].rho).sqrt()
285    }
286
287    /// [`Self::bridge_r`] as the exact real its stored operands make.
288    pub(crate) fn bridge_r_enclosed(&self) -> Ball {
289        let z = Ball::exact(self.tensions[0])
290            .scale(self.strings[0].rho)
291            .sqrt();
292        z.expect("a positive tension").scale(self.bridge_coupling)
293    }
294
295    pub(crate) fn springs(&self, i: usize) -> &[Felt] {
296        &self.felt[i]
297    }
298
299    /// With neither dashpot nor spring the felt leaves the string, to the bit.
300    pub(crate) fn pressing(&self) -> bool {
301        self.designed_for.0.to_bits() != 0 || self.designed_for.1.to_bits() != 0
302    }
303}
304
305impl Solver for ChaigneAskenfeltSite {
306    fn step(&mut self, args: &[f64]) -> Result<f64, SampleError> {
307        let now = (
308            args.first().copied().unwrap_or(self.own.0),
309            args.get(1).copied().unwrap_or(self.own.1),
310        );
311        for ((name, _), value) in VARYING.iter().zip([now.0, now.1]) {
312            if !(value >= 0.0 && value.is_finite()) {
313                return Err(SampleError::ArgumentOutOfRange {
314                    model: "chaigne_askenfelt",
315                    name,
316                    bits: value.to_bits(),
317                    sample: self.steps,
318                });
319            }
320        }
321        let bits = |(r, k): (f64, f64)| (r.to_bits(), k.to_bits());
322        if bits(now) != bits(self.designed_for) {
323            self.felt = designed(&self.under_felt, self.dt, now);
324            self.designed_for = now;
325        }
326        Ok(self.advance())
327    }
328
329    fn bytes(&self) -> usize {
330        let grids: usize = self.strings.iter().map(StringGrid::bytes).sum();
331        let strike: usize = self.strike.iter().map(|w| super::floats(w)).sum();
332        let felt: usize = self.felt.iter().map(std::mem::size_of_val).sum::<usize>()
333            + self
334                .under_felt
335                .iter()
336                .map(|(nodes, ..)| std::mem::size_of_val(nodes.as_slice()))
337                .sum::<usize>();
338        size_of::<Self>() + grids + strike + felt + super::floats(&self.tensions)
339    }
340}
341
342impl ChaigneAskenfeltSite {
343    /// A released anvil is not read.
344    fn hammer_forces(&mut self, forces: &mut [f64]) {
345        let mut y_h = [0.0f64; MAX_UNISON];
346        for (i, grid) in self.strings.iter().enumerate() {
347            if !self.detached[i] {
348                y_h[i] = read_at(&self.strike[i], &grid.y_now);
349            }
350        }
351        self.hammer.substeps(
352            self.dt,
353            &y_h[..self.strings.len()],
354            &mut self.detached,
355            forces,
356        );
357    }
358
359    fn step_single(&mut self) -> f64 {
360        let mut forces = [0.0f64; 1];
361        self.hammer_forces(&mut forces);
362        let pressing = self.pressing();
363        let grid = &mut self.strings[0];
364        let n = grid.n;
365        for i in 1..n {
366            grid.y_next[i] = stencil_update(grid, i, 0.0, 0.0);
367        }
368        spread(grid, &self.strike[0], forces[0], self.dt);
369        if pressing {
370            press(grid, &self.felt[0]);
371        }
372        grid.y_next[0] = 0.0;
373        grid.y_next[n] = 0.0;
374
375        let sample = self.tensions[0] * (grid.y_now[n] - grid.y_now[n - 1]) / grid.dx;
376
377        std::mem::swap(&mut grid.y_prev, &mut grid.y_now);
378        std::mem::swap(&mut grid.y_now, &mut grid.y_next);
379
380        sample
381    }
382
383    /// The hammer's reaction is summed over the strings, not divided.
384    fn step_unison(&mut self) -> f64 {
385        let mut forces = [0.0f64; MAX_UNISON];
386        let forces = &mut forces[..self.strings.len()];
387        self.hammer_forces(forces);
388        let pressing = self.pressing();
389
390        let (bridge_now, bridge_prev) = (self.bridge_now, self.bridge_prev);
391        let dt2 = self.dt * self.dt;
392        // `M a + R_B v = net string force`: tension implicit in p+, bending and b-loss at p^n.
393        let mut k_eff = 0.0;
394        let mut rhs_sum = 0.0;
395        let mut sample = 0.0;
396        for (i, &force) in forces.iter().enumerate() {
397            let grid = &mut self.strings[i];
398            let tension = self.tensions[i];
399            let n = grid.n;
400            for j in 1..n {
401                grid.y_next[j] = stencil_update(grid, j, bridge_now, bridge_prev);
402            }
403            spread(grid, &self.strike[i], force, self.dt);
404            if pressing {
405                press(grid, &self.felt[i]);
406            }
407            grid.y_next[0] = 0.0;
408            let w = grid.rho * grid.dx / dt2;
409            let (y, yp) = (&grid.y_now, &grid.y_prev);
410            k_eff += tension / grid.dx;
411            rhs_sum += tension * y[n - 1] / grid.dx
412                + w * grid.stiff_sq * (2.0 * y[n - 1] - y[n - 2] - bridge_now)
413                + w * grid.damp_b * ((y[n - 1] - yp[n - 1]) - (bridge_now - bridge_prev));
414            sample += tension * (grid.y_now[n] - grid.y_now[n - 1]) / grid.dx;
415        }
416
417        let r_over_dt = self.bridge_r() / self.dt;
418        let m_over_dt2 = self.bridge_mass / (self.dt * self.dt);
419        let bridge_next =
420            (rhs_sum + r_over_dt * bridge_now + m_over_dt2 * (2.0 * bridge_now - bridge_prev))
421                / (k_eff + r_over_dt + m_over_dt2);
422        for grid in &mut self.strings {
423            let n = grid.n;
424            grid.y_next[n] = bridge_next;
425        }
426
427        self.bridge_prev = bridge_now;
428        self.bridge_now = bridge_next;
429
430        for grid in &mut self.strings {
431            std::mem::swap(&mut grid.y_prev, &mut grid.y_now);
432            std::mem::swap(&mut grid.y_now, &mut grid.y_next);
433        }
434
435        sample
436    }
437}
sva_samples/physics/chaigne_askenfelt.rs

sva_samples/physics/
chaigne_askenfelt.rs
