1function G = cymbal_profile(P)
2% CYMBAL_PROFILE Meridian of the cymbal as a shell of revolution.
3%
4% The cymbal is a surface of revolution about the z axis: a central mounting
5% hole of radius hole_radius_m, a raised bell out to bell_radius_m, and a bow
6% falling from the bell to the rim, which sits at z = 0 and r = radius_m.
7% Arclength s runs from the hole (s = 0) to the rim (s = S), so r increases
8% with s and the parameter strike_r addresses the plate as s / S.
9%
10% The bow is a parabola in the rim-referred coordinate u = (a - r) / (a - rb),
11% z = bow_height_m u^2, which is flat at the rim and steepest where it meets
12% the bell. The bell is the quartic in r that matches the bow's height, slope
13% and curvature at r = rb and comes to a horizontal apex at r = 0. Matching
14% the curvature matters: the finite element uses the meridian curvature and
15% its arclength derivative, so a profile that is only C1 would put a delta
16% function in the bending strain at the bell/bow join.
17%
18% Everything the element routine needs is tabulated on a fine grid and
19% interpolated onto the Gauss points:
20% r, z meridian coordinates (m)
21% s arclength from the hole (m)
22% thk shell thickness (m), tapering linearly from bell to rim
23% drds dr/ds = cos(meridian angle)
24% dzds dz/ds
25% cs 1 / R_s, meridian curvature, positive where the surface curves
26% away from the outward normal (a sphere of radius R gives 1/R)
27% ct 1 / R_theta = -dz/ds / r, the circumferential curvature
28% csp, ctp d/ds of those two
29%
30% Returns G with those fields plus S (total arclength), snode (node
31% arclengths), and the node-sampled geometry.
33a = P.radius_m;
34rh = max(1e-4, min(P.hole_radius_m, 0.5 * a));
35rb = min(max(P.bell_radius_m, 2 * rh), 0.8 * a);
36hb = P.bell_height_m;
37hw = P.bow_height_m;
39M = 2001;
40r = linspace(rh, a, M);
42% ---- bow: z = hw u^2 with u = (a - r) / (a - rb) ----
43den = a - rb;
44u = (a - r) / den;
45zb = hw * u.^2;
46zb_r = -2 * hw * u / den;
47zb_rr = (2 * hw / den^2) * ones(1, M);
49% ---- bell: quartic matching z, z' and z'' at rb, horizontal apex at r = 0 ----
50mb = -2 * hw / den; % dz/dr of the bow at r = rb
51kb = 2 * hw / den^2; % d2z/dr2 of the bow at r = rb
52A0 = hw + hb; % apex height (z at r = 0)
53Amat = [rb^2, rb^3, rb^4; 2 * rb, 3 * rb^2, 4 * rb^3; 2, 6 * rb, 12 * rb^2];
54coef = Amat \ [-hb; mb; kb];
55B2 = coef(1); C3 = coef(2); E4 = coef(3);
56zc = A0 + B2 * r.^2 + C3 * r.^3 + E4 * r.^4;
57zc_r = 2 * B2 * r + 3 * C3 * r.^2 + 4 * E4 * r.^3;
58zc_rr = 2 * B2 + 6 * C3 * r + 12 * E4 * r.^2;
60inBell = r < rb;
61z = zb;
62z_r = zb_r;
63z_rr = zb_rr;
64z(inBell) = zc(inBell);
65z_r(inBell) = zc_r(inBell);
66z_rr(inBell) = zc_rr(inBell);
68% ---- arclength and the differential geometry of the meridian ----
69q = sqrt(1 + z_r.^2);
70sarc = cumtrapz(r, q);
71S = sarc(M);
72drds = 1 ./ q;
73dzds = z_r ./ q;
74d2zds2 = z_rr ./ q.^4;
75cs = -d2zds2 ./ drds;
76ct = -dzds ./ r;
77ctp = -d2zds2 ./ r + dzds .* drds ./ r.^2;
79% cs involves the third derivative of the profile; the fine grid is dense
80% enough that a central difference is accurate to ~1e-8 relative
81csp = zeros(1, M);
82csp(2:M - 1) = (cs(3:M) - cs(1:M - 2)) ./ (sarc(3:M) - sarc(1:M - 2));
83csp(1) = csp(2);
84csp(M) = csp(M - 1);
86% ---- thickness: constant through the bell, linear taper across the bow ----
87taper = min(1, max(0, (a - r) / den));
88thk = P.thickness_rim_m + (P.thickness_bell_m - P.thickness_rim_m) * taper;
90% ---- nodes, equally spaced in arclength ----
91nn = max(6, round(P.n_elements)) + 1;
92snode = linspace(0, S, nn);
94G.r = r;
95G.z = z;
96G.s = sarc;
97G.thk = thk;
98G.drds = drds;
99G.dzds = dzds;
100G.cs = cs;
101G.ct = ct;
102G.csp = csp;
103G.ctp = ctp;
104G.S = S;
105G.M = M;
106G.nn = nn;
107G.snode = snode;
108G.rnode = interp1(sarc, r, snode);
109G.znode = interp1(sarc, z, snode);
110G.drdsnode = interp1(sarc, drds, snode);
111G.dzdsnode = interp1(sarc, dzds, snode);
112G.thknode = interp1(sarc, thk, snode);
113G.area = 2 * pi * trapz(sarc, r);
114G.mass = 2 * pi * P.density * trapz(sarc, r .* thk);
115end