1function basis = cymbal_basis(P)
2% CYMBAL_BASIS Mode set of the cymbal: one shell eigensolve per nodal-diameter
3% order, plus damping and the radiation data the later stages need.
4%
5% Because the plate is a shell of revolution, the modes separate exactly in
6% the azimuth (see cymbal_shell.m), so the two-dimensional eigenproblem breaks
7% into n_diameters + 1 one-dimensional ones. Each is reduced to symmetric
8% standard form through the Cholesky factor of the consistent mass matrix and
9% diagonalized by svd rather than eig: the reduced matrix is symmetric
10% positive semidefinite, so its singular vectors are its eigenvectors, and svd
11% keeps them orthonormal for the closely spaced pairs this geometry produces.
12% The eigenvectors come back mass-normalized, phi' M phi = I. Rigid-body modes
13% (two at n = 0, two at n = 1, none above) sit at zero frequency and are
14% dropped with everything below min_mode_hz; modes above max_mode_hz are
15% truncated, and the broadband part they carry is reinstated by the direct
16% term in cymbal_radiate.m.
17%
18% Only the cos(n theta) member of each degenerate pair is retained. That is
19% exact for the model as posed: every contact sits on the meridian theta = 0,
20% where the sin member has a node. A real cymbal is not perfectly round, which
21% splits each pair into a slowly beating doublet; that splitting is outside
22% this model, and the README says so.
23%
24% Damping has three parts.
25% internal damp_rate / om + 0.5 (loss_factor + loss_slope om). The first
26% term is a frequency-independent decay RATE, which is what a
27% resistive path out of the plate gives and what keeps the low
28% modes of a cymbal ringing for tens of seconds; the other two are
29% the usual constant and frequency-proportional loss factors, and
30% they are what makes the shimmer die well before the hum does;
31% mount a point dashpot of mount_loss N s/m on the axial motion at the
32% inner edge, standing in for the felt washer and the stand. With
33% mass-normalized shapes its modal damping is c phi_axial^2 / 2 om,
34% which falls off steeply with n because the shapes of high
35% nodal-diameter modes vanish at the hole. This is what kills the
36% low n = 0 gong modes a mounted cymbal does not sustain;
37% radiation computed from the actual radiated power, below.
38%
39% Radiated power is evaluated in the far field, where it is defined, rather
40% than by summing element powers as the sibling triangle does (an incoherent
41% element sum is mesh dependent for a continuous surface: halving the element
42% size halves the estimate). The plate is treated as a sheet of normal dipoles
43% of volume dS * ldip, and for a cos(n theta) source on a surface of revolution
44% the azimuthal integral of the far-field kernel is a Bessel function, so the
45% far field is exactly
46% p(th_o, ph_o) = -(om^2 rho0 / (4 pi c R)) e^{-ikR} A_n(th_o) cos(n ph_o)
47% A_n = ldip sum_e r_e L_e W_e e^{i k z_e cos th_o}
48% [ 2 pi i^n J_n(al) r'_e cos th_o
49% - pi (i^{n+1} J_{n+1}(al) + i^{n-1} J_{n-1}(al)) z'_e sin th_o ]
50% al = k r_e sin th_o
51% and the power integral over the sphere collapses to one quadrature in
52% cos th_o. The cos(n ph_o) factor is the model's multipole structure made
53% explicit: a mode with n nodal diameters is a 2^n-pole, so high-n modes
54% radiate very little at low frequency, and none of them radiate on the axis.
55%
56% ldip is the one phenomenological length in the radiation model: the distance
57% over which the air short-circuits from one face of the plate to the other.
58% It is taken as half the shorter of the flexural and acoustic wavelengths,
59% capped at the plate radius, and scaled by dipole_scale. Below the plate's
60% coincidence frequency (about 18 kHz for 1 mm bronze) it is the flexural
61% wavelength that sets it, which is why a cymbal is such a poor radiator in
62% the low register and gets steadily better as the cascade fills the top.
64G = cymbal_profile(P);
65E = P.youngs_modulus;
66nu = P.poisson_ratio;
67rho = P.density;
68c = P.speed_of_sound;
69rho0 = P.air_density;
70nn = G.nn;
71nel = nn - 1;
72Lref = G.S / nel;
73nmax = max(0, round(P.n_diameters));
75% ---- radiation sample points: two per element, Hermite-evaluated ----
76nsub = 2;
77nsurf = nel * nsub;
78surfS = zeros(nsurf, 1);
79surfLen = zeros(nsurf, 1);
80surfEl = zeros(nsurf, 1);
81surfXi = zeros(nsurf, 1);
82for e = 1:nel
83 L = G.snode(e + 1) - G.snode(e);
84 for j = 1:nsub
85 idx = (e - 1) * nsub + j;
86 xi = (j - 0.5) / nsub;
87 surfS(idx) = G.snode(e) + xi * L;
88 surfLen(idx) = L / nsub;
89 surfEl(idx) = e;
90 surfXi(idx) = xi;
91 end
92end
93surfR = reshape(interp1(G.s, G.r, surfS), nsurf, 1);
94surfZ = reshape(interp1(G.s, G.z, surfS), nsurf, 1);
95surfDr = reshape(interp1(G.s, G.drds, surfS), nsurf, 1);
96surfDz = reshape(interp1(G.s, G.dzds, surfS), nsurf, 1);
98% ---- eigensolve per circumferential order ----
99fs = [];
100nd = [];
101NUall = [];
102NWall = [];
103NTall = [];
104SWall = [];
105omLo = 2 * pi * P.min_mode_hz;
106omHi = 2 * pi * P.max_mode_hz;
107for n = 0:nmax
108 [K, M] = cymbal_shell(G, n, E, nu, rho, Lref);
109 Rc = chol(M);
110 Am = (Rc') \ (K / Rc);
111 Am = (Am + Am') / 2;
112 [Uv, Sv, ~] = svd(Am);
113 om2 = flipud(diag(Sv));
114 Vv = fliplr(Uv);
115 om = sqrt(max(0, om2)) * P.stiffness_trim;
116 keep = find(om > omLo & om < omHi);
117 if length(keep) == 0
118 continue;
119 end
120 Phi = Rc \ Vv(:, keep);
121 nk = length(keep);
122 Uz = Phi(1:4:4 * nn, :); % meridional, nn x nk
123 Wz = Phi(3:4:4 * nn, :); % normal
124 Tz = Phi(4:4:4 * nn, :) * Lref; % dW/ds, unscaled back to physical
125 % normal displacement at the radiation sample points (Hermite)
126 SW = zeros(nsurf, nk);
127 for idx = 1:nsurf
128 e = surfEl(idx);
129 xi = surfXi(idx);
130 L = G.snode(e + 1) - G.snode(e);
131 h1 = 1 - 3 * xi^2 + 2 * xi^3;
132 h2 = L * (xi - 2 * xi^2 + xi^3);
133 h3 = 3 * xi^2 - 2 * xi^3;
134 h4 = L * (-xi^2 + xi^3);
135 SW(idx, :) = h1 * Wz(e, :) + h2 * Tz(e, :) + h3 * Wz(e + 1, :) + h4 * Tz(e + 1, :);
136 end
137 fs = [fs; om(keep) / (2 * pi)];
138 nd = [nd; n * ones(nk, 1)];
139 NUall = [NUall, Uz];
140 NWall = [NWall, Wz];
141 NTall = [NTall, Tz];
142 SWall = [SWall, SW];
143end
144N = length(fs);
145if N == 0
146 error('No cymbal modes between min_mode_hz and max_mode_hz -- check geometry/material.');
147end
148[fs, order] = sort(fs);
149fs = reshape(fs, N, 1);
150nd = reshape(nd(order), N, 1);
151NUall = NUall(:, order);
152NWall = NWall(:, order);
153NTall = NTall(:, order);
154SWall = SWall(:, order);
155omega = 2 * pi * fs;
157% ---- the acoustic short-circuit length, per mode ----
158tMean = 2 * pi * trapz(G.s, G.r .* G.thk) / G.area;
159cL = sqrt(E / (rho * (1 - nu^2)));
160kap = tMean / sqrt(12);
161kFlex = sqrt(omega / (kap * cL));
162kAc = omega / c;
163ldip = P.dipole_scale * min(P.radius_m, pi ./ max(kFlex, kAc));
165% ---- radiated power, and the radiation damping that follows from it ----
166nMu = 24;
167mu = reshape(((1:nMu) - 0.5) / nMu * 2 - 1, 1, nMu); % cos(theta_o) midpoints
168dmu = 2 / nMu;
169sth = sqrt(max(0, 1 - mu.^2));
170zetaRad = zeros(N, 1);
171sigmaEff = zeros(N, 1);
172for m = 1:N
173 om = omega(m);
174 k = om / c;
175 nn_d = nd(m);
176 alpha = (k * surfR) * sth; % nsurf x nMu
177 Jn = besselj(nn_d, alpha);
178 Jp1 = besselj(nn_d + 1, alpha);
179 if nn_d == 0
180 Jm1 = -Jp1;
181 else
182 Jm1 = besselj(nn_d - 1, alpha);
183 end
184 % i^(n+1) J_{n+1} + i^(n-1) J_{n-1} = i^n * i * (J_{n+1} - J_{n-1})
185 ipow = 1i^nn_d;
186 Tk = ipow * (2 * pi * Jn .* (surfDr * mu) ...
187 - 1i * pi * (Jp1 - Jm1) .* (surfDz * sth));
188 carrier = exp(1i * (k * surfZ) * mu);
189 wgt = ldip(m) * surfR .* surfLen .* SWall(:, m);
190 An = sum((wgt * ones(1, nMu)) .* Tk .* carrier, 1); % 1 x nMu
191 Kn = sum(abs(An).^2) * dmu;
192 if nn_d == 0
193 piN = 2 * pi;
194 else
195 piN = pi;
196 end
197 zetaRad(m) = om^3 * rho0 * piN * Kn / (32 * pi^2 * c^3);
198 % radiation efficiency, for the diagnostics only
199 vv = sum(surfR .* surfLen .* SWall(:, m).^2) * piN;
200 sigmaEff(m) = 2 * zetaRad(m) * om / max(1e-30, rho0 * c * vv);
201end
203% ---- mount: a point dashpot on the axial motion at the inner edge ----
204axHole = NUall(1, :).' * G.dzdsnode(1) + NWall(1, :).' * G.drdsnode(1);
205zetaMount = P.mount_loss * axHole.^2 ./ (2 * omega);
207zeta = P.damp_rate ./ omega ...
208 + 0.5 * (P.loss_factor + P.loss_slope * omega) ...
209 + P.radiation_loss * zetaRad + zetaMount;
211basis.n = N;
212basis.f = fs;
213basis.omega = omega;
214basis.zeta = zeta;
215basis.zetaRad = zetaRad;
216basis.zetaMount = zetaMount;
217basis.sigma = sigmaEff;
218basis.ndia = nd;
219basis.ldip = ldip;
220basis.NU = NUall;
221basis.NW = NWall;
222basis.NT = NTall;
223basis.rnode = reshape(G.rnode, nn, 1);
224basis.znode = reshape(G.znode, nn, 1);
225basis.snode = reshape(G.snode, nn, 1);
226basis.drdsnode = reshape(G.drdsnode, nn, 1);
227basis.dzdsnode = reshape(G.dzdsnode, nn, 1);
228basis.thknode = reshape(G.thknode, nn, 1);
229basis.csnode = reshape(interp1(G.s, G.cs, G.snode), nn, 1);
230basis.ctnode = reshape(interp1(G.s, G.ct, G.snode), nn, 1);
231basis.surfR = surfR;
232basis.surfZ = surfZ;
233basis.surfDr = surfDr;
234basis.surfDz = surfDz;
235basis.surfLen = surfLen;
236basis.SW = SWall;
237basis.S = G.S;
238basis.area = G.area;
239basis.mass = G.mass;
240basis.thkMean = tMean;
241end