function G = cymbal_profile(P)
% CYMBAL_PROFILE Meridian of the cymbal as a shell of revolution.
%
% The cymbal is a surface of revolution about the z axis: a central mounting
% hole of radius hole_radius_m, a raised bell out to bell_radius_m, and a bow
% falling from the bell to the rim, which sits at z = 0 and r = radius_m.
% Arclength s runs from the hole (s = 0) to the rim (s = S), so r increases
% with s and the parameter strike_r addresses the plate as s / S.
%
% The bow is a parabola in the rim-referred coordinate u = (a - r) / (a - rb),
% z = bow_height_m u^2, which is flat at the rim and steepest where it meets
% the bell. The bell is the quartic in r that matches the bow's height, slope
% and curvature at r = rb and comes to a horizontal apex at r = 0. Matching
% the curvature matters: the finite element uses the meridian curvature and
% its arclength derivative, so a profile that is only C1 would put a delta
% function in the bending strain at the bell/bow join.
%
% Everything the element routine needs is tabulated on a fine grid and
% interpolated onto the Gauss points:
%   r, z      meridian coordinates (m)
%   s         arclength from the hole (m)
%   thk       shell thickness (m), tapering linearly from bell to rim
%   drds      dr/ds = cos(meridian angle)
%   dzds      dz/ds
%   cs        1 / R_s, meridian curvature, positive where the surface curves
%             away from the outward normal (a sphere of radius R gives 1/R)
%   ct        1 / R_theta = -dz/ds / r, the circumferential curvature
%   csp, ctp  d/ds of those two
%
% Returns G with those fields plus S (total arclength), snode (node
% arclengths), and the node-sampled geometry.

a = P.radius_m;
rh = max(1e-4, min(P.hole_radius_m, 0.5 * a));
rb = min(max(P.bell_radius_m, 2 * rh), 0.8 * a);
hb = P.bell_height_m;
hw = P.bow_height_m;

M = 2001;
r = linspace(rh, a, M);

% ---- bow: z = hw u^2 with u = (a - r) / (a - rb) ----
den = a - rb;
u = (a - r) / den;
zb = hw * u.^2;
zb_r = -2 * hw * u / den;
zb_rr = (2 * hw / den^2) * ones(1, M);

% ---- bell: quartic matching z, z' and z'' at rb, horizontal apex at r = 0 ----
mb = -2 * hw / den;          % dz/dr of the bow at r = rb
kb = 2 * hw / den^2;         % d2z/dr2 of the bow at r = rb
A0 = hw + hb;                % apex height (z at r = 0)
Amat = [rb^2, rb^3, rb^4; 2 * rb, 3 * rb^2, 4 * rb^3; 2, 6 * rb, 12 * rb^2];
coef = Amat \ [-hb; mb; kb];
B2 = coef(1); C3 = coef(2); E4 = coef(3);
zc = A0 + B2 * r.^2 + C3 * r.^3 + E4 * r.^4;
zc_r = 2 * B2 * r + 3 * C3 * r.^2 + 4 * E4 * r.^3;
zc_rr = 2 * B2 + 6 * C3 * r + 12 * E4 * r.^2;

inBell = r < rb;
z = zb;
z_r = zb_r;
z_rr = zb_rr;
z(inBell) = zc(inBell);
z_r(inBell) = zc_r(inBell);
z_rr(inBell) = zc_rr(inBell);

% ---- arclength and the differential geometry of the meridian ----
q = sqrt(1 + z_r.^2);
sarc = cumtrapz(r, q);
S = sarc(M);
drds = 1 ./ q;
dzds = z_r ./ q;
d2zds2 = z_rr ./ q.^4;
cs = -d2zds2 ./ drds;
ct = -dzds ./ r;
ctp = -d2zds2 ./ r + dzds .* drds ./ r.^2;

% cs involves the third derivative of the profile; the fine grid is dense
% enough that a central difference is accurate to ~1e-8 relative
csp = zeros(1, M);
csp(2:M - 1) = (cs(3:M) - cs(1:M - 2)) ./ (sarc(3:M) - sarc(1:M - 2));
csp(1) = csp(2);
csp(M) = csp(M - 1);

% ---- thickness: constant through the bell, linear taper across the bow ----
taper = min(1, max(0, (a - r) / den));
thk = P.thickness_rim_m + (P.thickness_bell_m - P.thickness_rim_m) * taper;

% ---- nodes, equally spaced in arclength ----
nn = max(6, round(P.n_elements)) + 1;
snode = linspace(0, S, nn);

G.r = r;
G.z = z;
G.s = sarc;
G.thk = thk;
G.drds = drds;
G.dzds = dzds;
G.cs = cs;
G.ct = ct;
G.csp = csp;
G.ctp = ctp;
G.S = S;
G.M = M;
G.nn = nn;
G.snode = snode;
G.rnode = interp1(sarc, r, snode);
G.znode = interp1(sarc, z, snode);
G.drdsnode = interp1(sarc, drds, snode);
G.dzdsnode = interp1(sarc, dzds, snode);
G.thknode = interp1(sarc, thk, snode);
G.area = 2 * pi * trapz(sarc, r);
G.mass = 2 * pi * P.density * trapz(sarc, r .* thk);
end