1// The felt of the stand as a viscoelastic material: a spring in parallel
2// with a chain of Maxwell arms (a spring and a dashpot in series), whose
3// relaxation times are spaced two per decade from 0.5 Hz to above the audio
4// band and whose weights grow geometrically. Such a chain has a loss factor
5// that is nearly flat over the band, which is what a single loss factor
6// eta means, and unlike a dashpot it is causal and can be stepped in time:
7// each arm carries one state per node, its extension p, which relaxes
8// toward zero at its own rate while the node moves.
9//
10// Every quantity here is normalized to the stiffness at the reference
11// frequency (the bounce of the plate on the stand), where the storage
12// stiffness of the chain is 1, so that the felt modulus given by the user is
13// the modulus at the frequency where it matters most. Below the reference
14// the felt is softer and above it stiffer, as the Kramers-Kronig relations
15// require of any material with a loss; at rest its stiffness is kInf.
17export interface FeltChain {
18 /** Relaxation time of each arm [s]. */
19 tau: Float64Array;
20 /** Stiffness of each arm, per unit of the stiffness at the reference frequency. */
21 k: Float64Array;
22 /** Relaxed (static) stiffness, per unit of the stiffness at the reference frequency. */
23 kInf: number;
24 /** The loss factor the chain was built for. */
25 eta: number;
26}
28/** Lowest arm frequency [Hz]; the highest is FELT_ARMS_PER_DECADE arms per decade above it. */
29export const FELT_LOWEST_HZ = 0.5;
30export const FELT_ARMS_PER_DECADE = 2;
31/** Number of arms: 0.5 Hz to 158 kHz. */
32export const FELT_ARMS = 12;
34/**
35 * The chain for a loss factor eta, normalized at the frequency fRef [Hz]. A
36 * chain of arms with frequencies in the ratio r and stiffnesses in the ratio
37 * 1 + c approximates the power law K(i omega) ~ (i omega)^(2 delta / pi) of a
38 * material with the constant loss angle delta = atan(eta), whose modulus
39 * rises by the factor r^(2 delta / pi) per arm; c is that factor minus one.
40 */
41export function feltChain(eta: number, fRef: number): FeltChain {
42 const r = Math.pow(10, 1 / FELT_ARMS_PER_DECADE);
43 const c = eta > 0 ? Math.pow(r, (2 * Math.atan(eta)) / Math.PI) - 1 : 0;
44 const nArms = c > 0 ? FELT_ARMS : 0;
45 const tau = new Float64Array(nArms), k = new Float64Array(nArms);
46 let S = 1;
47 for (let i = 0; i < nArms; i++) {
48 tau[i] = 1 / (2 * Math.PI * FELT_LOWEST_HZ * Math.pow(r, i));
49 k[i] = c * S;
50 S *= 1 + c;
51 }
52 const raw: FeltChain = { tau, k, kInf: 1, eta };
53 const scale = 1 / feltModulus(raw, 2 * Math.PI * Math.max(fRef, 1e-6)).storage;
54 for (let i = 0; i < nArms; i++) k[i] *= scale;
55 return { tau, k, kInf: scale, eta };
56}
58/** Storage and loss stiffness of the chain at the angular frequency omega, per unit of the reference stiffness. */
59export function feltModulus(ch: FeltChain, omega: number): { storage: number; loss: number } {
60 let storage = ch.kInf, loss = 0;
61 for (let i = 0; i < ch.tau.length; i++) {
62 const x = omega * ch.tau[i], d = 1 + x * x;
63 storage += (ch.k[i] * x * x) / d;
64 loss += (ch.k[i] * x) / d;
65 }
66 return { storage, loss };
67}
69/**
70 * The constants of the exact time step of the arms for a step dt: with the
71 * node moving linearly by dw over the step, the extension of arm i goes from
72 * p to p decay_i + dw growth_i, where decay_i = exp(-dt / tau_i) and
73 * growth_i = (tau_i / dt)(1 - decay_i). The force of the arms on the node is
74 * -k_node sum_i k_i p_i.
75 */
76export function feltStepConstants(ch: FeltChain, dt: number): { decay: Float64Array; growth: Float64Array } {
77 const n = ch.tau.length;
78 const decay = new Float64Array(n), growth = new Float64Array(n);
79 for (let i = 0; i < n; i++) {
80 const x = dt / ch.tau[i];
81 decay[i] = Math.exp(-x);
82 growth[i] = -Math.expm1(-x) / x;
83 }
84 return { decay, growth };
85}