1"""Cymbal mid-surface height z0(x, y) and thickness h(x, y).
3The shell is described on the flat projected domain, the disk of radius R. The mid-surface height
4z0 is the sum of a bow, a bell, a low-order rim warp and a set of hammering dents; the thickness
5tapers linearly along the bow and steps up smoothly into the bell. All pieces are written with
6jax.numpy so that the slope grad z0, which enters the membrane strains of the shallow shell, is
7obtained by automatic differentiation rather than by finite differences. Imperfections are drawn
8from a seeded generator so that a run is reproducible. See docs/model.md, section "Geometry".
9"""
10from __future__ import annotations
12import math
14import jax
15import jax.numpy as jnp
16import numpy as np
18from .config import Imperfections, Profile
21class CymbalProfile:
22 """Callable description of the cymbal surface, including the seeded imperfections.
24 ``z0(x, y)`` and ``thickness(x, y)`` take scalars and are differentiable; ``z0_at``,
25 ``grad_z0_at`` and ``thickness_at`` take an (n, 2) array of points.
26 """
28 def __init__(self, profile: Profile, imperfections: Imperfections):
29 self.profile = profile
30 self.imperfections = imperfections
31 p, im = profile, imperfections
32 rng = np.random.default_rng(im.seed)
34 # hammer dents: uniform over the annulus between the shoulder and the rim, random sign
35 r0, r1 = p.bell_shoulder_radius + im.dent_margin, p.radius - im.dent_margin
36 n = int(im.n_dents)
37 rr = np.sqrt(rng.uniform(r0**2, r1**2, n))
38 tt = rng.uniform(0.0, 2 * np.pi, n)
39 self.dent_x = jnp.asarray(rr * np.cos(tt))
40 self.dent_y = jnp.asarray(rr * np.sin(tt))
41 self.dent_amp = jnp.asarray(im.dent_depth * rng.choice([-1.0, 1.0], n) * rng.uniform(0.5, 1.5, n))
42 self.dent_sigma = jnp.asarray(im.dent_sigma * rng.uniform(0.7, 1.3, n))
44 # rim warp: azimuthal harmonics growing linearly from the shoulder to the rim
45 modes = tuple(int(m) for m in im.rim_warp_modes)
46 self.warp_m = jnp.asarray(modes, dtype=float)
47 self.warp_a = jnp.asarray([im.rim_warp * 2.0 / m for m in modes], dtype=float)
48 self.warp_phase = jnp.asarray(rng.uniform(0.0, 2 * np.pi, len(modes)))
50 # shoulder ovality and bell offset
51 self.oval_phase = float(rng.uniform(0.0, 2 * np.pi))
52 bo = float(rng.uniform(0.0, 2 * np.pi))
53 self.bell_cx = im.bell_offset * math.cos(bo)
54 self.bell_cy = im.bell_offset * math.sin(bo)
56 self._z0_v = jax.vmap(self.z0)
57 self._h_v = jax.vmap(self.thickness)
58 self._g_v = jax.vmap(jax.grad(lambda pt: self.z0(pt[0], pt[1])))
60 # ---- building blocks
62 def _polar(x, y):
63 r = jnp.sqrt(x * x + y * y + 1e-30)
64 phi = jnp.arctan2(y, jnp.where(r < 1e-6, 1e-6, x)) # finite gradient at the centre
65 return r, phi
67 def _shoulder_radius(self, phi):
68 p, im = self.profile, self.imperfections
69 return p.bell_shoulder_radius * (1.0 + im.shoulder_ovality * jnp.cos(2 * phi + self.oval_phase))
71 def _bell_blend(self, r_bell, r_sh):
72 """1 on the bell plateau, 0 beyond the shoulder, C1 cosine ramp in between."""
73 t = jnp.clip((r_bell - self.profile.bell_top_radius) / (r_sh - self.profile.bell_top_radius), 0.0, 1.0)
74 return 0.5 * (1.0 + jnp.cos(jnp.pi * t))
76 def _bow_coordinate(self, x, y):
77 """s in [0, 1]: 0 at the rim, 1 at the shoulder; plus the bell radius and shoulder radius."""
78 p = self.profile
79 r, phi = self._polar(x, y)
80 r_sh = self._shoulder_radius(phi)
81 r_bell = jnp.sqrt((x - self.bell_cx) ** 2 + (y - self.bell_cy) ** 2 + 1e-30)
82 s = jnp.clip((p.radius - r) / (p.radius - r_sh), 0.0, 1.0)
83 return s, phi, r_bell, r_sh
85 # ---- public scalar functions
86 def z0(self, x, y):
87 """Mid-surface height z0(x, y) [m]."""
88 p = self.profile
89 s, phi, r_bell, r_sh = self._bow_coordinate(x, y)
90 u = jnp.maximum(1.0 - s, 1e-12) # 1 at the rim, 0 at the shoulder
91 bow = p.bow_height * (1.0 - u ** p.bow_power)
92 bell = p.bell_height * self._bell_blend(r_bell, r_sh)
93 warp = u * jnp.sum(self.warp_a * jnp.cos(self.warp_m * phi + self.warp_phase))
94 dents = jnp.sum(self.dent_amp * jnp.exp(
95 -((x - self.dent_x) ** 2 + (y - self.dent_y) ** 2) / (2 * self.dent_sigma ** 2)))
96 return bow + bell + warp + dents
98 def thickness(self, x, y):
99 """Thickness h(x, y) [m]: linear taper along the bow, smooth step up into the bell."""
100 p = self.profile
101 s, phi, r_bell, r_sh = self._bow_coordinate(x, y)
102 bow = p.h_rim + (p.h_shoulder - p.h_rim) * s
103 return bow + (p.h_bell - p.h_shoulder) * self._bell_blend(r_bell, r_sh)
105 # ---- vectorized evaluation on (n, 2) point arrays
106 def z0_at(self, xy) -> np.ndarray:
107 xy = jnp.asarray(xy, dtype=jnp.float64)
108 return np.asarray(self._z0_v(xy[:, 0], xy[:, 1]))
110 def thickness_at(self, xy) -> np.ndarray:
111 xy = jnp.asarray(xy, dtype=jnp.float64)
112 return np.asarray(self._h_v(xy[:, 0], xy[:, 1]))
114 def grad_z0_at(self, xy) -> np.ndarray:
115 """Analytic slope (dz0/dx, dz0/dy) at each point, shape (n, 2)."""
116 return np.asarray(self._g_v(jnp.asarray(xy, dtype=jnp.float64)))