31:T435,Elastic waves traveling through a flat plate are called [Lamb waves](https://en.wikipedia.org/wiki/Lamb_waves). They exist in discrete symmetric and antisymmetric modes. The combination of the governing PDE together with stress-free boundary conditions gives the dispersion relation, ${tex.block`\displaystyle {\frac {\tanh(\beta d/2)}{\tanh(\alpha d/2)}}=\left[{\frac {4\alpha \beta k^{2}}{(k^{2}+\beta ^{2})^{2}}}\right]^{\pm 1}`} where symmetric modes correspond to exponent ${tex`+1`} and antisymmetric modes to ${tex`-1`}. Here, ${tex`d`} is the plate thickness and ${tex`\alpha`} and ${tex`\beta`} are given by ${tex.block`\alpha^2 = k^2 - \frac{\omega^2}{c_l^2}`} ${tex.block`\beta^2 = k^2 - \frac{\omega^2}{c_t^2}.`} In the above equations, ${tex`k`} is the wavenumber and ${tex`\omega`} is the temporal frequency. ${tex`c_l`} and ${tex`c_t`} are the longitudinal and transverse wave speeds (p- and s-wave speeds), respectively. And that's everything! From here it's as simple as finding roots to determine the modes. Which turns out to be relatively complicated. 32:T101b,drawMode = new Map( ((regl) => [ { type: "symmetric", func: `vec2 lambDisplacement(vec2 position, float time, vec2 k, float omega, out float fade) { vec2 wt = vec2(omega * time, 0.0); vec2 psi = cexp(cmul(vec2(0, 1), k * position.x - wt)); fade = smoothstep(4.0, 1.0, length(psi)); vec2 k2 = csqr(k); vec2 a2 = cinv(csqr(cl)) * (omega * omega) - k2; vec2 b2 = cinv(csqr(ct)) * (omega * omega) - k2; vec2 a = csqrt(a2); vec2 b = csqrt(b2); vec2 B = cmul(csqr(k) - b2, csin(b * 0.5)); vec2 C = 2.0 * cmul(k, a, csin(a * 0.5)); float mag = max(length(B), length(C)) * length(k) * omega; B /= mag; C /= mag; float y = position.y; vec2 ux = cmul( cmul(B, k, ccos(a * y)) + cmul(C, b, ccos(b * y)), psi, vec2(0, 1)); vec2 uy = cmul(-cmul(B, a, csin(a * y)) + cmul(C, k, csin(b * y)), psi); return vec2(ux.x, uy.x); }` }, { type: "antisymmetric", func: `vec2 lambDisplacement(vec2 position, float time, vec2 k, float omega, out float fade) { vec2 wt = vec2(omega * time, 0.0); vec2 psi = cexp(cmul(vec2(0, 1), k * position.x - wt)); fade = smoothstep(4.0, 1.0, length(psi)); vec2 k2 = csqr(k); vec2 a2 = cinv(csqr(cl)) * (omega * omega) - k2; vec2 b2 = cinv(csqr(ct)) * (omega * omega) - k2; vec2 a = csqrt(a2); vec2 b = csqrt(b2); vec2 A = 2.0 * cmul(k, b, csin(b * 0.5)); vec2 D = cmul(csqr(k) - b2, csin(a * 0.5)); float mag = max(length(A), length(D)) * length(k) * omega; A /= mag; D /= mag; float y = position.y; vec2 ux = cmul(cmul(A, k, csin(a * y)) - cmul(D, b, csin(b * y)), psi, vec2(0, 1)); vec2 uy = cmul(cmul(A, a, ccos(a * y)) + cmul(D, k, ccos(b * y)), psi); return vec2(ux.x, uy.x); }` } ].map(({ type, func }) => [ type, regl({ vert: ` precision highp float; ${glslComplex} attribute vec2 position; uniform vec2 aspect, k, cl, ct; uniform float time, magnitude, omega, speed; uniform bool symmetric; varying vec2 uv; varying float fade; ${func} void main () { uv = position * aspect; vec2 displacement = lambDisplacement(uv, time * 2.0 * speed / omega, k, omega, fade); gl_Position = vec4(position + magnitude * displacement / aspect, 0, 1); }`, frag: ` #extension GL_OES_standard_derivatives : enable precision highp float; varying vec2 uv; varying float fade; float gridFactor (vec2 parameter, float width, float feather) { float w1 = width - feather * 0.5; vec2 d = fwidth(parameter); vec2 looped = 0.5 - abs(mod(parameter, 1.0) - 0.5); vec2 a2 = smoothstep(d * w1, d * (w1 + feather), looped); return min(a2.x, a2.y); } void main () { gl_FragColor = vec4(vec3(0), (1.0 - gridFactor(uv * 16.0, 0.75, 2.0)) * fade); }`, attributes: { position: regl.prop("positionsBuffer") }, uniforms: { aspect: ({ viewportWidth, viewportHeight }) => [ viewportWidth / viewportHeight, 1 ], time: regl.context("time"), k: regl.prop("wavenumber"), cl: regl.prop("cl"), ct: regl.prop("ct"), omega: regl.prop("omega"), magnitude: regl.prop("mag"), speed: regl.prop("speed"), symmetric: (ctx, props) => props.type === "symmetric" }, blend: { enable: true, func: { srcRGB: "src alpha", srcAlpha: 1, dstRGB: "one minus src alpha", dstAlpha: 1 }, equation: { rgb: "add", alpha: "add" } }, depth: { enable: false }, primitive: "triangles", elements: regl.prop("cellsBuffer") }) ]))(regl2) )33:T661,function meshFromFunction(surfaceFn, opts) { let i, j, u, v; opts = opts || {}; const res = opts.resolution === undefined ? 11 : opts.resolution; const nbUFaces = Array.isArray(opts.resolution) ? opts.resolution[0] : res; const nbVFaces = Array.isArray(opts.resolution) ? opts.resolution[1] : res; const uDomain = opts.uDomain === undefined ? [0, 1] : opts.uDomain; const vDomain = opts.vDomain === undefined ? [0, 1] : opts.vDomain; let nbBoundaryAdjustedUFaces = nbUFaces; let nbBoundaryAdjustedVFaces = nbVFaces; if (!opts.vPeriodic) nbBoundaryAdjustedUFaces += 1; if (!opts.uPeriodic) nbBoundaryAdjustedVFaces += 1; const positions = []; const cells = []; for (i = 0; i < nbBoundaryAdjustedUFaces; i++) { u = uDomain[0] + ((uDomain[1] - uDomain[0]) * i) / nbUFaces; for (j = 0; j < nbBoundaryAdjustedVFaces; j++) { v = vDomain[0] + ((vDomain[1] - vDomain[0]) * j) / nbVFaces; positions.push(surfaceFn(u, v)); } } var faceIndex = 0; for (i = 0; i < nbUFaces; i++) { var iPlusOne = i + 1; if (opts.vPeriodic) iPlusOne = iPlusOne % nbUFaces; for (j = 0; j < nbVFaces; j++) { var jPlusOne = j + 1; if (opts.uPeriodic) jPlusOne = jPlusOne % nbVFaces; cells.push([ i * nbBoundaryAdjustedVFaces + j, iPlusOne * nbBoundaryAdjustedVFaces + j, iPlusOne * nbBoundaryAdjustedVFaces + jPlusOne ]); cells.push([ i * nbBoundaryAdjustedVFaces + j, iPlusOne * nbBoundaryAdjustedVFaces + jPlusOne, i * nbBoundaryAdjustedVFaces + jPlusOne ]); } } return { positions, cells }; }34:T1229,computeSelectedPoint = { const DEPTH_BIAS = 10; const el = regl._gl.canvas; function onMove(event) { if (mutable selectedMode.locked) return; const x = (2 * event.offsetX) / el.offsetWidth - 1; const y = 1 - (2 * event.offsetY) / el.offsetHeight; const dx = (40 / el.offsetWidth) * 2; const dy = (40 / el.offsetHeight) * 2; const antiNear = ~viewOpts.indexOf("antisymmetric") ? index.antisymmetric.index.range(x - dx, y - dy, x + dx, y + dy) : []; const symmNear = ~viewOpts.indexOf("symmetric") ? index.symmetric.index.range(x - dx, y - dy, x + dx, y + dy) : []; const data = new Float32Array((antiNear.length + symmNear.length) * 3); let c = 0; let minDist = Infinity; let minSet = null; let minIdx = -1; for (let i = 0; i < antiNear.length; i++) { const idx = antiNear[i] * 3; const kx = antisymmetricLociData[idx]; const ky = antisymmetricLociData[idx + 1]; if (!~viewOpts.indexOf("pure real") && Math.abs(ky) < 1e-4) continue; if (!~viewOpts.indexOf("pure imaginary") && Math.abs(kx) < 1e-4) continue; if ( !~viewOpts.indexOf("complex") && Math.abs(kx) > 1e-4 && Math.abs(ky) > 1e-4 ) continue; const px = index.antisymmetric.projectedVertices[idx]; const py = index.antisymmetric.projectedVertices[idx + 1]; const pz = index.antisymmetric.projectedVertices[idx + 2]; const dx = px - x; const dy = py - y; const r = Math.hypot(dx, dy) + pz * DEPTH_BIAS; if (r < minDist) { minDist = r; minSet = "antisymmetric"; minIdx = idx; } //data[c] = antisymmetricLociData[idx]; //data[c + 1] = antisymmetricLociData[idx + 2]; //data[c + 2] = antisymmetricLociData[idx + 1]; c += 3; } for (let i = 0; i < symmNear.length; i++) { const idx = symmNear[i] * 3; const kx = symmetricLociData[idx]; const ky = symmetricLociData[idx + 1]; if (!~viewOpts.indexOf("pure real") && Math.abs(kx) < 1e-4) continue; if (!~viewOpts.indexOf("pure imaginary") && Math.abs(ky) < 1e-4) continue; if ( !~viewOpts.indexOf("complex") && Math.abs(kx) > 1e-4 && Math.abs(ky) > 1e-4 ) continue; const px = index.symmetric.projectedVertices[idx]; const py = index.symmetric.projectedVertices[idx + 1]; const pz = index.symmetric.projectedVertices[idx + 2]; const dx = px - x; const dy = py - y; const r = Math.hypot(dx, dy) + pz * DEPTH_BIAS; if (r < minDist) { minDist = r; minSet = "symmetric"; minIdx = idx; } //data[c] = symmetricLociData[idx]; //data[c + 1] = symmetricLociData[idx + 2]; //data[c + 2] = symmetricLociData[idx + 1]; c += 3; } if (minIdx !== -1) { const lociData = minSet === "antisymmetric" ? antisymmetricLociData : symmetricLociData; data[0] = lociData[minIdx]; data[1] = lociData[minIdx + 2]; data[2] = lociData[minIdx + 1]; selectedPointsBuffer._count = 1; } else { selectedPointsBuffer._count = 0; } //selectedPointsBuffer._count = antiNear.length + symmNear.length; let modeChanged = false; const freq = data[1]; const kx = data[0]; const ky = data[2]; if ( minSet !== mutable selectedMode.type || freq !== mutable selectedMode.frequency || kx !== mutable selectedMode.wavenumber[0] || ky !== mutable selectedMode.wavenumber[1] ) { mutable selectedMode = mutable selectedMode; mutable selectedMode.type = minSet; mutable selectedMode.frequency = freq; mutable selectedMode.wavenumber[0] = kx; mutable selectedMode.wavenumber[1] = ky; } selectedPointsBuffer(data); camera.taint(); } const p = [NaN, NaN]; function onClick(event) { p[0] = event.offsetX; p[1] = event.offsetY; } function onMouseUp(event) { const dx = p[0] - event.offsetX; const dy = p[1] - event.offsetY; const r = Math.hypot(dx, dy); if (r > 5) return; mutable selectedMode.locked = !mutable selectedMode.locked; mutable selectedMode = mutable selectedMode; } el.addEventListener("mousemove", onMove); el.addEventListener("mousedown", onClick); el.addEventListener("mouseup", onMouseUp); el.addEventListener("wheel", onMove); invalidation.then(() => { el.removeEventListener("mousemove", onMove); el.removeEventListener("mousedown", onClick); el.removeEventListener("mouseup", onMouseUp); el.removeEventListener("wheel", onMove); }); }35:T6f1,lambDispersion = function (ω, parity) { const tmp = ccreateD(); const k2 = ccreateD(); const alpha2 = ccreateD(); const beta2 = ccreateD(); const alpha = ccreateD(); const beta = ccreateD(); const alphahalf = ccreateD(); const betahalf = ccreateD(); const sinha = ccreateD(); const cosha = ccreateD(); const sinhb = ccreateD(); const coshb = ccreateD(); const ab4k2 = ccreateD(); const k2b22 = ccreateD(); const γ2 = csqrS([], cdivSS([], cl.getValue(), ct.getValue())); const cl2inv = cinvS([], csqrS([], cl.getValue())); const ct2inv = cinvS([], csqrS([], ct.getValue())); const ω2 = csqrS([], ω); const lhs = ccreateD(); const rhs = ccreateD(); const symm = parity === "symmetric"; return function fLamb(out, kx, ky) { // k2 := k * k csetValuesD(k2, kx, ky, 1, 0); csqrD(k2, k2); // alpha2 := k2 - w2/cl2 csubDS(alpha2, k2, cmulSS(tmp, ω2, cl2inv)); // beta2 := k2 - w2/ct2 csubDS(beta2, k2, cmulSS(tmp, ω2, ct2inv)); // Use a sqrt which places the branch cut opposite of most of the roots. This increases // root-finding performance since we otherwise can only overcome the branch cut by // excessively subdividing. csqrtAltD(alpha, alpha2); csqrtAltD(beta, beta2); cscaleD(alphahalf, alpha, 0.5); cscaleD(betahalf, beta, 0.5); csinhcoshD(sinha, cosha, alphahalf); csinhcoshD(sinhb, coshb, betahalf); // k2b22 := (k2 + beta2)^2 csqrD(k2b22, caddDD(tmp, k2, beta2)); // ab4k2 := 4 * alpha * beta * k2 cscaleD(ab4k2, cmulDD(tmp, cmulDD(tmp, alpha, beta), k2), 4.0); cmulDD(lhs, cmulDD(tmp, sinhb, cosha), symm ? k2b22 : ab4k2); cmulDD(rhs, cmulDD(tmp, sinha, coshb), symm ? ab4k2 : k2b22); return csubDD(out, lhs, rhs); }; }36:T69e,lambDiffK = function (parity) { const tmp = ccreateD(); const k2 = ccreateD(); const alpha2 = ccreateD(); const beta2 = ccreateD(); const alpha = ccreateD(); const beta = ccreateD(); const alphahalf = ccreateD(); const betahalf = ccreateD(); const sinha = ccreateD(); const cosha = ccreateD(); const sinhb = ccreateD(); const coshb = ccreateD(); const ab4k2 = ccreateD(); const k2b22 = ccreateD(); const cl2inv = cinvS([], csqrS([], cl.getValue())); const ct2inv = cinvS([], csqrS([], ct.getValue())); const ω2 = ccreateS(); const lhs = ccreateD(); const rhs = ccreateD(); const symm = parity === "symmetric"; return function fLamb(out, wx, wy, kx, ky) { // k2 := k * k csetValuesD(k2, kx, ky, 1, 0); csqrD(k2, k2); csqrS(ω2, [wx, wy]); // alpha2 := k2 - w2/cl2 csubDS(alpha2, k2, cmulSS(tmp, ω2, cl2inv)); // beta2 := k2 - w2/ct2 csubDS(beta2, k2, cmulSS(tmp, ω2, ct2inv)); // Use a sqrt which places the branch cut opposite of most of the roots. This increases // performance since we overcome the branch cut by excessively subdividing. csqrtAltD(alpha, alpha2); csqrtAltD(beta, beta2); cscaleD(alphahalf, alpha, 0.5); cscaleD(betahalf, beta, 0.5); csinhcoshD(sinha, cosha, alphahalf); csinhcoshD(sinhb, coshb, betahalf); // k2b22 := (k2 + beta2)^2 csqrD(k2b22, caddDD(tmp, k2, beta2)); // ab4k2 := 4 * alpha * beta * k2 cscaleD(ab4k2, cmulDD(tmp, cmulDD(tmp, alpha, beta), k2), 4.0); cmulDD(lhs, cmulDD(tmp, sinhb, cosha), symm ? k2b22 : ab4k2); cmulDD(rhs, cmulDD(tmp, sinha, coshb), symm ? ab4k2 : k2b22); return csubDD(out, lhs, rhs); }; }37:T6d2,lambDiffW = function (parity) { const tmp = ccreateD(); const k2 = ccreateS(); const alpha2 = ccreateD(); const beta2 = ccreateD(); const alpha = ccreateD(); const beta = ccreateD(); const alphahalf = ccreateD(); const betahalf = ccreateD(); const sinha = ccreateD(); const cosha = ccreateD(); const sinhb = ccreateD(); const coshb = ccreateD(); const ab4k2 = ccreateD(); const k2b22 = ccreateD(); const cl2inv = cinvS([], csqrS([], cl.getValue())); const ct2inv = cinvS([], csqrS([], ct.getValue())); const ω = ccreateD(); const ω2 = ccreateD(); const lhs = ccreateD(); const rhs = ccreateD(); const symm = parity === "symmetric"; return function fLamb(out, ωx, ωy, kx, ky) { // k2 := k * k csetValuesS(k2, kx, ky); csqrS(k2, k2); csetValuesD(ω, ωx, ωy, 1, 0); csqrD(ω2, ω); // alpha2 := k2 - w2/cl2 csubSD(alpha2, k2, cmulDS(tmp, ω2, cl2inv)); // beta2 := k2 - w2/ct2 csubSD(beta2, k2, cmulDS(tmp, ω2, ct2inv)); // Use a sqrt which places the branch cut opposite of most of the roots. This increases // performance since we overcome the branch cut by excessively subdividing. csqrtAltD(alpha, alpha2); csqrtAltD(beta, beta2); cscaleD(alphahalf, alpha, 0.5); cscaleD(betahalf, beta, 0.5); csinhcoshD(sinha, cosha, alphahalf); csinhcoshD(sinhb, coshb, betahalf); // k2b22 := (k2 + beta2)^2 csqrD(k2b22, caddSD(tmp, k2, beta2)); // ab4k2 := 4 * alpha * beta * k2 cscaleD(ab4k2, cmulSD(tmp, k2, cmulDD(tmp, alpha, beta)), 4.0); cmulDD(lhs, cmulDD(tmp, sinhb, cosha), symm ? k2b22 : ab4k2); cmulDD(rhs, cmulDD(tmp, sinha, coshb), symm ? ab4k2 : k2b22); return csubDD(out, lhs, rhs); }; }38:T1f5f,findComplexRoots = { function newtonsIdentities(er, ei, M, p) { var i, j, j2, sgn, sgn2; // Use Newton's Identities to construct a polynomial, the roots of // which match our complex analytic function: er[0] = 1; ei[0] = 0; er[1] = p[0]; ei[1] = p[1]; for (i = 2, sgn = -1; i <= M; i++, sgn = -sgn) { er[i] = p[2 * i - 2] * sgn; ei[i] = p[2 * i - 1] * sgn; for (j = 1, sgn2 = -sgn; j < i; j++, sgn2 = -sgn2) { j2 = 2 * (i - j - 1); er[i] += (er[j] * p[j2] - ei[j] * p[j2 + 1]) * sgn2; ei[i] += (er[j] * p[j2 + 1] + ei[j] * p[j2]) * sgn2; } er[i] /= i; ei[i] /= i; } sgn = M % 2 === 0 ? 1 : -1; for (i = 0; i <= M; i++, sgn = -sgn) { er[i] *= sgn; ei[i] *= sgn; } er.reverse(); ei.reverse(); } var out = new Array(2); var params = {}; function cmod(a, b) { var c; var aa = Math.abs(a); var ab = Math.abs(b); if (aa === 0 && ab === 0) return 0; if (aa >= ab) { c = b / a; return aa * Math.sqrt(1 + c * c); } c = a / b; return ab * Math.sqrt(1 + c * c); } function cdiv(out, a, b, c, d) { var e, f; if (Math.abs(c) >= Math.abs(d)) { e = d / c; f = 1 / (c + d * e); out[0] = (a + b * e) * f; out[1] = (b - a * e) * f; } else { e = c / d; f = 1 / (c * e + d); out[0] = (a * e + b) * f; out[1] = (b * e - a) * f; } return out; } var dzrsub = [ 0, 0.81, 0.5727564927611035, 0, -0.5727564927611035, -0.81, -0.5727564927611035, 0, 0.5727564927611035 ]; var dzisub = [ 0, 0, 0.5727564927611035, 0.81, 0.5727564927611035, 0, -0.5727564927611035, -0.81, -0.5727564927611035 ]; var rsub = [ 0.5, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667, 0.4166666666666667 ]; function s0Integrand(theta) { var c, s, a, b, rtmp, itmp; c = Math.cos(theta); s = Math.sin(theta); a = params.z0r + c * params.r; b = params.z0i + s * params.r; // tmp <- f(z) params.f(out, a, b); rtmp = out[0]; itmp = out[1]; // out <- f'(z) out[0] = out[2]; out[1] = out[3]; // out /= f(z) cdiv(out, out[0], out[1], rtmp, itmp); // re(tmp * dz / (2 pi i)) return out[0] * c - out[1] * s; } function smIntegrand(out, theta) { var c, s, a, b, rtmp, itmp; c = Math.cos(theta); s = Math.sin(theta); // out <- f(z) params.f(out, params.z0r + c * params.r, params.z0i + s * params.r); // out <- f'(z) / f(z) cdiv(out, out[2], out[3], out[0], out[1]); // out *= dz (/ (2 pi) but leave this part until outside the loop) rtmp = out[0]; out[0] = out[0] * c - out[1] * s; out[1] = rtmp * s + out[1] * c; // out *= z rtmp = out[0]; out[0] = out[0] * c - out[1] * s; out[1] = rtmp * s + out[1] * c; // Recursively multiply by z to compute successive moments z^N: // out[i] = out[i - 1] * z for (var i = 2; i < 2 * params.M; i += 2) { out[i] = out[i - 2] * c - out[i - 1] * s; out[i + 1] = out[i - 2] * s + out[i - 1] * c; } return out; } function locateZeros( qr, qi, f, zr, zi, r, z0r, z0i, r0, tol, maxIntegrationDepth, s0IntegrationTolerance, s0IntegrationDepth, curDepth, maxRootsPerCircle, searchCircleOutput ) { var p, s0, er, ei, i, j, qri, qii, qmod, qC, M, rr; var qNr, qNi, tmp, qNew, success; // Exit if we've reached the recursion limit: if (curDepth <= 0) { return false; } // Exit early if this circle is completely outside the search region if (cmod(zr - z0r, zi - z0i) - r > r0) { return true; } params.z0r = zr; params.z0i = zi; params.r = r; params.f = f; s0 = adaptiveSimpson( s0Integrand, 0, Math.PI * 2, s0IntegrationTolerance, s0IntegrationDepth ) * r * (0.5 / Math.PI); M = Math.round(s0); qC = 0; for (i = qr.length - 1; i >= 0; i--) { // Count this zero if it's inside the circle: if (cmod(qr[i] - zr, qi[i] - zi) < r) qC++; } if (M - qC <= 0.9 || s0 < 0.9) { } else if (Math.abs(s0 - M) > 1e-1 || M - qC > maxRootsPerCircle) { // If this happens, then we had a pretty inaccurate integration: for (i = 0, success = true; i < 9; i++) { rr = r * rsub[i]; success = locateZeros( qr, qi, f, zr + r * dzrsub[i], zi + r * dzisub[i], rr, z0r, z0i, r0, tol, maxIntegrationDepth, s0IntegrationTolerance, s0IntegrationDepth, curDepth - 1, maxRootsPerCircle, searchCircleOutput ) && success; } } else { // Compute the zeros: params.M = M - qC; var p = vectorAdaptiveSimpson( [], smIntegrand, 0, Math.PI * 2, tol, maxIntegrationDepth, 2 * params.M ); // post-multiply the integral by this constant factor var f = (r * 0.5) / Math.PI; for (var i = 2 * params.M - 1; i >= 0; i--) { p[i] *= f; } // Deflate the result by known zeros inside the current circle: for (i = qr.length - 1; i >= 0; i--) { // Scale to the unit circle in which we're working qri = (qr[i] - zr) / r; qii = (qi[i] - zi) / r; if (cmod(qri, qii) < 1) { // Initialize the powers of qN var qNr = qri; var qNi = qii; for (j = 0, qNr = 1, qNi = 0; j < p.length; j += 2) { tmp = qNr; qNr = tmp * qri - qNi * qii; qNi = tmp * qii + qNi * qri; // Subtract off the power of the zero: p[j] -= qNr; p[j + 1] -= qNi; } } } // Apply Newton's Identities, constructing a polynomial sharing the roots // of the analytic function: er = []; ei = []; newtonsIdentities(er, ei, params.M, p); // Compute the roots of the constructed polynomial: qNew = findRoots(er, ei, 100 * params.M * params.M, tol); // Append and filter to zeros inside the original contour: for (i = qNew[0].length - 1; i >= 0; i--) { var qrRescaled = zr + qNew[0][i] * r; var qiRescaled = zi + qNew[1][i] * r; //if (cmod(qrRescaled - z0r, qiRescaled - z0i) < r0) { qr.push(qrRescaled); qi.push(qiRescaled); //} } } if (searchCircleOutput) { searchCircleOutput.push({ z: [zr, zi], r: r, text: "s₀ = " + s0 .toFixed(5) .toString() .replace(/(\.[0-9]*[1-9]|\.0)0*$/, "$1") .replace(/\.0$/, "") .replace(/^-0$/, "0") }); } return [qr, qi]; } return function delvesLyness( f, z0, r, { s0IntegrationDepth = 13, s0IntegrationTolerance = 1e-6, maxIntegrationDepth = 12, integrationTolerance = 1e-8, maxSearchDepth = 6, maxRootsPerCircle = 3, searchOutput = null } = {} ) { if (searchOutput) searchOutput.length = 0; var foundRoots = locateZeros( [], [], f, z0[0], z0[1], r, z0[0], z0[1], r, integrationTolerance, maxIntegrationDepth, s0IntegrationTolerance, s0IntegrationDepth, maxSearchDepth, maxRootsPerCircle, searchOutput ); var output = []; for (var i = 0; i < foundRoots[0].length; i++) { if (cmod(foundRoots[0][i] - z0[0], foundRoots[1][i] - z0[1]) <= r) { output.push([foundRoots[0][i], foundRoots[1][i]]); } } return output; }; }39:T6c9,function traceRoots( roots, parity, { tolerance = 1e-8, maxIterations = 3, fmax = 4 } = {} ) { const t0 = performance.now(); const wr0 = 2.0 * Math.PI * f0 * d; const wi0 = 0; const diffW = lambDiffW(parity); const diffK = lambDiffK(parity); const dw = ccreateS(); const dk = ccreateS(); const deltaK = ccreateS(); const integrationOptions = { tolerance: 1e-6 }; const ans = roots.map((root) => { let wr = f0 * 2.0 * Math.PI * d; let wi = 0; let kr = root[0]; let ki = root[1]; const path = [[kr, ki, f0]]; const state = { y: [kr, ki, wr], t: 0, dt: 0.01 }; for (let i = 0; i < 400; i++) { const deriv = (yp, y, t) => { const kr = y[0]; const ki = y[1]; const wr = y[2]; const wi = 0; diffK(dk, wr, wi, kr, ki); diffW(dw, wr, wi, kr, ki); cdivSS(deltaK, [-dw[2], -dw[3]], [dk[2], dk[3]]); const l = Math.sqrt(1 + deltaK[0] * deltaK[0] + deltaK[1] * deltaK[1]); const ds = 0.001 / l; yp[0] = deltaK[0] * ds; yp[1] = deltaK[1] * ds; yp[2] = ds; }; ode45(state, deriv, integrationOptions); const cur = [state.y[0], state.y[1]]; complexNewtonRaphson( cur, (out, kr, ki) => diffK(out, state.y[2], 0, kr, ki), cur, { tolerance, maxIterations } ); state.y[0] = cur[0]; state.y[1] = cur[1]; const f = state.y[2] / (2.0 * Math.PI * d); path.push([state.y[0], state.y[1], f]); if (f > fmax) break; } path.unshift([0, 0, -1]); path.push([0, 0, -1]); return path; }); const t1 = performance.now(); mutable traceTimeMS = t1 - t0; return ans; }3a:T7c8,drawLines = reglLines(regl, { vert: ` precision highp float; #pragma lines: attribute vec3 position; #pragma lines: attribute float orientation; #pragma lines: position = projectPoint(position); #pragma lines: width = pointWidth(); #pragma lines: varying float computedWidth = pointWidth(); uniform mat4 uProjectionView, model, modelViewProjection; uniform float width; uniform vec3 types; vec4 projectPoint (vec3 point) { const float tol = 1e-2; if (point.z < 0.0) return vec4(0); if ((abs(point.x) < tol && types.y == 0.0) || (abs(point.y) < tol && types.x == 0.0) || (abs(point.x) > tol && abs(point.y) > tol && types.z == 0.0)) { return vec4(0); } return modelViewProjection * vec4(point.xzy, 1); } float pointWidth () { return width; }`, frag: ` precision highp float; varying vec3 lineCoord; varying float computedWidth; uniform vec2 borderWidth; uniform vec4 borderColor, color; float linearstep(float x0, float x1, float x) { return clamp((x - x0) / (x1 - x0), 0.0, 1.0); } void main() { float sdf = length(lineCoord.xy); vec2 borderThreshold = 1.0 - borderWidth / computedWidth; float isBorder = linearstep(borderThreshold.x, borderThreshold.y, sdf); gl_FragColor = vec4(mix(color.rgb, borderColor.rgb, isBorder * borderColor.a * color.a), 1); }`, uniforms: { modelViewProjection: regl.prop("modelViewProjection"), borderColor: regl.prop("borderColor"), width: (ctx, props) => props.width * ctx.pixelRatio, color: regl.prop("color"), borderWidth: (ctx, props) => [ 2.0 * props.borderWidth * ctx.pixelRatio + 0.75, 2.0 * props.borderWidth * ctx.pixelRatio - 0.75 ], types: (ctx, { viewOpts }) => [ ~viewOpts.indexOf("pure real"), ~viewOpts.indexOf("pure imaginary"), ~viewOpts.indexOf("complex") ] }, cull: { enable: false } })3b:T849,function renderFrame({ dTheta = 0, tick = 0 } = {}) { if (dTheta) camera.params.dTheta = dTheta; const lineProps = { insertCaps: true, join: "round", cap: "round", width: 3, borderWidth: 0.75, borderColor: [1, 1, 1, 1], model, viewOpts }; camera(({ dirty }) => { if (!dirty) return; mat4multiply( modelViewProjection, camera.state.projection, mat4multiply(modelViewProjection, camera.state.view, model) ); updateIndex(modelViewProjection); camera.params.rotationCenter = camera.params.center; regl.clear({ color: [1, 1, 1, 1], depth: 1 }); drawFancyAxes(); if ( false && index.antisymmetric && index.antisymmetric.projectedVerticesBuffer ) { drawPoints([ { position: index.antisymmetric.projectedVerticesBuffer, count: index.antisymmetric.projectedVertices.length / 3, color: [45 / 255, 160 / 255, 211 / 255] }, { position: index.symmetric.projectedVerticesBuffer, count: index.symmetric.projectedVertices.length / 3, color: [211 / 255, 45 / 255, 160 / 255] } ]); } if (~viewOpts.indexOf("symmetric")) { drawLines({ ...lineProps, modelViewProjection, color: [211 / 255, 45 / 255, 160 / 255, 1], vertexCount: symmetricLociData.length / 3, vertexAttributes: { position: symmetricLociBuffer } }); } if (~viewOpts.indexOf("antisymmetric")) { drawLines({ ...lineProps, modelViewProjection, color: [45 / 255, 160 / 255, 211 / 255, 1], vertexCount: antisymmetricLociData.length / 3, vertexAttributes: { position: antisymmetricLociBuffer } }); } drawPoints({ position: selectedPointsBuffer, color: selectedMode.type === "antisymmetric" ? [45 / 255, 160 / 255, 211 / 255] : [211 / 255, 45 / 255, 160 / 255], count: selectedPointsBuffer._count || 0, matrix: modelViewProjection }); }); }3c:Te81,drawFancyAxes = { const boxVertices = regl.buffer( [0, 1, 2].map(() => [ [-1, -1, -1], [1, -1, -1], [1, 1, -1], [-1, 1, -1], [-1, -1, 1], [1, -1, 1], [1, 1, 1], [-1, 1, 1] ]) ); const boxFaceNormal = regl.buffer( [0, 1, 2].map((i) => [...Array(8).keys()].map(() => [0, 1, 2].map((j) => (i === (5 - j) % 3 ? 1 : 0)) ) ) ); const boxFaceTangent = regl.buffer( [0, 1, 2].map((i) => [...Array(8).keys()].map(() => [0, 1, 2].map((j) => (i === (3 - j) % 3 ? 1 : 0)) ) ) ); const boxFaceBitangent = regl.buffer( [0, 1, 2].map((i) => [...Array(8).keys()].map(() => [0, 1, 2].map((j) => (i === (4 - j) % 3 ? 1 : 0)) ) ) ); // prettier-ignore const boxFaceElements = regl.elements([ [0, 1, 2], [0, 2, 3], [6, 4, 7], [6, 5, 4], [0, 4, 5].map(i => i + 8), [0, 5, 1].map(i => i + 8), [2, 7, 3].map(i => i + 8), [2, 6, 7].map(i => i + 8), [1, 5, 2].map(i => i + 16), [2, 5, 6].map(i => i + 16), [0, 3, 7].map(i => i + 16), [0, 7, 4].map(i => i + 16), ]); invalidation.then(() => { boxVertices.destroy(); boxFaceElements.destroy(); boxFaceNormal.destroy(); boxFaceTangent.destroy(); boxFaceBitangent.destroy(); }); const transform = mat4create(); const normalMatrix = mat3create(); const drawFaces = regl({ vert: ` precision highp float; uniform mat4 uProjectionView, transform; uniform mat3 normalMatrix; attribute vec3 position, tangent, bitangent, normal; varying vec2 coord; varying vec3 p, n, t, b; void main () { p = (transform * vec4(position, 1)).xyz; n = normalMatrix * normal; t = normalMatrix * tangent; b = normalMatrix * bitangent; coord = vec2(dot(position, b), dot(position, t)); gl_Position = uProjectionView * vec4(p, 1); }`, frag: ` #extension GL_OES_standard_derivatives : enable precision mediump float; varying vec2 coord; varying vec3 p, n, t, b; uniform vec3 eye; float gridFactor (vec2 parameter, float width) { vec2 d = fwidth(parameter); vec2 looped = 0.5 - abs(mod(parameter, 1.0) - 0.5); vec2 a2 = smoothstep(d * (width + 0.5), d * (width - 0.5), looped); return max(a2.x, a2.y); } void main () { float vDotN = dot(normalize(p - eye), n); float fade = smoothstep(0.01, 0.05, abs(vDotN)); float alpha = fade * 0.25 * ( gridFactor(coord * 10.0, 0.5) + 0.2 * gridFactor(coord * 10.0 * 10.0, 0.5) ); gl_FragColor = vec4(vec3(0), alpha); if (alpha < 0.001) discard; }`, attributes: { position: boxVertices, normal: boxFaceNormal, tangent: boxFaceTangent, bitangent: boxFaceBitangent }, uniforms: { transform: mat4scale( mat4create(), mat4fromTranslation(transform, [0, 10, 0]), [10, 10, 10] ), normalMatrix: mat3normalFromMat4( normalMatrix, mat4scale( mat4create(), mat4translate(mat4create(), mat4create(), [0, 10, 0]), [10, 10, 10] ) ), eye: regl.context("eye") }, depth: { enable: true, mask: false }, cull: { enable: true, face: "back" }, blend: { enable: true, func: { srcRGB: "src alpha", dstRGB: "one minus src alpha", srcAlpha: 1, dstAlpha: 1 }, equation: { rgb: "add", alpha: "add" } }, primitive: "triangles", elements: boxFaceElements }); function labelMaker() {} return function () { drawFaces(); }; }3d:[] 2f:{"title":"Tracing Lamb Wave Modes in the Complex Plane","cells":[{"id":0,"value":"# Tracing Lamb Wave Modes in the Complex Plane","pinned":false,"hidden":false,"mode":"md"},{"id":1890,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":1905,"value":"First, the above function is meromorphic. It has poles, zeros, and branch cuts. Gross. We can cross-multiply and rewrite the function to remove poles. For symmetric modes,\n\n${tex.block`\\displaystyle \n\n\\begin{aligned}\nf(\\omega, k) = &\\;\n{\\sinh(\\beta d/2) \\cosh(\\alpha d/2) (k^{2}+\\beta ^{2})^{2}}\\\\\n& -\n\\sinh(\\alpha d/2) \\cosh(\\beta d/2) 4\\alpha \\beta k^{2}\n\\end{aligned}\n`}\n\n We can construct an equivalent relation for antisymmetric modes. We'll still have branch cuts via ${tex`\\alpha`} and ${tex`\\beta`}, but this is a significant improvement.\n\nThe above equation is solved when\n\n${tex.block`f(\\omega, k) = 0.`}","pinned":false,"hidden":false,"mode":"md"},{"id":2000,"value":"## Computation method\n\nThe plots above illustrate solutions of this equation in the complex wavenumber plane, tabulated over a range of frequencies. All math is performed using automatic differentiation to compute derivatives together with values. This is implemented in [@rreusser/complex-auto-differentiation](https://observablehq.com/@rreusser/complex-auto-differentiation). Derivative-free methods exist, but this makes some tasks like Newton iteration somewhat easier.\n","pinned":false,"hidden":false,"mode":"md"},{"id":2002,"value":"\n### Roots for a single frequency slice\n\nFirst, we compute roots for a single frequency slice.\n\n- Roots are computed using the Delves-Lyness method via contour integration. I've written up an explanation in [Finding Roots in the Complex Plane](https://observablehq.com/@rreusser/finding-roots-in-the-complex-plane)\n- Roots are refined using Newton iteration.\n","pinned":false,"hidden":false,"mode":"md"},{"id":3241,"value":"\n### Tracing roots over frequency\n\nOnce computed, roots are traced through varying frequency. If we ran the above method over and over, independently for each frequency slice, we'd find two things. First, it would take a long time. Second, we'd have discrete sets of roots at each frequency. We could try to pair them greedily and construct a continuous locus of roots, but this is slow and error-prone. Instead, we trace them.\n","pinned":false,"hidden":false,"mode":"md"},{"id":3243,"value":"\nGiven a root at ${tex`(\\omega_0, k_0),`} the basic method is to perfom a first-order Taylor series expansion,\n","pinned":false,"hidden":false,"mode":"md"},{"id":3245,"value":"\n${tex.block`f(\\omega_0+\\Delta \\omega, k_0 + \\Delta k) = \\frac{\\partial f}{\\partial \\omega} \\Delta \\omega + \\frac{\\partial f}{\\partial k} \\Delta k.`}\n","pinned":false,"hidden":false,"mode":"md"},{"id":3247,"value":"\nFor brevity, evaluation of derivatives at ${tex`(\\omega_0, k_0)`} is implied. When tracing roots, we seek to maintain ${tex`f(\\omega, k) = 0,`} so that we set the ${tex`f`} term above to zero and rearrange, yielding\n","pinned":false,"hidden":false,"mode":"md"},{"id":3249,"value":"\n${tex.block`\\frac{\\partial f}{\\partial \\omega} \\Delta \\omega = -\\frac{\\partial f}{\\partial k} \\Delta k.`}\n","pinned":false,"hidden":false,"mode":"md"},{"id":3251,"value":"\nThe direction ${tex`(\\Delta \\omega, \\Delta k)`} represents the tangent direction at which ${tex`f(\\omega, k) = 0`} is maintained. Thus, we frame this as an ODE and integrate.","pinned":false,"hidden":false,"mode":"md"},{"id":3253,"value":"\nI've used my adaptive ODE45 implementation from [@rreusser/integration](https://observablehq.com/@rreusser/integration).\n","pinned":false,"hidden":false,"mode":"md"},{"id":1961,"value":"\nThere is some trickiness, for example when two roots meet along the complex axis and then immediately diverge along the real axis. If we're not careful, we'll follow the wrong branch for one or the other. To work around this, I've added a very small amount of viscoelasticity to bias roots toward one direction or the other. \n","pinned":false,"hidden":false,"mode":"md"},{"id":3255,"value":"Finally, since exact equality with zero is not maintained, I've used a couple passes of complex Newton iteration at the very end. Newton iteration in the complex plane is also quite error-prone, yielding the well-known [Newton fractal](https://en.wikipedia.org/wiki/Newton_fractal). Although Newton iteration finds roots just like it does for real-valued functions, it falls into fractal patterns if you start iteration too far away from the neighborhood of a root. Fortunately, errors in the integration step are quite small, so that Newton iteration works quite reliably.","pinned":false,"hidden":false,"mode":"md"},{"id":171,"value":"## Material properties\n\nWithout loss of generality since the final result differs only by dimensional scaling, I've chosen unit properties, with the exception of [Poisson's ratio](https://en.wikipedia.org/wiki/Poisson%27s_ratio), ${tex`\\nu`}, which is non-dimensional.","pinned":false,"hidden":false,"mode":"md"},{"id":5,"value":"viewof E = scifmt`E = ${1} Pa`","pinned":false,"hidden":false,"mode":"ojs"},{"id":47,"value":"viewof ρ = fmt`ρ = ${1} kg/m³`","pinned":false,"hidden":false,"mode":"ojs"},{"id":15,"value":"viewof ν = fmt`ν = ${0.3}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":112,"value":"viewof d = fmt`d = ${1} m`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1882,"value":"## Roots for a single frequency slice","pinned":false,"hidden":false,"mode":"md"},{"id":3260,"value":"The plot below shows the initial roots over complex-valued wavenumber ${tex`k`} at a single frequency slice ${tex`\\omega`}. The roots are computed with contour integration and refined with Newton iteration.","pinned":false,"hidden":false,"mode":"md"},{"id":211,"value":"viewof f = Inputs.range([0.005, 4], {\n step: 0.0001,\n value: 1,\n label: \"Frequency, f, Hz\"\n})\n/*viewof f = Scrubber(\n [...Array(200).keys()].map((f) => 0.01 + f * 0.01),\n { autoplay: false }\n)*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":2006,"value":"viewof radius = Inputs.range([1, 60], {\n label: \"Initial search radius\",\n value: 40\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":161,"value":"viewof parity = Inputs.radio([\"symmetric\", \"antisymmetric\"], {\n value: \"symmetric\",\n label: \"Parity\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":374,"value":"viewof pl = {\n const pl = plot([\n DomainColoringLayer(fLambDispersion, \"z\", {\n ...domainColoringConfig,\n shadeRange: [20, 55]\n }),\n ...[wd, cl, ct],\n ...roots[parity].map((r) => ComplexVariable({ value: r }))\n ]);\n const el = html`

${pl}

Constant-frequency slice of the Lamb wave dispersion relation.

`;\n el.value = pl.value;\n return el;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1884,"value":"## Roots traced through varying frequency","pinned":false,"hidden":false,"mode":"md"},{"id":3263,"value":"Since the method above is somewhat expensive, we use the above roots as seeds and proceed to tracing them as frequency ${tex`\\omega`} varies, rather than re-computing roots from scratch for each frequency slice. The result is the roots locus below over complex wavenumber ${tex`k`} in the horizontal plane and frequency ${tex`\\omega`} in the vertical dimension.","pinned":false,"hidden":false,"mode":"md"},{"id":1336,"value":"viewof viewOpts = Inputs.checkbox(\n [\"symmetric\", \"antisymmetric\", \"pure real\", \"pure imaginary\", \"complex\"],\n {\n value: [\n \"symmetric\",\n \"antisymmetric\",\n \"pure real\",\n \"pure imaginary\",\n \"complex\"\n ]\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":181,"value":"viewof θ = Inputs.range([1e-8, Math.PI], {\n step: 0.00001,\n value: 1e-8,\n label: \"Viscoelastic loss angle, radians\",\n transform: Math.log\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":890,"value":"viewof regl = {\n const regl = reglCanvas(this && this.querySelector(\"canvas\"), {\n width: viewport.width,\n height: viewport.height,\n pixelRatio: viewport.pixelRatio,\n attributes: {\n preserveDrawingBuffer: true\n },\n extensions: [\n \"OES_standard_derivatives\",\n \"ANGLE_instanced_arrays\",\n \"OES_texture_float\"\n ],\n optionalExtensions: [\"OES_vertex_array_object\"]\n });\n regl.value._gl.canvas.style.cursor = \"crosshair\";\n const el = html`

${regl}

Symmetric (red) and antisymmetric (blue) Lamb wave modes. Click to select/deselect a mode.

`;\n el.value = regl.value;\n return el;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2490,"value":"viewof regl2 = {\n const regl = reglCanvas(null, {\n width: viewport.width,\n height: Math.round(viewport.width * 0.2),\n pixelRatio: viewport.pixelRatio,\n attributes: {\n preserveDrawingBuffer: true\n },\n extensions: [\n \"OES_standard_derivatives\",\n \"ANGLE_instanced_arrays\",\n \"OES_texture_float\",\n \"OES_element_index_uint\"\n ],\n optionalExtensions: [\"OES_vertex_array_object\"]\n });\n\n const el = html`

${regl}

Displacement visualized along the plate cross section.

`;\n el.value = regl.value;\n return el;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2751,"value":"viewof mag = Inputs.range([0.01, 1], {\n transform: Math.log,\n label: \"Magnitude\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":3047,"value":"viewof speed = Inputs.range([0.1, 10], {\n transform: Math.log,\n label: \"Speed\",\n value: 1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2765,"value":"viewof animate = Inputs.checkbox([\"animate\"], { value: [\"animate\"] })","pinned":false,"hidden":false,"mode":"ojs"},{"id":2412,"value":"mutable selectedMode = ({\n type: null,\n frequency: NaN,\n wavenumber: [NaN, NaN],\n locked: false\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2493,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":2510,"value":"drawLoop2 = {\n if (!~animate.indexOf(\"animate\")) return;\n let frame = regl2.frame(() => {\n try {\n if ([\"symmetric\", \"antisymmetric\"].indexOf(selectedMode.type) === -1)\n return;\n regl2.clear({ color: [1, 1, 1, 1], depth: 1 });\n drawMode.get(selectedMode.type)({\n ...mesh,\n wavenumber: selectedMode.wavenumber,\n omega: selectedMode.frequency * 2.0 * Math.PI,\n cl: cl.getValue(),\n ct: ct.getValue(),\n mag,\n speed\n });\n } catch (e) {\n console.error(e);\n if (frame) frame.cancel();\n frame = null;\n throw e;\n }\n });\n invalidation.then(() => {\n if (frame) frame.cancel();\n frame = null;\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2504,"value":"mesh = {\n const mesh = meshFromFunction((u, v) => [u, v], {\n uDomain: [-1.5, 1.5],\n vDomain: [-0.5, 0.5],\n resolution: [201, 40],\n uPeriodic: false,\n vPeriodic: false\n });\n mesh.positionsBuffer = regl2.buffer(mesh.positions);\n mesh.cellsBuffer = regl2.elements(mesh.cells);\n invalidation.then(() => {\n mesh.positionsBuffer.destroy();\n mesh.cellsBuffer.destroy();\n });\n return mesh;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2500,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":2188,"value":"function updateIndex(matrix) {\n if (\n mat4equals(index.matrix, modelViewProjection)\n //|| !index.antisymmetric.projectedVerticesBuffer\n ) {\n return;\n }\n\n mat4copy(index.matrix, matrix);\n index.antisymmetric = makeIndex(antisymmetricLociData, matrix);\n index.symmetric = makeIndex(symmetricLociData, matrix);\n\n index.antisymmetric.projectedVerticesBuffer = (\n index.antisymmetric.projectedVerticesBuffer || regl.buffer\n )(index.antisymmetric.projectedVertices);\n\n index.symmetric.projectedVerticesBuffer = (\n index.symmetric.projectedVerticesBuffer || regl.buffer\n )(index.symmetric.projectedVertices);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2352,"value":"selectedPointsBuffer = regl.buffer([0, 0, 0])","pinned":false,"hidden":false,"mode":"ojs"},{"id":2112,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":2077,"value":"function makeIndex(verts, matrix) {\n const index = new KDBush(verts.length / 3);\n const p = new Float32Array(3);\n const pp = new Float32Array(3);\n const projectedVertices = new Float32Array(verts.length);\n for (let i = 0; i < verts.length; i += 3) {\n p[0] = verts[i];\n p[1] = verts[i + 2];\n p[2] = verts[i + 1];\n vec3transformMat4(pp, p, matrix);\n if (isFinite(pp[0]) && pp[2] <= 1 && pp[2] >= -1) {\n index.add(pp[0], pp[1]);\n projectedVertices[i] = pp[0];\n projectedVertices[i + 1] = pp[1];\n projectedVertices[i + 2] = pp[2];\n } else {\n index.add(-2, -2);\n projectedVertices[i] = -2;\n projectedVertices[i + 1] = -2;\n }\n }\n index.finish();\n return { index, projectedVertices };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2183,"value":"index = ({\n matrix: mat4create(),\n antisymmetric: {},\n symmetric: {}\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2081,"value":"KDBush = (await import(\"https://cdn.skypack.dev/kdbush\")).default","pinned":false,"hidden":false,"mode":"ojs"},{"id":1260,"value":"f0 = 0.006 // Root seed frequency","pinned":false,"hidden":false,"mode":"ojs"},{"id":1414,"value":"mutable rootsTimeMS = 0","pinned":false,"hidden":false,"mode":"ojs"},{"id":1848,"value":"mutable traceTimeMS = 0","pinned":false,"hidden":false,"mode":"ojs"},{"id":168,"value":"## Viscoelastic properties","pinned":false,"hidden":false,"mode":"md"},{"id":391,"value":"viewof wd = ComplexVariable({\n name: \"wd\",\n value: [2.0 * Math.PI * f * d, 0],\n visible: false\n}).show()","pinned":false,"hidden":false,"mode":"ojs"},{"id":395,"value":"viewof cl = {\n const Ev = cscaleS([], [Math.cos(-θ), Math.sin(-θ)], E);\n const cl = cscaleS(\n [],\n csqrtS([], Ev),\n Math.sqrt((1 - ν) / ρ / (1 + ν) / (1 - 2 * ν))\n );\n return ComplexVariable({ name: \"c_l\", value: cl, visible: false }).show();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":413,"value":"viewof ct = {\n const Ev = cscaleS([], [Math.cos(-θ), Math.sin(-θ)], E);\n const ct = cscaleS([], csqrtS([], Ev), 1 / Math.sqrt(2 * ρ * (1 + ν)));\n return ComplexVariable({ name: \"c_t\", value: ct, visible: false }).show();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1955,"value":"## Implementation","pinned":false,"hidden":false,"mode":"md"},{"id":105,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":1487,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":1511,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":380,"value":"fLambDispersion = {\n const ɑ2 = `z^2 - (wd/c_l)^2`;\n const β2 = `z^2 - (wd/c_t)^2`;\n const ɑ = `(i * sqrt(-(${ɑ2})))`;\n const β = `(i * sqrt(-(${β2})))`;\n const k2b22 = `((z^2 + ${β2})^2)`;\n const ab4k2 = `(4 * ${ɑ} * ${β} * z^2)`;\n return `\nsinh(${β}/2) * cosh(${ɑ}/2) * ${parity === \"symmetric\" ? k2b22 : `-${ab4k2}`} -\nsinh(${ɑ}/2) * cosh(${β}/2) * ${parity === \"symmetric\" ? ab4k2 : `-${k2b22}`}`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":307,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":622,"value":"function complexNewtonRaphson(\n out,\n f,\n z0,\n { tolerance = 1e-5, maxIterations = 50 } = {}\n) {\n var iteration = maxIterations;\n var delta = 1,\n previousDelta = 0;\n var tmp = [0, 0, 0, 0];\n var zre = z0[0];\n var zim = z0[1];\n while (delta / previousDelta > tolerance && --iteration > 0) {\n previousDelta = delta;\n\n // Compute the function and its derivative tmp <- [f_re, f_im, f'_re, f'_im]\n f(tmp, zre, zim);\n\n // Divide tmp <- f / f'\n cdivSS(tmp, [tmp[0], tmp[1]], [tmp[2], tmp[3]]);\n\n // Update the current guess\n zre -= tmp[0];\n zim -= tmp[1];\n\n // Check convergence\n delta = Math.hypot(tmp[0], tmp[1]);\n }\n // Return the answer if converged\n if (iteration > 0 || delta / previousDelta <= tolerance) {\n out[0] = zre;\n out[1] = zim;\n } else {\n // Well, even if not converged, still just return the best we got\n out[0] = zre;\n out[1] = zim;\n }\n return out;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":332,"value":"function computeRoots(f0, { r0 = 3, z0 = [0, 0] } = {}) {\n const wd = [2.0 * Math.PI * f0 * d, 0];\n function computeRoots(parity) {\n const t0 = performance.now();\n const dispersion = lambDispersion(wd, parity);\n let roots = findComplexRoots(dispersion, z0, r0, {\n s0IntegrationDepth: 10,\n s0IntegrationTolerance: 1e-4,\n maxIntegrationDepth: 12,\n integrationTolerance: 1e-7,\n maxSearchDepth: 12,\n maxRootsPerCircle: 2\n });\n\n if (true) {\n for (let i = 0; i < roots.length; i++) {\n complexNewtonRaphson(roots[i], dispersion, roots[i], {\n tolerance: 1e-5,\n multiplicity: 2\n });\n }\n }\n roots = roots.filter(([a, b]) => !isNaN(a) && !isNaN(b));\n roots.sort(([ar, ai], [br, bi]) => ar - br);\n const t1 = performance.now();\n mutable rootsTimeMS = t1 - t0;\n return roots;\n }\n return {\n symmetric: computeRoots(\"symmetric\"),\n antisymmetric: computeRoots(\"Antisymmetric\")\n };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1200,"value":"roots = computeRoots(f, { r0: radius })","pinned":false,"hidden":false,"mode":"ojs"},{"id":1208,"value":"rootSeeds = computeRoots(f0, { r0: radius })","pinned":false,"hidden":false,"mode":"ojs"},{"id":1496,"value":"$39","pinned":false,"hidden":false,"mode":"ojs"},{"id":1111,"value":"symmetricRootsLoci = traceRoots(rootSeeds.symmetric, \"symmetric\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":960,"value":"symmetricLociData = new Float32Array(symmetricRootsLoci.flat().flat())","pinned":false,"hidden":false,"mode":"ojs"},{"id":932,"value":"symmetricLociBuffer = {\n const buf = regl.buffer(symmetricLociData);\n invalidation.then(() => buf.destroy());\n return buf;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1119,"value":"antisymmetricRootsLoci = traceRoots(rootSeeds.antisymmetric, \"antisymmetric\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1121,"value":"antisymmetricLociData = new Float32Array(antisymmetricRootsLoci.flat().flat())","pinned":false,"hidden":false,"mode":"ojs"},{"id":1124,"value":"antisymmetricLociBuffer = {\n const buf = regl.buffer(antisymmetricLociData);\n invalidation.then(() => buf.destroy());\n return buf;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2151,"value":"model = mat4fromScaling(mat4create(), [0.1, 0.8, 0.1])","pinned":false,"hidden":false,"mode":"ojs"},{"id":2163,"value":"modelViewProjection = mat4create()","pinned":false,"hidden":false,"mode":"ojs"},{"id":372,"value":"plot = PlotContext({\n yrange: [-20, 20],\n width: viewport.width,\n height: Math.max(400, Math.floor(viewport.width * 0.8))\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":896,"value":"camera = {\n const camera = createReglCamera(regl, {\n distance: 3,\n near: 0.01,\n far: 100,\n phi: 0.1,\n theta: 0.6,\n center: [0, 0.75, 0]\n });\n camera.interactions = createInteractions(camera);\n return camera;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":929,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs"},{"id":2202,"value":"drawPoints = regl({\n vert: `\nprecision highp float;\nattribute vec3 position;\nuniform mat4 matrix;\nuniform float pointSize;\nvoid main () {\n gl_Position = matrix * vec4(position, 1);\n gl_PointSize = pointSize;\n}`,\n frag: `\nprecision highp float;\nuniform vec3 color;\nuniform float pointSize;\nvoid main () {\nfloat r = length(gl_PointCoord.xy - 0.5) * pointSize * 2.0;\nfloat alpha = smoothstep(pointSize, pointSize - 1.0, r);\n gl_FragColor = vec4(color, alpha);\n\n}`,\n attributes: {\n position: regl.prop(\"position\")\n },\n uniforms: {\n color: regl.prop(\"color\"),\n matrix: (ctx, props) => props.matrix || mat4create(),\n pointSize: (ctx, props) => ctx.pixelRatio * (props.pointSize || 10)\n },\n blend: {\n enable: true,\n func: {\n srcRGB: \"src alpha\",\n srcAlpha: 1,\n dstRGB: \"one minus src alpha\",\n dstAlpha: 1\n },\n equation: {\n rgb: \"add\",\n alpha: \"add\"\n }\n },\n count: regl.prop(\"count\"),\n primitive: \"point\",\n depth: {\n enable: false\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":909,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs"},{"id":907,"value":"drawLoop = {\n camera.taint();\n let frame = regl.frame(({ tick }) => {\n try {\n renderFrame({ tick });\n } catch (e) {\n console.error(e);\n frame && frame.cancel();\n frame = null;\n }\n });\n invalidation.then(() => {\n frame && frame.cancel();\n frame = null;\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":912,"value":"$3c","pinned":false,"hidden":false,"mode":"ojs"},{"id":21,"value":"## Imports","pinned":false,"hidden":false,"mode":"md"},{"id":1813,"value":"import { ode45 } from \"@rreusser/integration\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":900,"value":"reglLines = require(\"https://unpkg.com/regl-gpu-lines@2.2.0/dist/regl-gpu-lines.js\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":893,"value":"import { createReglCamera, createInteractions } from \"@rreusser/regl-tools\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":879,"value":"import { reglCanvas } from \"@rreusser/regl-canvas\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":3180,"value":"figureWidth = width","pinned":false,"hidden":false,"mode":"ojs"},{"id":882,"value":"viewport = ({\n width: figureWidth,\n height: Math.floor(Math.max(400, figureWidth * 0.7)),\n pixelRatio: devicePixelRatio\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":362,"value":"import {\n PlotContext,\n DomainColoringLayer,\n ComplexVariable,\n domainColoringConfig\n} from \"@rreusser/complex-function-plotter\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":317,"value":"import { adaptiveSimpson, vectorAdaptiveSimpson } from \"@rreusser/integration\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":321,"value":"findRoots = (await import(\"https://cdn.skypack.dev/durand-kerner\")).default","pinned":false,"hidden":false,"mode":"ojs"},{"id":158,"value":"import {\n ccreateD,\n ccreateS,\n csetValuesD,\n csetValuesS,\n csqrD,\n csubDS,\n csubSD,\n cmulSS,\n cmulDS,\n cmulSD,\n csqrtD,\n csqrtAltD,\n caddDD,\n caddSD,\n cmulDD,\n csubDD,\n cscaleD,\n cscaleS,\n csqrtS,\n csqrS,\n cinvS,\n cdivSS,\n csinhcoshD\n} from \"@rreusser/complex-auto-differentiation\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":9,"value":"import { fmt as createFmt, sci } from \"@rreusser/fmt\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":2090,"value":"import { vec3transformMat4 } from \"@rreusser/gl-vec3\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":914,"value":"import {\n mat4fromScaling,\n mat4scale,\n mat4create,\n mat4fromTranslation,\n mat4translate,\n mat4invert,\n mat4rotateX,\n mat4rotateY,\n mat4rotateZ,\n mat4multiply,\n mat4rotate,\n mat4copy,\n mat4equals\n} from \"@rreusser/gl-mat4\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":2570,"value":"import { glslComplex } from \"@rreusser/glsl-complex\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":917,"value":"import { mat3normalFromMat4, mat3create } from \"@rreusser/gl-mat3\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1549,"value":"import { vec2normalize, vec2length, vec2scale } from \"@rreusser/gl-vec2\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":85,"value":"fmt = createFmt(\".3~f\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":72,"value":"scifmt = createFmt(sci(\".3f\"))","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W3d","stdlib":"1","titling":"auto"}