function [uh, wh, thk, drds, dzds, cs, ct] = cymbal_at(basis, sfrac)
% CYMBAL_AT Mode shapes and local geometry at an arclength fraction.
%
% sfrac addresses the meridian from the hole (0) to the rim (1); it is the
% coordinate the strike positions are expressed in. The meridional shape U is
% linear within an element, so it is interpolated linearly; the normal shape W
% is cubic Hermite, and is interpolated with its own shape functions from the
% nodal (W, dW/ds) pair, which matters because at the top of the retained band
% there are only about four nodes per flexural wavelength.
%
% Returns row vectors uh, wh (1 x n modes) and the local thickness, meridian
% direction cosines and principal curvatures at that point.

nn = length(basis.snode);
S = basis.S;
s = min(0.999999, max(0, sfrac)) * S;
idx = sum(basis.snode <= s);
idx = min(max(idx, 1), nn - 1);
L = basis.snode(idx + 1) - basis.snode(idx);
xi = (s - basis.snode(idx)) / L;

uh = basis.NU(idx, :) + xi * (basis.NU(idx + 1, :) - basis.NU(idx, :));
h1 = 1 - 3 * xi^2 + 2 * xi^3;
h2 = L * (xi - 2 * xi^2 + xi^3);
h3 = 3 * xi^2 - 2 * xi^3;
h4 = L * (-xi^2 + xi^3);
wh = h1 * basis.NW(idx, :) + h2 * basis.NT(idx, :) ...
   + h3 * basis.NW(idx + 1, :) + h4 * basis.NT(idx + 1, :);

thk = basis.thknode(idx) + xi * (basis.thknode(idx + 1) - basis.thknode(idx));
drds = basis.drdsnode(idx) + xi * (basis.drdsnode(idx + 1) - basis.drdsnode(idx));
dzds = basis.dzdsnode(idx) + xi * (basis.dzdsnode(idx + 1) - basis.dzdsnode(idx));
cs = basis.csnode(idx) + xi * (basis.csnode(idx + 1) - basis.csnode(idx));
ct = basis.ctnode(idx) + xi * (basis.ctnode(idx + 1) - basis.ctnode(idx));
end