fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 1// Time integration of the bead system.
2//
3// The conservative part (springs, external force pulses) is integrated with
4// velocity Verlet, optionally composed into Yoshida's fourth-order scheme.
5// Damping (friction and dashpots) is linear in the velocities, dv/dt = -A v,
6// and is applied as a separate exact step v <- exp(-tau A) v, once per audio
7// sample (Strang splitting over the sample interval Delta = 1/fs). The
8// splitting error is second order in Delta and scales with the damping rate,
9// which is small compared with the oscillation frequencies; the conservative
10// part is second order (Verlet) or fourth order (Yoshida) in the step h.
a8b5c65Dome example: chaos that keeps its pitch. Constant load with static relaxation (FIRE), full 3D linear modes for curved nets, spectrum-based default frequency range, screening scriptJeremy Magland 12import { loadForces, springBlocks, springForces } from "./forces";
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 13import type { BeadSystem, Method, Modes, SimResult } from "./types";
15const YOSHIDA_W1 = 1 / (2 - Math.cbrt(2));
16const YOSHIDA_W0 = -Math.cbrt(2) * YOSHIDA_W1;
18export function stepCoefficients(method: Method): number[] {
19 return method === "yoshida4" ? [YOSHIDA_W1, YOSHIDA_W0, YOSHIDA_W1] : [1];
20}
22// Largest h*omega for which the method is stable on x'' = -omega^2 x.
23export function stabilityLimit(method: Method): number {
24 if (method === "verlet") return 2;
25 const cs = stepCoefficients(method);
26 const stable = (h: number) => {
27 // transfer matrix of one full step for omega = 1
28 let m = [1, 0, 0, 1]; // [x<-x, x<-v, v<-x, v<-v]
29 const mul = (a: number[], b: number[]) => [
30 a[0] * b[0] + a[1] * b[2],
31 a[0] * b[1] + a[1] * b[3],
32 a[2] * b[0] + a[3] * b[2],
33 a[2] * b[1] + a[3] * b[3],
34 ];
35 for (const c of cs) {
36 const kick = [1, 0, (-c * h) / 2, 1];
37 const drift = [1, c * h, 0, 1];
38 m = mul(kick, mul(drift, mul(kick, m)));
39 }
40 return Math.abs(m[0] + m[3]) < 2;
41 };
42 let h = 1e-3;
43 while (stable(h + 1e-3)) h += 1e-3;
44 return h;
45}
47// Error of the numerical oscillation frequency for x'' = -omega^2 x, in cents
48// (positive: the simulated tone is sharp).
49export function frequencyErrorCents(method: Method, omega: number, h: number): number {
50 let m = [1, 0, 0, 1];
51 const mul = (a: number[], b: number[]) => [
52 a[0] * b[0] + a[1] * b[2],
53 a[0] * b[1] + a[1] * b[3],
54 a[2] * b[0] + a[3] * b[2],
55 a[2] * b[1] + a[3] * b[3],
56 ];
57 for (const c of stepCoefficients(method)) {
58 const kick = [1, 0, (-c * h * omega * omega) / 2, 1];
59 const drift = [1, c * h, 0, 1];
60 m = mul(kick, mul(drift, mul(kick, m)));
61 }
62 const half = (m[0] + m[3]) / 2;
63 if (Math.abs(half) >= 1) return NaN;
64 const omegaNum = Math.acos(half) / h;
65 return 1200 * Math.log2(omegaNum / omega);
66}
68export interface SimOptions {
69 modes?: Modes | null;
70 onProgress?: (t: number) => void;
95f2216Minimal chaos example and perturbed-twin chaos test with separation plot and growth-rate fitJeremy Magland 71 // called at every frame record with the full-precision state (do not keep the arrays)
72 onFrame?: (frame: number, t: number, x: Float64Array, v: Float64Array) => void;
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 73 maxAnimBytes?: number;
74}
76export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
77 const t0 = performance.now();
78 const warnings: string[] = [];
79 const { n, mass, pinned } = sys;
80 const N3 = 3 * n;
81 const x = new Float64Array(sys.pos0);
82 const v = new Float64Array(N3);
83 const F = new Float64Array(N3);
84 const invm = new Float64Array(n);
85 for (let i = 0; i < n; i++) invm[i] = pinned[i] ? 0 : 1 / mass[i];
87 // ---- forces ------------------------------------------------------------
88 const { si, sj, sk, sL0, sk3 } = sys;
89 const ns = si.length;
a8b5c65Dome example: chaos that keeps its pitch. Constant load with static relaxation (FIRE), full 3D linear modes for curved nets, spectrum-based default frequency range, screening scriptJeremy Magland 90 // spring forces plus the constant load; returns their potential energy
91 const internalForces = (): number => springForces(sys, x, F) + loadForces(sys, x, F);
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 92 const pulses = sys.hits.filter(h => h.duration > 0);
93 const externalForces = (t: number) => {
94 for (const p of pulses) {
95 const u = (t - p.t) / p.duration;
96 if (u <= 0 || u >= 1) continue;
97 const w = 0.5 * (1 - Math.cos(2 * Math.PI * u));
98 for (const b of p.beads) {
99 F[3 * b] += w * p.force[0];
100 F[3 * b + 1] += w * p.force[1];
101 F[3 * b + 2] += w * p.force[2];
102 }
103 }
104 };
106 // ---- rest state checks ---------------------------------------------------
107 let pe = internalForces();
108 const peRest = pe; // energy stored by pre-tension; energies are reported relative to rest
109 let residualForce = 0;
110 for (let i = 0; i < n; i++) {
111 if (pinned[i]) continue;
112 const f = Math.hypot(F[3 * i], F[3 * i + 1], F[3 * i + 2]);
113 if (f > residualForce) residualForce = f;
114 }
115 let maxTension = 0;
116 for (let s = 0; s < ns; s++) {
117 const i3 = 3 * si[s];
118 const j3 = 3 * sj[s];
119 const L = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
120 const e = L - sL0[s];
121 maxTension = Math.max(maxTension, Math.abs(sk[s] * e + sk3[s] * e * e * e));
122 }
123 if (residualForce > 1e-6 * Math.max(maxTension, 1e-30)) {
124 warnings.push(
125 `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
a8b5c65Dome example: chaos that keeps its pitch. Constant load with static relaxation (FIRE), full 3D linear modes for curved nets, spectrum-based default frequency range, screening scriptJeremy Magland 126 `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike. Set sys.relax = true to find the equilibrium first.`
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 127 );
128 }
130 // ---- time step -----------------------------------------------------------
131 const omegaMax = estimateOmegaMax(sys);
132 const fs = sys.sampleRate;
133 const hLimit = stabilityLimit(sys.method) / omegaMax;
134 let substeps: number;
135 if (sys.dt > 0) {
136 substeps = Math.max(1, Math.round(1 / (fs * sys.dt)));
137 if (1 / (fs * substeps) > hLimit) {
138 warnings.push(
139 `dt = ${(1 / (fs * substeps)).toExponential(3)} s exceeds the stability limit ${hLimit.toExponential(3)} s.`
140 );
141 }
142 } else {
143 substeps = Math.max(1, Math.ceil(1 / (fs * sys.safety * hLimit)));
144 }
145 const h = 1 / (fs * substeps);
146 const nSamples = Math.max(1, Math.round(sys.duration * fs));
147 const nSteps = nSamples * substeps;
148 const cs = stepCoefficients(sys.method);
150 // ---- damping: v <- exp(-tau A) v ------------------------------------------
151 const { di, dj, dc, friction } = sys;
152 const nd = di.length;
153 let normA = 0; // infinity-norm bound of A
154 {
155 const row = new Float64Array(n);
156 for (let i = 0; i < n; i++) row[i] = friction[i];
157 for (let d = 0; d < nd; d++) {
158 row[di[d]] += 2 * dc[d];
159 row[dj[d]] += 2 * dc[d];
160 }
161 for (let i = 0; i < n; i++) normA = Math.max(normA, row[i] * invm[i]);
162 }
163 const hasDamping = normA > 0;
0330c41Damping by decay times: bp_decay helper, per-mode decay-time plot, examples convertedJeremy Magland 164 const frictionOnly = dc.every(c => c === 0); // then exp(-tau A) is diagonal
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 165 const term = new Float64Array(N3);
166 const tmp = new Float64Array(N3);
167 const fa = new Float64Array(n);
168 for (let i = 0; i < n; i++) fa[i] = friction[i] * invm[i];
169 const dai = new Float64Array(nd);
170 const daj = new Float64Array(nd);
171 for (let d = 0; d < nd; d++) {
172 dai[d] = dc[d] * invm[di[d]];
173 daj[d] = dc[d] * invm[dj[d]];
174 }
175 const applyA = (w: Float64Array, out: Float64Array) => {
176 for (let i = 0; i < n; i++) {
177 const a = fa[i];
178 out[3 * i] = a * w[3 * i];
179 out[3 * i + 1] = a * w[3 * i + 1];
180 out[3 * i + 2] = a * w[3 * i + 2];
181 }
182 for (let d = 0; d < nd; d++) {
183 const i3 = 3 * di[d];
184 const j3 = 3 * dj[d];
185 const ai = dai[d];
186 const aj = daj[d];
187 const r0 = w[i3] - w[j3];
188 const r1 = w[i3 + 1] - w[j3 + 1];
189 const r2 = w[i3 + 2] - w[j3 + 2];
190 out[i3] += ai * r0;
191 out[i3 + 1] += ai * r1;
192 out[i3 + 2] += ai * r2;
193 out[j3] -= aj * r0;
194 out[j3 + 1] -= aj * r1;
195 out[j3 + 2] -= aj * r2;
196 }
197 };
198 const kinetic = () => {
199 let ke = 0;
200 for (let i = 0; i < n; i++) {
201 ke += 0.5 * mass[i] * (v[3 * i] ** 2 + v[3 * i + 1] ** 2 + v[3 * i + 2] ** 2);
202 }
203 return ke;
204 };
205 let dissipated = 0;
206 const applyDamping = (tau: number) => {
207 if (!hasDamping || tau === 0) return;
208 const keBefore = kinetic();
209 if (frictionOnly) {
210 for (let i = 0; i < n; i++) {
211 const f = Math.exp(-tau * friction[i] * invm[i]);
212 v[3 * i] *= f;
213 v[3 * i + 1] *= f;
214 v[3 * i + 2] *= f;
215 }
216 } else {
217 // Taylor series with scaling; |tau_s| * normA <= 0.5 so terms shrink fast.
218 const nsub = Math.max(1, Math.ceil((Math.abs(tau) * normA) / 0.5));
219 const ts = tau / nsub;
220 for (let r = 0; r < nsub; r++) {
221 let vmax = 0;
222 for (let q = 0; q < N3; q++) {
223 term[q] = v[q];
224 vmax = Math.max(vmax, Math.abs(v[q]));
225 }
226 if (vmax === 0) break;
227 for (let k = 1; k < 40; k++) {
228 applyA(term, tmp);
229 const c = -ts / k;
230 let tmax = 0;
231 for (let q = 0; q < N3; q++) {
232 const val = c * tmp[q];
233 term[q] = val;
234 v[q] += val;
235 const a = Math.abs(val);
236 if (a > tmax) tmax = a;
237 }
238 if (tmax <= 1e-17 * vmax) break;
239 }
240 }
241 }
242 dissipated += keBefore - kinetic();
243 };
245 // ---- recording -------------------------------------------------------------
246 const audio = new Float32Array(nSamples);
247 const lb = sys.listenBeads;
248 const [ldx, ldy, ldz] = sys.listenDir;
249 const pos0 = sys.pos0;
250 const listen = () => {
251 let s = 0;
252 for (const b of lb) {
253 s +=
254 (x[3 * b] - pos0[3 * b]) * ldx +
255 (x[3 * b + 1] - pos0[3 * b + 1]) * ldy +
256 (x[3 * b + 2] - pos0[3 * b + 2]) * ldz;
257 }
258 return s / lb.length;
259 };
260 const frameStride = Math.max(1, Math.round(fs / sys.frameRate));
261 const nFrames = Math.floor((nSamples - 1) / frameStride) + 1;
262 const maxAnimBytes = opts.maxAnimBytes ?? 150e6;
263 const animStride = frameStride * Math.max(1, Math.ceil((nFrames * N3 * 4) / maxAnimBytes));
264 const frameTimes = new Float64Array(nFrames);
265 const energy = new Float64Array(nFrames);
266 const dissipatedRec = new Float64Array(nFrames);
267 const animTimes: number[] = [];
268 const animFrames: Float32Array[] = [];
270 const modes = opts.modes ?? null;
271 const nModes = modes ? modes.freqs.length : 0;
272 const modeEnergy = modes ? new Float32Array(nFrames * nModes) : null;
273 let bz: Float64Array | null = null;
274 let bv: Float64Array | null = null;
275 let om2: Float64Array | null = null;
276 if (modes) {
a8b5c65Dome example: chaos that keeps its pitch. Constant load with static relaxation (FIRE), full 3D linear modes for curved nets, spectrum-based default frequency range, screening scriptJeremy Magland 277 bz = new Float64Array(modes.dof.length);
278 bv = new Float64Array(modes.dof.length);
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 279 om2 = new Float64Array(nModes);
280 for (let k = 0; k < nModes; k++) om2[k] = (2 * Math.PI * modes.freqs[k]) ** 2;
281 }
282 const recordModes = (f: number) => {
283 if (!modes || !modeEnergy || !bz || !bv || !om2) return;
a8b5c65Dome example: chaos that keeps its pitch. Constant load with static relaxation (FIRE), full 3D linear modes for curved nets, spectrum-based default frequency range, screening scriptJeremy Magland 284 const dof = modes.dof;
285 const nr = dof.length;
286 for (let r = 0; r < nr; r++) {
287 const q = dof[r];
288 const sm = Math.sqrt(mass[Math.floor(q / 3)]);
289 bz[r] = sm * (x[q] - pos0[q]);
290 bv[r] = sm * v[q];
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 291 }
292 const psi = modes.psi;
293 for (let k = 0; k < nModes; k++) {
294 let q = 0;
295 let qd = 0;
a8b5c65Dome example: chaos that keeps its pitch. Constant load with static relaxation (FIRE), full 3D linear modes for curved nets, spectrum-based default frequency range, screening scriptJeremy Magland 296 const off = k * nr;
297 for (let r = 0; r < nr; r++) {
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 298 q += psi[off + r] * bz[r];
299 qd += psi[off + r] * bv[r];
300 }
301 modeEnergy[f * nModes + k] = 0.5 * (qd * qd + om2[k] * q * q);
302 }
303 };
305 // ---- impulses, by audio sample ----------------------------------------------
306 const impulsesAt = new Map<number, typeof sys.hits>();
307 for (const hit of sys.hits) {
308 if (hit.duration > 0) continue;
309 const s = Math.round(hit.t * fs);
310 if (s >= nSamples) continue;
311 if (!impulsesAt.has(s)) impulsesAt.set(s, []);
312 impulsesAt.get(s)!.push(hit);
313 }
314 const tExcitationEnd = Math.max(0, ...sys.hits.map(hh => hh.t + hh.duration));
316 // ---- main loop ---------------------------------------------------------------
317 // Each audio sample interval Delta = substeps * h is one macro step:
318 // D(Delta/2) [substeps conservative steps] D(Delta/2). The damping halves of
319 // neighboring macro steps are merged unless something needs the velocities
320 // at the boundary (a frame record or an impulse).
321 const Delta = 1 / fs;
322 pe = internalForces();
323 externalForces(0);
324 let pendingD = 0;
325 let aborted: string | null = null;
326 let frame = 0;
327 let lastProgress = performance.now();
328 for (let s = 0; s <= nSteps; s++) {
329 const t = s * h;
330 const sample = s % substeps === 0 ? s / substeps : -1;
331 const isFrame = sample >= 0 && sample % frameStride === 0 && frame < nFrames;
332 const imps = sample >= 0 ? impulsesAt.get(sample) : undefined;
333 if ((isFrame || imps) && pendingD !== 0) {
334 applyDamping(pendingD);
335 pendingD = 0;
336 }
337 if (imps) {
338 for (const hit of imps) {
339 for (const b of hit.beads) {
340 if (pinned[b]) continue;
341 v[3 * b] += hit.velocity[0];
342 v[3 * b + 1] += hit.velocity[1];
343 v[3 * b + 2] += hit.velocity[2];
344 }
345 }
346 }
347 if (sample >= 0 && sample < nSamples) audio[sample] = listen();
348 if (isFrame) {
349 frameTimes[frame] = t;
350 energy[frame] = kinetic() + pe - peRest;
351 dissipatedRec[frame] = dissipated;
352 recordModes(frame);
95f2216Minimal chaos example and perturbed-twin chaos test with separation plot and growth-rate fitJeremy Magland 353 opts.onFrame?.(frame, t, x, v);
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 354 if (sample % animStride === 0) {
355 const d = new Float32Array(N3);
356 for (let q = 0; q < N3; q++) d[q] = x[q] - pos0[q];
357 animTimes.push(t);
358 animFrames.push(d);
359 }
360 frame++;
361 }
362 if (s === nSteps) break;
363 if (sample >= 0) {
364 pendingD += Delta / 2;
365 applyDamping(pendingD);
366 pendingD = 0;
367 }
369 let tj = t;
370 for (const c of cs) {
371 const kh = (c * h) / 2;
372 for (let i = 0; i < n; i++) {
373 const im = invm[i] * kh;
374 v[3 * i] += im * F[3 * i];
375 v[3 * i + 1] += im * F[3 * i + 1];
376 v[3 * i + 2] += im * F[3 * i + 2];
377 }
378 const dh = c * h;
379 for (let q = 0; q < N3; q++) x[q] += dh * v[q];
380 tj += c * h;
381 pe = internalForces();
382 externalForces(tj);
383 for (let i = 0; i < n; i++) {
384 const im = invm[i] * kh;
385 v[3 * i] += im * F[3 * i];
386 v[3 * i + 1] += im * F[3 * i + 1];
387 v[3 * i + 2] += im * F[3 * i + 2];
388 }
389 }
390 if ((s + 1) % substeps === 0) pendingD = Delta / 2;
392 if ((s & 1023) === 0) {
393 if (!Number.isFinite(pe)) {
394 aborted = `The solution blew up at t = ${t.toFixed(4)} s. Reduce dt (or safety) or the strike strength.`;
395 break;
396 }
397 const now = performance.now();
398 if (opts.onProgress && now - lastProgress > 200) {
399 opts.onProgress(t);
400 lastProgress = now;
401 }
402 }
403 }
405 // energy drift after the last excitation
406 let energyDrift = 0;
407 let ref = -1;
408 for (let f = 0; f < frame; f++) {
409 if (frameTimes[f] < tExcitationEnd) continue;
410 const total = energy[f] + dissipatedRec[f];
411 if (ref < 0) ref = total;
412 else if (ref > 0) energyDrift = Math.max(energyDrift, Math.abs(total - ref) / ref);
413 }
415 return {
416 audio,
417 sampleRate: fs,
418 dt: h,
419 substeps,
420 omegaMax,
421 frameTimes: frameTimes.slice(0, frame),
422 energy: energy.slice(0, frame),
423 dissipated: dissipatedRec.slice(0, frame),
424 animTimes: Float64Array.from(animTimes),
425 animFrames,
426 modeEnergy: modeEnergy ? modeEnergy.slice(0, frame * nModes) : null,
427 residualForce,
428 maxTension,
429 energyDrift,
430 freqErrors: [1000, 2000, 5000, 10000, 20000]
431 .filter(hz => hz < fs / 2)
432 .map(hz => ({ hz, cents: frequencyErrorCents(sys.method, 2 * Math.PI * hz, h) })),
433 elapsedMs: performance.now() - t0,
434 aborted,
435 warnings,
436 };
437}
439// Highest frequency (rad/s) of the system linearized about the rest positions,
440// by power iteration on M^{-1/2} H M^{-1/2}, with a small margin.
441export function estimateOmegaMax(sys: BeadSystem): number {
a8b5c65Dome example: chaos that keeps its pitch. Constant load with static relaxation (FIRE), full 3D linear modes for curved nets, spectrum-based default frequency range, screening scriptJeremy Magland 442 const { n, si, sj, mass, pinned } = sys;
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 443 const ns = si.length;
a8b5c65Dome example: chaos that keeps its pitch. Constant load with static relaxation (FIRE), full 3D linear modes for curved nets, spectrum-based default frequency range, screening scriptJeremy Magland 444 const B = springBlocks(sys);
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 445 const ism = new Float64Array(n);
446 for (let i = 0; i < n; i++) ism[i] = pinned[i] ? 0 : 1 / Math.sqrt(mass[i]);
447 let w = new Float64Array(3 * n);
448 let seed = 12345;
449 for (let q = 0; q < 3 * n; q++) {
450 seed = (seed * 1103515245 + 12345) & 0x7fffffff;
451 w[q] = ism[Math.floor(q / 3)] * (seed / 0x7fffffff - 0.5);
452 }
453 const u = new Float64Array(3 * n);
454 const y = new Float64Array(3 * n);
455 let lambda = 0;
456 for (let it = 0; it < 500; it++) {
457 let norm = 0;
458 for (let q = 0; q < 3 * n; q++) norm += w[q] * w[q];
459 norm = Math.sqrt(norm);
460 if (norm === 0) return 1;
461 for (let q = 0; q < 3 * n; q++) u[q] = (w[q] / norm) * ism[Math.floor(q / 3)];
462 y.fill(0);
463 for (let s = 0; s < ns; s++) {
464 const i3 = 3 * si[s];
465 const j3 = 3 * sj[s];
466 const r0 = u[i3] - u[j3];
467 const r1 = u[i3 + 1] - u[j3 + 1];
468 const r2 = u[i3 + 2] - u[j3 + 2];
469 const b = 6 * s;
470 const f0 = B[b] * r0 + B[b + 3] * r1 + B[b + 4] * r2;
471 const f1 = B[b + 3] * r0 + B[b + 1] * r1 + B[b + 5] * r2;
472 const f2 = B[b + 4] * r0 + B[b + 5] * r1 + B[b + 2] * r2;
473 y[i3] += f0;
474 y[i3 + 1] += f1;
475 y[i3 + 2] += f2;
476 y[j3] -= f0;
477 y[j3 + 1] -= f1;
478 y[j3 + 2] -= f2;
479 }
480 let rq = 0;
481 const wNext = new Float64Array(3 * n);
482 for (let q = 0; q < 3 * n; q++) {
483 wNext[q] = y[q] * ism[Math.floor(q / 3)];
484 rq += (w[q] / norm) * wNext[q];
485 }
486 w = wNext;
487 if (it > 20 && Math.abs(rq - lambda) < 1e-8 * Math.abs(rq)) {
488 lambda = rq;
489 break;
490 }
491 lambda = rq;
492 }
493 return Math.sqrt(Math.max(lambda, 1e-300)) * 1.02;
494}