// Predefined design scripts, selectable from the interface. Each is a
// MATLAB-syntax script that fills in the `design` struct; see design.ts for
// the fields. A script is self-contained: the shape functions it uses are
// ordinary MATLAB functions written at the bottom of the script itself, so
// that reading the script tells the whole story of the design and editing it
// cannot go wrong for want of a hidden helper.

export interface Preset {
  id: string;
  label: string;
  source: string;
}

const flatRectangle = `% A flat rectangular sheet of stainless steel, hanging freely (no mount).
% The script is an ordinary MATLAB program that fills in a struct named
% design. Lengths are meters; the constants at the top are the only units.
% A field can be a number or a function handle @(x, y, r, theta) of the
% point, measured from the center of the bounding box; the handles are
% applied to whole grid matrices, so use elementwise operators (.*, ./, .^)
% inside them. Anything else a script needs it defines itself, as MATLAB
% functions at the bottom; the other presets show how.
inch = 0.0254;                  % m
mm = 1e-3;                      % m

design.name = 'Flat steel rectangle';
design.material.E = 170e9;      % Young's modulus [Pa]
design.material.nu = 0.29;      % Poisson ratio
design.material.rho = 7900;     % density [kg/m^3]
design.bounds.width = 12*inch;  % bounding box [m]
design.bounds.height = 18*inch;
design.thickness = 0.79*mm;
design.profile = 0;             % stress-free height [m]: flat
design.mount = [];              % free plate
design.loss.sigma0 = 1.15;      % frequency-independent loss [1/s]
design.loss.sigma1 = 5e-4;      % loss proportional to k^2 [m^2/s]
design.loss.sigmaIn = 20;       % extra loss of the in-plane motion [1/s]
`;

const bowedRectangle = `% The same steel sheet with a 3 mm parabolic bow across its width, the
% kind of bow that rolling and handling leave in real sheet metal. The bow
% makes the sheet a shallow shell, which oil-cans and wobbles when struck.
inch = 0.0254; mm = 1e-3;
W = 12*inch; H = 18*inch;

design.name = 'Bowed steel rectangle';
design.material.E = 170e9;
design.material.nu = 0.29;
design.material.rho = 7900;
design.bounds.width = W;
design.bounds.height = H;
design.thickness = 0.79*mm;
design.profile = @(x, y, r, theta) bow(x, W, 3*mm, 0.3);
design.mount = [];
design.loss.sigma0 = 1.15;
design.loss.sigma1 = 5e-4;
design.loss.sigmaIn = 20;

function v = bow(s, span, sag, offset)
% A parabola of the given sag across a span centered on 0; s runs along the
% span, and the offset in -1..1 moves the crest toward one end.
t = 2*s/span;
v = sag * (1 - t.^2) .* (1 + offset*t);
end
`;

const flatDisc = `% A flat 16 inch steel disc with a center hole, sitting on a stand: the
% nodes within the mount radius are tied to the ground by a spring (the felt)
% and a dashpot, spread over the mount area.
inch = 0.0254; mm = 1e-3;
D = 16*inch; hole = 12*mm;

design.name = 'Flat steel disc';
design.material.E = 170e9;
design.material.nu = 0.29;
design.material.rho = 7900;
design.bounds.width = D;
design.bounds.height = D;
design.inside = @(x, y, r, theta) (r <= D/2) & (r >= hole/2);
design.thickness = 1.0*mm;
design.profile = 0;
design.mount.radius = 20*mm;        % felt radius
design.mount.stiffness = 2e4;       % [N/m]
design.mount.damping = 5;           % [N s/m]
design.loss.sigma0 = 1.15;
design.loss.sigma1 = 5e-4;
design.loss.sigmaIn = 20;
`;

const crash16 = `% A 16 inch crash: B20 bronze, tapered from the bell to the edge, with a
% shallow bow and a cosine bell. The shallow-shell model behind the
% simulation is strained where the slope of the profile exceeds about 0.1
% to 0.2, so the bow and the bell are a little flatter than on a real crash.
inch = 0.0254; mm = 1e-3;
D = 16*inch; R = D/2; hole = 12*mm;
bellR = 60*mm; bellH = 12*mm;   % base radius and height of the bell
bowH = 16*mm;                   % rise from the edge to the foot of the bell
tBell = 1.4*mm; tEdge = 0.8*mm; % thickness at the bell and at the edge

design.name = '16 inch crash';
design.material.E = 110e9;
design.material.nu = 0.34;
design.material.rho = 8800;
design.bounds.width = D;
design.bounds.height = D;
design.inside = @(x, y, r, theta) (r <= R) & (r >= hole/2);
% Linear taper from the foot of the bell to the edge.
design.thickness = @(x, y, r, theta) tBell + (tEdge - tBell) * ramp(bellR, R, r);
% A paraboloidal bow plus the bell; dome has zero slope at its rim.
design.profile = @(x, y, r, theta) bowH * (1 - (r/R).^2) + dome(r, bellR, bellH);
design.mount.radius = 20*mm;
design.mount.stiffness = 2e4;
design.mount.damping = 5;
design.loss.sigma0 = 0.6;
design.loss.sigma1 = 2.5e-4;
design.loss.sigmaIn = 20;

function v = ramp(a, b, x)
% 0 at x = a, 1 at x = b, clamped to that range outside.
v = min(max((x - a) ./ (b - a), 0), 1);
end

function v = dome(d, radius, height)
% A cosine dome of the given height and base radius, with zero slope both at
% its rim and at its apex; d is the distance from its center.
v = 0.5 * height * (1 + cos(pi * d / radius)) .* (d < radius);
end
`;

const crashDented = `% A 16 inch crash that has been through a lot: the same bronze, taper,
% bow, and bell as the clean crash, plus a dozen dents scattered over the
% bow and a faint hammering texture. The dents are Gaussian bumps of a few
% centimeters (broad compared with the 8 mm grid) and a millimeter or two
% in height, up or down; change the seed for a different set. Dents break
% the symmetry of the disc, split the degenerate mode pairs, and give the
% cymbal a less pure, more trashy tone.
inch = 0.0254; mm = 1e-3;
D = 16*inch; R = D/2; hole = 12*mm;
bellR = 60*mm; bellH = 12*mm; bowH = 16*mm;
tBell = 1.4*mm; tEdge = 0.8*mm;

% Twelve dents from the seeded generator at the bottom (the same on every run).
nd = 12;
u = reshape(seededrand(7, 4*nd), 4, nd);            % rows: radius, angle, size, height
rr = bellR + 25*mm + (R - bellR - 45*mm) * u(1,:);  % on the bow, clear of the bell and the edge
th = 2*pi * u(2,:);
cx = rr .* cos(th); cy = rr .* sin(th);
rad = (18 + 14*u(3,:)) * mm;
hgt = sign(u(4,:) - 0.5) .* (0.8 + 1.4*u(4,:)) * mm;
hammer = noise2d(3, 40*mm);                         % unit RMS, correlation length 40 mm

design.name = '16 inch crash, dented';
design.material.E = 110e9;
design.material.nu = 0.34;
design.material.rho = 8800;
design.bounds.width = D;
design.bounds.height = D;
design.inside = @(x, y, r, theta) (r <= R) & (r >= hole/2);
design.thickness = @(x, y, r, theta) tBell + (tEdge - tBell) * ramp(bellR, R, r);
design.profile = @(x, y, r, theta) bowH * (1 - (r/R).^2) + dome(r, bellR, bellH) ...
  + dents(x, y, cx, cy, rad, hgt) ...
  + 0.3*mm * hammer(x, y) .* smoothstep(bellR, bellR + 30*mm, r);   % no hammering on the bell
design.mount.radius = 20*mm;
design.mount.stiffness = 2e4;
design.mount.damping = 5;
design.loss.sigma0 = 0.6;
design.loss.sigma1 = 2.5e-4;
design.loss.sigmaIn = 20;

function v = ramp(a, b, x)
% 0 at x = a, 1 at x = b, clamped to that range outside.
v = min(max((x - a) ./ (b - a), 0), 1);
end

function v = smoothstep(e0, e1, x)
% 0 below e0, 1 above e1, a smooth cubic between, with zero slope at both.
t = ramp(e0, e1, x);
v = t.^2 .* (3 - 2*t);
end

function v = dome(d, radius, height)
% A cosine dome of the given height and base radius, with zero slope both at
% its rim and at its apex; d is the distance from its center.
v = 0.5 * height * (1 + cos(pi * d / radius)) .* (d < radius);
end

function v = dents(x, y, cx, cy, radius, height)
% The sum of Gaussian bumps centered at (cx(k), cy(k)), each of the given
% height and falling to 1/e of it at the given radius.
v = zeros(size(x));
for k = 1:numel(cx)
  d2 = (x - cx(k)).^2 + (y - cy(k)).^2;
  v = v + height(k) * exp(-d2 / radius(k)^2);
end
end

function u = seededrand(seed, n)
% n uniform numbers in [0, 1) from a small seeded generator, the same on
% every run and every machine.
s = floor(seed);
u = zeros(1, n);
for k = 1:n
  s = mod(s * 1664525 + 1013904223, 4294967296);
  u(k) = s / 4294967296;
end
end

function f = noise2d(seed, len)
% A smooth random field f(x, y) of unit RMS with the covariance
% exp(-d^2 / len^2), the same for a given seed on any grid: a sum of plane
% waves whose wavevectors are drawn from the matching Gaussian spectrum.
n = 64;
u = reshape(seededrand(seed, 5*n), 5, n);
kx = sqrt(-2*log(1 - u(1,:))) .* cos(2*pi*u(2,:)) * sqrt(2) / len;
ky = sqrt(-2*log(1 - u(3,:))) .* cos(2*pi*u(4,:)) * sqrt(2) / len;
ph = 2*pi*u(5,:);
f = @(x, y) noisesum(kx, ky, ph, x, y);
end

function v = noisesum(kx, ky, ph, x, y)
v = zeros(size(x));
for k = 1:numel(kx)
  v = v + cos(kx(k) * x + ky(k) * y + ph(k));
end
v = v / sqrt(numel(kx) / 2);
end
`;

const ride20 = `% A 20 inch ride: heavier bronze, a larger bell, a flatter bow. The taper
% is stronger than on the crash, and the mount is a little stiffer.
inch = 0.0254; mm = 1e-3;
D = 20*inch; R = D/2; hole = 12*mm;
bellR = 70*mm; bellH = 14*mm;
bowH = 15*mm;
tBell = 2.0*mm; tEdge = 1.1*mm;

design.name = '20 inch ride';
design.material.E = 110e9;
design.material.nu = 0.34;
design.material.rho = 8800;
design.bounds.width = D;
design.bounds.height = D;
design.inside = @(x, y, r, theta) (r <= R) & (r >= hole/2);
design.thickness = @(x, y, r, theta) tBell + (tEdge - tBell) * smoothstep(bellR, R, r);
design.profile = @(x, y, r, theta) bowH * (1 - (r/R).^2) + dome(r, bellR, bellH);
design.mount.radius = 22*mm;
design.mount.stiffness = 3e4;
design.mount.damping = 8;
design.loss.sigma0 = 0.5;
design.loss.sigma1 = 2e-4;
design.loss.sigmaIn = 20;

function v = smoothstep(e0, e1, x)
% 0 below e0, 1 above e1, a smooth cubic between, with zero slope at both.
t = min(max((x - e0) ./ (e1 - e0), 0), 1);
v = t.^2 .* (3 - 2*t);
end

function v = dome(d, radius, height)
% A cosine dome of the given height and base radius, with zero slope both at
% its rim and at its apex; d is the distance from its center.
v = 0.5 * height * (1 + cos(pi * d / radius)) .* (d < radius);
end
`;

const ellipseTwist = `% A steel ellipse with a saddle (twist) and a bow, hanging freely. The
% twist raises two opposite quadrants and lowers the other two, which
% breaks the symmetry of the disc modes.
inch = 0.0254; mm = 1e-3;
W = 14*inch; H = 20*inch;

design.name = 'Steel ellipse with a twist';
design.material.E = 170e9;
design.material.nu = 0.29;
design.material.rho = 7900;
design.bounds.width = W;
design.bounds.height = H;
% 2x/W and 2y/H run from -1 to 1 over the bounding box.
design.inside = @(x, y, r, theta) (2*x/W).^2 + (2*y/H).^2 <= 1;
design.thickness = 0.79*mm;
design.profile = @(x, y, r, theta) 2*mm * (2*x/W) .* (2*y/H) + bow(y, H, 2*mm);
design.mount = [];
design.loss.sigma0 = 1.15;
design.loss.sigma1 = 5e-4;
design.loss.sigmaIn = 20;

function v = bow(s, span, sag)
% A parabola of the given sag across a span centered on 0.
v = sag * (1 - (2*s/span).^2);
end
`;

const china18 = `% An 18 inch china: a conical bow that falls from the bell toward a flange,
% and an upturned edge that rises again over the outer 45 mm. The bell is
% square-topped (a plateau with rounded shoulders). The edge slope is well
% beyond the shallow-shell range, so this preset is the least trustworthy.
inch = 0.0254; mm = 1e-3;
D = 18*inch; R = D/2; hole = 12*mm;
bellR = 65*mm; bellH = 14*mm;
flangeR = R - 45*mm; bowH = 14*mm; edgeH = 12*mm;
tBell = 1.6*mm; tEdge = 0.9*mm;

design.name = '18 inch china';
design.material.E = 110e9;
design.material.nu = 0.34;
design.material.rho = 8800;
design.bounds.width = D;
design.bounds.height = D;
design.inside = @(x, y, r, theta) (r <= R) & (r >= hole/2);
design.thickness = @(x, y, r, theta) tBell + (tEdge - tBell) * smoothstep(bellR, flangeR, r);
design.profile = @(x, y, r, theta) bowH * max(0, 1 - r/flangeR) ...        % cone down to the flange
  + edgeH * smoothstep(flangeR, R, r) ...                                   % upturned edge
  + bellH * (1 - smoothstep(bellR - 20*mm, bellR, r));                      % flat-topped bell
design.mount.radius = 20*mm;
design.mount.stiffness = 2e4;
design.mount.damping = 5;
design.loss.sigma0 = 0.7;
design.loss.sigma1 = 3e-4;
design.loss.sigmaIn = 20;

function v = smoothstep(e0, e1, x)
% 0 below e0, 1 above e1, a smooth cubic between, with zero slope at both.
t = min(max((x - e0) ./ (e1 - e0), 0), 1);
v = t.^2 .* (3 - 2*t);
end
`;

export const PRESETS: Preset[] = [
  { id: "flat-rectangle", label: "Flat steel rectangle", source: flatRectangle },
  { id: "bowed-rectangle", label: "Bowed steel rectangle", source: bowedRectangle },
  { id: "flat-disc", label: "Flat steel disc on a stand", source: flatDisc },
  { id: "crash-16", label: "16 inch crash", source: crash16 },
  { id: "crash-16-dented", label: "16 inch crash, dented", source: crashDented },
  { id: "ride-20", label: "20 inch ride", source: ride20 },
  { id: "ellipse-twist", label: "Steel ellipse with a twist", source: ellipseTwist },
  { id: "china-18", label: "18 inch china", source: china18 },
];

export const DEFAULT_PRESET = PRESETS[0];