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

sva_samples/physics/
chaigne_askenfelt.rs
