c5866daReimplement the algorithm in MATLAB syntax, run it with numblJeremy Magland 1function basis = woodblock_basis(P)
2% WOODBLOCK_BASIS Modal basis: free-free Euler-Bernoulli bar.
3% Mass-normalized bending modes and their couplings (woodblock.md section 3).
4% basis.f modal frequencies (Hz), n x 1
5% basis.phiC mode shape at the strike point (1/sqrt(kg))
6% basis.A cavity volume coupling (m^2/sqrt(kg))
7% basis.shapes mass-normalized shapes sampled on basis.s, nGrid x n
9L = P.length_m; W = P.width_m; H = P.wall_thickness_m;
10mass = P.participating_mass_kg;
11if mass == 0
12 mass = P.density * W * H * L;
13end
14cLong = sqrt(P.youngs_modulus / P.density);
15% f_n = (beta_n L)^2 / (2 pi L^2) * sqrt(E I / (rho A)), I = W H^3/12
16prefactor = cLong * H / (2 * pi * sqrt(12) * L^2);
18nGrid = 4096;
19s = linspace(0, 1, nGrid).';
21% Roots of cos(bL)cosh(bL) = 1; asymptotic beyond the 4th
22blTab = [4.73004074 7.85320462 10.99560784 14.13716549];
23f = zeros(P.n_modes, 1);
24phiC = zeros(P.n_modes, 1);
25A = zeros(P.n_modes, 1);
26shapes = zeros(nGrid, P.n_modes);
27f1 = prefactor * blTab(1)^2 * P.hollow_softening;
28n = 0;
29for m = 1:P.n_modes
30 if m <= 4
31 bl = blTab(m);
32 else
33 bl = (2 * m + 1) * pi / 2;
34 end
35 fm = prefactor * bl^2 * P.hollow_softening;
36 if m > 1 && P.inharmonicity ~= 1
37 fm = f1 * (fm / f1)^P.inharmonicity;
38 end
39 if fm > P.max_mode_hz
40 break;
41 end
42 n = n + 1;
43 sigma = (cosh(bl) - cos(bl)) / (sinh(bl) - sin(bl));
44 x = bl * s;
45 shape = (cosh(x) + cos(x)) - sigma * (sinh(x) + sin(x));
46 rmsval = sqrt(trapz(s, shape.^2));
47 shape = shape / (rmsval * sqrt(mass));
48 f(n) = fm;
49 phiC(n) = interp1(s, shape, P.strike_position);
50 mask = (s >= P.cavity_span_lo) & (s <= P.cavity_span_hi);
51 A(n) = W * L * trapz(s, shape .* mask);
52 shapes(:, n) = shape;
53end
54if n == 0
55 error('No modes below max_mode_hz -- check geometry/material.');
56end
57f = f(1:n);
58omega = 2 * pi * f;
59% assemble by assignment (struct() with array values would make a struct array)
60basis.n = n;
61basis.f = f;
62basis.omega = omega;
63basis.zeta = 0.5 * (P.loss_factor + P.loss_slope * omega);
64basis.phiC = phiC(1:n);
65basis.A = A(1:n) * P.cavity_coupling;
66basis.mass = mass;
67basis.s = s;
68basis.shapes = shapes(:, 1:n);
69end