1function [acc, dinfo, fhist] = cymbal2_sim(P, S, MB)
2% CYMBAL2_SIM Time-domain integration of the full system: von Karman
3% shallow spherical shell (implicit in the stiff linear terms, explicit in
4% the nonlinearity), washer support built into the boundary rows, mallet
5% head + stick coupled through the Hertz contact. Returns the normal
6% acceleration of the shell at the pickup node at the simulation rate,
7% diagnostics, and the contact force history.
9G = S.G;
10nrt = G.nrt;
11Nr = G.Nr;
12nth = G.nth;
13ht = S.ht;
14rhoh = S.rhoh;
15Rs = S.Rs;
16EhRs = S.Eh / Rs;
17Eh2 = S.Eh / 2;
18ls = S.ls; js = S.js; lp = S.lp; jp = S.jp;
19fsc = S.forceScale;
20nlr0 = S.eomRow0;
22K = max(1, round(P.strike_count));
24% Stroke schedule; per-stroke peak hand speeds with a seeded spread. Each
25% stroke is launched early by its predicted hand travel time (acceleration
26% plus coast over the rest clearance), so the contacts land near the
27% nominal beats instead of drifting with velocity. A nonpositive
28% strike_gap_m asks for the clearance to be sized to the stroke: contact
29% arrives just after the hand reaches full speed, which keeps the tip's
30% whip phase, and with it the contact speed, consistent across strokes.
31tstr = P.pre_delay_ms * 1e-3 + (0:K - 1) / max(1e-3, P.strike_rate_hz);
32vstr = zeros(1, K);
33state = P.seed * 2654435761;
34for i = 1:K
35 state = mod(state * 1664525 + 1013904223, 4294967296);
36 x = state / 4294967296;
37 vstr(i) = max(0.05, P.strike_velocity * (1 + P.velocity_spread * (2 * x - 1)));
38end
39gap = P.strike_gap_m;
40if gap <= 0
41 gap = 1.05 * max(vstr)^2 / (2 * P.strike_accel) + 0.005;
42end
43travel = gap ./ vstr + vstr / (2 * P.strike_accel);
44tlaunch = tstr - travel;
45off = -min(0, min(tlaunch));
46tlaunch = tlaunch + off;
48duration = off + (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;
49nsteps = floor(duration * P.sample_rate * P.oversample);
51KM = P.head_stiffness;
52alM = P.head_alpha;
53tipi = MB.tipi;
54ps = MB.ps;
56% Hand state machine: the hand accelerates downward, coasts at the stroke
57% speed until the head makes contact (or a miss depth), holds briefly,
58% then retracts at half speed and brakes so it arrives back at rest,
59% ready for the next scheduled stroke. This replaces the paper's
60% constant-acceleration hand, which never reverses. The machine runs one
61% step ahead because the beam coupling needs u_H at k + 1.
62aH = P.strike_accel;
63aHr = P.retract_accel;
64holdSteps = max(1, round(P.strike_hold_ms * 1e-3 / ht));
65missDepth = gap + 0.06;
66uHm = 0; uH0 = 0; uHp = 0;
67vH = 0;
68phase = 0;
69strokeIdx = 1;
70holdLeft = 0;
71vtar = 0;
73D1p = G.D1; D2p = G.D2;
74T1t = G.T1t; T2 = G.T2;
75F = G.F; Gi = G.Gi;
76Rd = G.Rd; R2d = G.R2d;
77AiB = S.AiB; AiC = S.AiC; Ai = S.Ai; Fi = S.Fi;
78sNL = ht^2 / rhoh;
80U = zeros(nrt, nth);
81U1 = zeros(nrt, nth);
82Uh3 = zeros(nrt, 1, nth);
83Uh13 = zeros(nrt, 1, nth);
84d = zeros(MB.nfree, 1);
85d1 = zeros(MB.nfree, 1);
87acc = zeros(1, nsteps);
88fhist = zeros(1, nsteps);
89shellActive = 0;
90peakDisp = 0;
91peakXi = 0;
92inContact = 0;
93cvel = zeros(1, K);
95for k = 1:nsteps
96 % ---- mallet: contact force from the state at step k ----
97 xi = U(ls, js) - d(tipi) - gap;
98 if xi > 0
99 fM1 = KM * xi^alM;
100 shellActive = 1;
101 if inContact == 0 && strokeIdx <= K && cvel(strokeIdx) == 0
102 cvel(strokeIdx) = (d(tipi) - d1(tipi)) / ht; % tip speed at contact
103 end
104 inContact = 1;
105 if xi > peakXi
106 peakXi = xi;
107 end
108 else
109 fM1 = 0;
110 inContact = 0;
111 end
112 fhist(k) = fM1;
113 hA = -(uHp - 2 * uH0 + uHm);
114 hB = -(ht / 2) * (uHp - uHm);
115 hC = -ht^2 * (ps * uH0 + ((1 - ps) / 2) * (uHp + uHm));
116 dnew = MB.ASiB * d + MB.ASiC * d1 + MB.vM * hA + MB.vR * hB + MB.vK * hC ...
117 + (ht^2 * fM1) * MB.tipcol;
119 % ---- advance the hand one step ----
120 if phase == 0
121 if strokeIdx <= K && k * ht >= tlaunch(strokeIdx)
122 phase = 1;
123 vtar = vstr(strokeIdx);
124 end
125 elseif phase == 1
126 vH = max(vH - aH * ht, -vtar);
127 if inContact ~= 0
128 phase = 3;
129 holdLeft = holdSteps;
130 elseif vH <= -vtar
131 phase = 2;
132 end
133 elseif phase == 2
134 if inContact ~= 0
135 phase = 3;
136 holdLeft = holdSteps;
137 elseif uHp <= -missDepth
138 phase = 4;
139 end
140 elseif phase == 3
141 holdLeft = holdLeft - 1;
142 if holdLeft <= 0
143 phase = 4;
144 end
145 elseif phase == 4
146 vH = min(vH + aHr * ht, 0.5 * vtar);
147 if vH > 0 && uHp >= -vH^2 / (2 * aHr)
148 phase = 5;
149 end
150 else
151 vH = vH - aHr * ht;
152 if vH <= 0 || uHp >= 0
153 vH = 0;
154 phase = 0;
155 strokeIdx = strokeIdx + 1;
156 end
157 end
158 uHm = uH0;
159 uH0 = uHp;
160 uHp = uHp + vH * ht;
162 if shellActive ~= 0
163 % ---- shell derivatives and von Karman brackets at step k ----
164 Ur = D1p * U;
165 Urr = D2p * U;
166 v1 = Rd .* Ur + R2d .* (U * T2);
167 v2 = Rd .* (Ur * T1t) - R2d .* (U * T1t);
168 LapU = Urr + v1;
169 rhsA = EhRs * LapU - Eh2 * 2 * (Urr .* v1 - v2 .* v2);
170 rhsA(1:3, :) = 0;
171 rhsA(Nr:nrt, :) = 0;
172 Ph3 = pagemtimes(Fi, reshape(rhsA * F, nrt, 1, nth));
173 Phi = reshape(Ph3, nrt, nth) * Gi;
174 Pr = D1p * Phi;
175 Prr = D2p * Phi;
176 w1 = Rd .* Pr + R2d .* (Phi * T2);
177 w2 = Rd .* (Pr * T1t) - R2d .* (Phi * T1t);
178 NL = Prr .* v1 + Urr .* w1 - 2 * w2 .* v2 - (Prr + w1) / Rs;
179 NL(ls, js) = NL(ls, js) - fM1 * fsc;
180 NL(1:nlr0 - 1, :) = 0;
181 NL(Nr + 1:nrt, :) = 0;
182 NLh = pagemtimes(Ai, reshape((sNL * NL) * F, nrt, 1, nth));
183 Uh3n = pagemtimes(AiB, Uh3) + pagemtimes(AiC, Uh13) + NLh;
184 Unew = reshape(Uh3n, nrt, nth) * Gi;
186 acc(k) = (Unew(lp, jp) - 2 * U(lp, jp) + U1(lp, jp)) / ht^2;
187 U1 = U;
188 U = Unew;
189 Uh13 = Uh3;
190 Uh3 = Uh3n;
191 if mod(k, 512) == 0
192 pd = max(abs(U(:, js)));
193 if pd > peakDisp
194 peakDisp = pd;
195 end
196 if ~isfinite(pd)
197 error('cymbal2: simulation diverged at t = %.4f s (reduce strike_velocity, or raise n_theta / oversample)', k * ht);
198 end
199 end
200 end
202 d1 = d;
203 d = dnew;
204end
206% Contact bookkeeping from the force history
207inC = fhist > 0;
208dinfo.peakForceN = max(fhist);
209dinfo.contactSteps = sum(inC);
210dinfo.contactMs = 1000 * ht * dinfo.contactSteps;
211dinfo.contactEvents = sum(diff([0, inC]) > 0);
212dinfo.peakXiMm = peakXi * 1000;
213dinfo.peakDispM = peakDisp;
214dinfo.peakDispThickness = peakDisp / P.thickness_m;
215dinfo.peakAccel = max(abs(acc));
216dinfo.nsteps = nsteps;
217dinfo.simRate = P.sample_rate * P.oversample;
218dinfo.tstr = tstr;
219dinfo.vstr = vstr;
220dinfo.meanV = mean(vstr);
221dinfo.contactVel = cvel; % tip speed at first contact of each stroke
222dinfo.strokesPlayed = min(strokeIdx - 1 + (phase > 0), K);
223end