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::physics::Solver;
7
8use crate::physics::bound::Bound::*;
9use crate::physics::bound::all;
10use crate::physics::hammer::Hammer;
11
12#[derive(Clone, Debug, PartialEq)]
13pub struct ChaigneAskenfeltParams {
14    pub f0: f64,
15    pub b: f64,
16    pub strike_pos: f64,
17    pub vel: f64,
18    pub hammer_mass: f64,
19    pub hammer_k: f64,
20    pub hammer_p: f64,
21    pub damp_dc: f64,
22    pub damp_freq: f64,
23    pub unison_count: f64,
24    pub detune: f64,
25    pub bridge_coupling: f64,
26}
27
28impl ChaigneAskenfeltParams {
29    /// Chaigne & Askenfelt's own C4 reference set, at the fundamental asked for.
30    pub fn at(f0: f64) -> ChaigneAskenfeltParams {
31        ChaigneAskenfeltParams {
32            f0,
33            b: 0.00021,
34            strike_pos: 0.125,
35            vel: 3.2,
36            hammer_mass: 2.9e-3,
37            hammer_k: 2.6646e8,
38            hammer_p: 2.5,
39            damp_dc: 0.6,
40            damp_freq: 1.6e-4,
41            unison_count: 1.0,
42            detune: 2f64.powf(3.0 / 1200.0),
43            bridge_coupling: 1000.0,
44        }
45    }
46}
47
48impl ChaigneAskenfeltParams {
49    pub fn valid(&self) -> bool {
50        all(&[
51            (self.f0, Positive),
52            (self.b, NonNegative),
53            (self.strike_pos, OpenUnit),
54            (self.vel, Positive),
55            (self.hammer_mass, Positive),
56            (self.hammer_k, Positive),
57            (self.hammer_p, Positive),
58            (self.damp_dc, NonNegative),
59            (self.damp_freq, NonNegative),
60            (self.unison_count, Within(1.0, MAX_UNISON as f64)),
61            (self.detune, AtLeast(1.0)),
62            (self.bridge_coupling, Positive),
63        ])
64    }
65}
66
67const WIRE_DENSITY_KG_M3: f64 = 7850.0;
68/// A generic wire radius, tuned against this module's bridge-force and contact-time tests.
69const WIRE_RADIUS_M: f64 = 0.6e-3;
70/// A generic string tension, tuned alongside [`WIRE_RADIUS_M`].
71const STRING_TENSION_N: f64 = 1500.0;
72/// `SharedBridge`: terminated at the site's common bridge, not a rigid pin. Set uniformly
73/// across every string in a site; `step()` reads only `strings[0]`'s to dispatch.
74#[derive(Clone, Copy)]
75enum Termination {
76    Rigid,
77    SharedBridge,
78}
79
80pub(crate) struct StringGrid {
81    pub(crate) y_now: Vec<f64>,
82    pub(crate) y_prev: Vec<f64>,
83    pub(crate) y_next: Vec<f64>,
84    pub(crate) n: usize,
85    pub(crate) dx: f64,
86    pub(crate) rho: f64,
87    pub(crate) courant_sq: f64,
88    pub(crate) stiff_sq: f64,
89    pub(crate) damp_a: f64,
90    pub(crate) damp_b: f64,
91    far_termination: Termination,
92}
93
94/// The ghost a 5-point biharmonic stencil needs past `0` and `n`, far end `pin`.
95pub(crate) fn ghost_pinned(y: &[f64], n: usize, idx: isize, pin: f64) -> f64 {
96    if idx < 0 {
97        -y[(-idx) as usize]
98    } else if idx as usize > n {
99        2.0 * pin - y[(2 * n as isize - idx) as usize]
100    } else {
101        y[idx as usize]
102    }
103}
104
105/// `pin` is the far-end ghost: `0.0` rigid, or a shared bridge's own state.
106pub(crate) fn stencil_update(grid: &StringGrid, j: usize, pin_now: f64, pin_prev: f64) -> f64 {
107    let n = grid.n;
108    let jj = j as isize;
109    let lap_now = ghost_pinned(&grid.y_now, n, jj + 1, pin_now) - 2.0 * grid.y_now[j]
110        + ghost_pinned(&grid.y_now, n, jj - 1, pin_now);
111    let biharm = ghost_pinned(&grid.y_now, n, jj + 2, pin_now)
112        - 4.0 * ghost_pinned(&grid.y_now, n, jj + 1, pin_now)
113        + 6.0 * grid.y_now[j]
114        - 4.0 * ghost_pinned(&grid.y_now, n, jj - 1, pin_now)
115        + ghost_pinned(&grid.y_now, n, jj - 2, pin_now);
116    let lap_prev = ghost_pinned(&grid.y_prev, n, jj + 1, pin_prev) - 2.0 * grid.y_prev[j]
117        + ghost_pinned(&grid.y_prev, n, jj - 1, pin_prev);
118    2.0 * grid.y_now[j] - grid.y_prev[j] + grid.courant_sq * lap_now
119        - grid.stiff_sq * biharm
120        - grid.damp_a * (grid.y_now[j] - grid.y_prev[j])
121        + grid.damp_b * (lap_now - lap_prev)
122}
123
124pub(crate) fn stiff_string_grid(
125    rho: f64,
126    c: f64,
127    length: f64,
128    kappa: f64,
129    damp_dc: f64,
130    damp_freq: f64,
131    sr: f64,
132) -> StringGrid {
133    let dt = 1.0 / sr;
134    let (c2, dt2, k2) = (c * c, dt * dt, kappa * kappa);
135    let dx_bound = ((c2 * dt2 + (c2 * c2 * dt2 * dt2 + 16.0 * k2 * dt2).sqrt()) / 2.0).sqrt();
136    let n = ((length / dx_bound).floor() as usize).max(4);
137    let dx = length / n as f64;
138
139    StringGrid {
140        y_now: vec![0.0; n + 1],
141        y_prev: vec![0.0; n + 1],
142        y_next: vec![0.0; n + 1],
143        n,
144        dx,
145        rho,
146        courant_sq: c2 * dt2 / (dx * dx),
147        stiff_sq: k2 * dt2 / dx.powi(4),
148        damp_a: 2.0 * damp_dc * dt,
149        damp_b: 2.0 * damp_freq * dt / (dx * dx),
150        far_termination: Termination::Rigid,
151    }
152}
153
154/// Derives a concrete grid from `(f0, b)`: fixes a generic wire's `c`, scales this note's own
155/// length `L = c/(2 f0)`, its stiffness `kappa = c L sqrt(b)/pi`, and a CFL-stable `dx`.
156fn build_grid(params: &ChaigneAskenfeltParams, f0: f64, sr: f64) -> StringGrid {
157    let rho = std::f64::consts::PI * WIRE_RADIUS_M * WIRE_RADIUS_M * WIRE_DENSITY_KG_M3;
158    let c = (STRING_TENSION_N / rho).sqrt();
159    let length = c / (2.0 * f0);
160    let kappa = c * length * params.b.sqrt() / std::f64::consts::PI;
161    stiff_string_grid(rho, c, length, kappa, params.damp_dc, params.damp_freq, sr)
162}
163
164/// Strings placed symmetrically in log frequency around `f0` at `+-sqrt(detune)`.
165fn unison_frequencies(f0: f64, detune: f64, unison_count: usize) -> Vec<f64> {
166    let spread = detune.sqrt();
167    match unison_count {
168        1 => vec![f0],
169        2 => vec![f0 / spread, f0 * spread],
170        _ => vec![f0 / spread, f0, f0 * spread],
171    }
172}
173
174/// Bounds `valid()`, so the per-string force buffer can be a stack array.
175const MAX_UNISON: usize = 3;
176
177pub struct ChaigneAskenfeltSite {
178    strings: Vec<StringGrid>,
179    hammer: Hammer,
180    detached: bool,
181    contact_index: usize,
182    contact_indices: Vec<usize>,
183    strings_detached: Vec<bool>,
184    bridge_now: f64,
185    bridge_prev: f64,
186    bridge_coupling: f64,
187    tension: f64,
188    dt: f64,
189}
190
191impl ChaigneAskenfeltSite {
192    pub fn new(params: &ChaigneAskenfeltParams, sr: f64) -> ChaigneAskenfeltSite {
193        let unison_count = params.unison_count.round().clamp(1.0, MAX_UNISON as f64) as usize;
194        let freqs = unison_frequencies(params.f0, params.detune, unison_count);
195        let mut strings: Vec<StringGrid> =
196            freqs.iter().map(|&f0| build_grid(params, f0, sr)).collect();
197        if strings.len() > 1 {
198            for grid in &mut strings {
199                grid.far_termination = Termination::SharedBridge;
200            }
201        }
202        let contact_indices: Vec<usize> = strings
203            .iter()
204            .map(|g| {
205                (params.strike_pos * g.n as f64)
206                    .round()
207                    .clamp(1.0, (g.n - 1) as f64) as usize
208            })
209            .collect();
210        let contact_index = contact_indices[0];
211        let strings_detached = vec![false; strings.len()];
212        let dt = 1.0 / sr;
213        ChaigneAskenfeltSite {
214            strings,
215            hammer: Hammer::new(
216                params.hammer_mass,
217                params.hammer_k,
218                params.hammer_p,
219                params.vel,
220                dt,
221            ),
222            detached: false,
223            contact_index,
224            contact_indices,
225            strings_detached,
226            bridge_now: 0.0,
227            bridge_prev: 0.0,
228            bridge_coupling: params.bridge_coupling,
229            tension: STRING_TENSION_N,
230            dt,
231        }
232    }
233}
234
235impl Solver for ChaigneAskenfeltSite {
236    fn step(&mut self) -> f64 {
237        match self.strings[0].far_termination {
238            Termination::Rigid => self.step_single(),
239            Termination::SharedBridge => self.step_unison(),
240        }
241    }
242}
243
244impl ChaigneAskenfeltSite {
245    /// The contact point is frozen at this sample's start.
246    fn step_single(&mut self) -> f64 {
247        let grid = &mut self.strings[0];
248        let n = grid.n;
249        let y_h = grid.y_now[self.contact_index];
250
251        let mut forces = [0.0f64; 1];
252        self.hammer.substeps(
253            self.dt,
254            &[y_h],
255            std::slice::from_mut(&mut self.detached),
256            &mut forces,
257        );
258        let force = forces[0];
259
260        for i in 1..n {
261            let mut next = stencil_update(grid, i, 0.0, 0.0);
262            if i == self.contact_index {
263                next += (self.dt * self.dt / (grid.rho * grid.dx)) * force;
264            }
265            grid.y_next[i] = next;
266        }
267        grid.y_next[0] = 0.0;
268        grid.y_next[n] = 0.0;
269
270        let sample = self.tension * (grid.y_now[n] - grid.y_now[n - 1]) / grid.dx;
271
272        std::mem::swap(&mut grid.y_prev, &mut grid.y_now);
273        std::mem::swap(&mut grid.y_now, &mut grid.y_next);
274
275        sample
276    }
277
278    /// One hammer against every contact point: the reaction is summed, not divided.
279    fn step_unison(&mut self) -> f64 {
280        let count = self.strings.len();
281
282        let y_h: Vec<f64> = (0..count)
283            .map(|i| self.strings[i].y_now[self.contact_indices[i]])
284            .collect();
285
286        let mut forces = [0.0f64; MAX_UNISON];
287        let forces = &mut forces[..count];
288        self.hammer
289            .substeps(self.dt, &y_h, &mut self.strings_detached, forces);
290
291        let (bridge_now, bridge_prev) = (self.bridge_now, self.bridge_prev);
292        // Massless resistive bridge: R_B*v = net string force, solved as a stable weighted average.
293        let mut k_eff = 0.0;
294        let mut rhs_sum = 0.0;
295        let mut sample = 0.0;
296        for (i, &force) in forces.iter().enumerate() {
297            let grid = &mut self.strings[i];
298            let n = grid.n;
299            for j in 1..n {
300                let mut next = stencil_update(grid, j, bridge_now, bridge_prev);
301                if j == self.contact_indices[i] {
302                    next += (self.dt * self.dt / (grid.rho * grid.dx)) * force;
303                }
304                grid.y_next[j] = next;
305            }
306            grid.y_next[0] = 0.0;
307            k_eff += self.tension / grid.dx;
308            rhs_sum += self.tension * grid.y_now[n - 1] / grid.dx;
309            sample += self.tension * (grid.y_now[n] - grid.y_now[n - 1]) / grid.dx;
310        }
311
312        let z_string = (self.tension * self.strings[0].rho).sqrt();
313        let r_bridge = self.bridge_coupling * z_string;
314        let r_over_dt = r_bridge / self.dt;
315        let bridge_next = (rhs_sum + r_over_dt * bridge_now) / (k_eff + r_over_dt);
316        for grid in &mut self.strings {
317            let n = grid.n;
318            grid.y_next[n] = bridge_next;
319        }
320
321        self.bridge_prev = bridge_now;
322        self.bridge_now = bridge_next;
323
324        for grid in &mut self.strings {
325            std::mem::swap(&mut grid.y_prev, &mut grid.y_now);
326            std::mem::swap(&mut grid.y_now, &mut grid.y_next);
327        }
328
329        sample
330    }
331}
sva_samples/physics/chaigne_askenfelt.rs

sva_samples/physics/
chaigne_askenfelt.rs
