sva_formula/modal/
mod.rs

1// Concern: the geometries and the hammer a bank is built from, and the atoms it expands to | Non-concern: any geometry's own eigenfrequencies (the siblings) | IO: (&ModalBank) -> Vec<SpectralAtom>
2
3pub mod bar;
4pub mod bore;
5pub mod damping;
6mod grid;
7pub mod membrane;
8pub mod room;
9pub mod string;
10
11use std::f64::consts::TAU;
12
13use crate::closed_form::{Edge, Excitation, ModalBank, Mode};
14use crate::complex::C64;
15use crate::modal::damping::Damping;
16use crate::origin::Origin;
17use crate::spectral_sum::atom::{Exp, Factors, Indicator, Poly, Singular, SpectralAtom};
18use crate::spectral_sum::image::raised_cosine;
19
20/// An impulse rings each mode from `t0`; a hammer's force rings it by convolution.
21pub fn atoms(bank: &ModalBank, origin: Origin) -> Vec<SpectralAtom> {
22    let rings = bank.modes.iter().flat_map(|m| {
23        [1.0, -1.0].map(|sign| {
24            let onset = C64::new(0.0, sign * m.phase).exp().scale(m.amp * 0.5);
25            (onset, C64::new(-1.0 / m.tau, sign * m.omega))
26        })
27    });
28    match bank.excite {
29        Excitation::Impulse { t0 } => rings
30            .map(|(a, alpha)| {
31                onset_atom(
32                    a,
33                    alpha,
34                    t0,
35                    (Edge::at(t0), Edge::PosInf),
36                    Poly::ONE,
37                    origin,
38                )
39            })
40            .collect(),
41        Excitation::HammerPulse { f0, t0, contact } => {
42            struck(rings.collect(), Force { f0, t0, contact }, origin)
43        }
44    }
45}
46
47/// The felt's force.
48struct Force {
49    f0: f64,
50    t0: f64,
51    contact: f64,
52}
53
54impl Force {
55    fn lines(&self) -> [(C64, C64); 3] {
56        let w = TAU / self.contact;
57        let side = C64::real(-self.f0 * 0.25);
58        [
59            (C64::real(self.f0 * 0.5), C64::ZERO),
60            (side, C64::new(0.0, w)),
61            (side, C64::new(0.0, -w)),
62        ]
63    }
64}
65
66/// Series terms where `|d|T <= 1/2`; the rest is under 1e-23 of `|a q| T`.
67const NEAR: usize = 18;
68
69/// `a e^(alpha u)` driven by `q e^(beta u)` over `T`: `a q (e^(beta u) - e^(alpha u)) / d`
70/// while pressed, each term under `2|aq|T` past `|d|T = 1/2`, a series within; then a ring.
71fn struck(rings: Vec<(C64, C64)>, force: Force, origin: Origin) -> Vec<SpectralAtom> {
72    let (t0, contact) = (force.t0, force.contact);
73    if contact <= 0.0 || !contact.is_finite() {
74        return Vec::new();
75    }
76    let t1 = t0 + contact;
77    let pressing = (Edge::at(t0), Edge::at(t1));
78    let lines = force.lines();
79    let mut driven = [C64::ZERO; 3];
80    let mut out = Vec::with_capacity(3 * rings.len() + lines.len());
81    for (a, alpha) in rings {
82        let ringing = alpha.scale(contact).exp();
83        let (mut head, mut released) = (C64::ZERO, C64::ZERO);
84        let mut near = [C64::ZERO; NEAR];
85        for (k, (q, beta)) in lines.iter().enumerate() {
86            let aq = a * *q;
87            let d = *beta - alpha;
88            released = released
89                + aq.scale(contact) * apart(beta.scale(contact), ringing, d.scale(contact));
90            if d.abs() * contact > 0.5 {
91                let over = aq / d;
92                driven[k] = driven[k] + over;
93                head = head - over;
94                continue;
95            }
96            let mut term = aq;
97            for (n, held) in near.iter_mut().enumerate() {
98                *held = *held + term;
99                term = (term * d).over((n + 2) as f64);
100            }
101        }
102        out.push(onset_atom(head, alpha, t0, pressing, Poly::ONE, origin));
103        for (n, c) in near.into_iter().enumerate() {
104            let power = Poly {
105                degree: (n + 1) as u16,
106                at: t0,
107            };
108            out.push(onset_atom(c, alpha, t0, pressing, power, origin));
109        }
110        let after = (Edge::at(t1), Edge::PosInf);
111        out.push(onset_atom(released, alpha, t1, after, Poly::ONE, origin));
112    }
113    for ((_, beta), c) in lines.iter().zip(driven) {
114        out.push(onset_atom(c, *beta, t0, pressing, Poly::ONE, origin));
115    }
116    out.retain(|atom| !atom.c.is_zero());
117    out
118}
119
120/// `(e^(beta T) - e^(alpha T)) / (d T)`, uncancelled near `d = 0`.
121fn apart(beta_t: C64, ringing: C64, d_t: C64) -> C64 {
122    if d_t.re == 0.0 && d_t.im == 0.0 {
123        return ringing;
124    }
125    let lifted = match d_t.re.abs() <= 1.0 && d_t.im.abs() <= 1.0 {
126        true => ringing * d_t.exp_m1(),
127        false => beta_t.exp() - ringing,
128    };
129    lifted / d_t
130}
131
132/// `c e^(alpha (t - t0))` over `[l, r)`, from `t0` where from zero it overflows.
133fn onset_atom(
134    c: C64,
135    alpha: C64,
136    t0: f64,
137    (l, r): (Edge, Edge),
138    poly: Poly,
139    origin: Origin,
140) -> SpectralAtom {
141    let from_zero = c * (-alpha.scale(t0)).exp();
142    let (c, exp) = match (alpha.re == 0.0 && alpha.im == 0.0, from_zero.is_finite()) {
143        (true, _) => (c, None),
144        (false, true) => (from_zero, Some(Exp::at(alpha.re, alpha.im))),
145        (false, false) => (
146            c * C64::new(0.0, -alpha.im * t0).exp(),
147            Some(Exp {
148                sigma: alpha.re,
149                omega: alpha.im,
150                mu: t0,
151            }),
152        ),
153    };
154    SpectralAtom::new(
155        c,
156        Factors {
157            poly,
158            exp,
159            ind: Some(Indicator { l, r }),
160            ..Factors::NONE
161        },
162        Singular::Regular,
163        origin,
164    )
165}
166
167/// Which geometry a modal builtin names.
168pub enum Geometry {
169    String(string::StringGeom),
170    Membrane(membrane::MembraneGeom),
171    Bar(bar::BarGeom),
172    Bore(bore::BoreGeom),
173    Room(room::RoomGeom),
174}
175
176pub fn mode(hz: f64, amp: f64, damping: Damping) -> Mode {
177    Mode {
178        omega: TAU * hz,
179        tau: damping.tau(hz),
180        amp,
181        phase: 0.0,
182    }
183}
184
185pub fn modes(g: &Geometry, count: usize) -> Vec<Mode> {
186    match g {
187        Geometry::String(s) => string::modes(s, count),
188        Geometry::Membrane(m) => membrane::modes(m, count),
189        Geometry::Bar(b) => bar::modes(b, count),
190        Geometry::Bore(b) => bore::modes(b, count),
191        Geometry::Room(r) => room::modes(r, count),
192    }
193}
194
195/// Hertzian felt is `p = 1.5`, piano felt about 2.5; Chaigne & Askenfelt, JASA 95 (1994) 1112.
196pub struct HammerGeom {
197    pub reference_velocity: f64,
198    pub reference_contact: f64,
199    pub stiffness_exponent: f64,
200    pub force: f64,
201}
202
203pub fn contact_time(velocity: f64, h: &HammerGeom) -> f64 {
204    let p = h.stiffness_exponent;
205    h.reference_contact * (velocity / h.reference_velocity).powf((1.0 - p) / (1.0 + p))
206}
207
208/// `F0 * Tc` is the hammer's momentum, so the peak force carries what `Tc` does not.
209pub fn strike(velocity: f64, h: &HammerGeom, at: f64) -> Excitation {
210    let p = h.stiffness_exponent;
211    let ratio = velocity / h.reference_velocity;
212    Excitation::HammerPulse {
213        f0: h.force * ratio.powf(2.0 * p / (p + 1.0)),
214        t0: at,
215        contact: contact_time(velocity, h),
216    }
217}
218
219/// A whole turn of a raised cosine, upside down, over the contact; none over no contact.
220pub fn pulse(f0: f64, t0: f64, contact: f64, origin: Origin) -> Vec<SpectralAtom> {
221    if contact <= 0.0 || !contact.is_finite() {
222        return Vec::new();
223    }
224    let [(level, _), (side, turn), _] = Force { f0, t0, contact }.lines();
225    raised_cosine(t0, t0 + contact, turn.im, level.re, side.re, origin)
226}
sva_formula/modal/mod.rs

sva_formula/modal/
mod.rs
