このシミュレーションは、皿の水面を細かなマスに分け、GPUで全部のマスを同時に計算している。1回の刻みは2段に分かれる。まずマスごとに水面の曲がり具合を求め、次にそれを使って、上下の速さと高さを少しだけ進める。これを1コマに数十回くり返す。周波数や水深、液体を変えるたびに、まずCPUで分散の式としきい値を解き、GPUに渡す係数をそろえる。
液体と容器の数字
長さはmm、時間は秒で計算する。表面張力σ ρ σ / ρ
定数と単位(src/sims/faraday-waves.ts の 26〜39 行) const G = 9810 ; // 重力加速度 mm/s² const D = 100 ; // 容器の直径(四角は 1 辺)mm const LIQ: Record<string, { label: string; rho: number; sigma: number; nu: number }> = { water: { label: "水" , rho: 998.2 , sigma: 0.0727 , nu: 1.0 }, // 20 °C。σ N/m・ρ kg/m³・ν mm²/s oil: { label: "シリコーン油" , rho: 950 , sigma: 0.0206 , nu: 20 }, // ジメチルシリコーンオイル 20 mm²/s(25 °C) }; const NLAM = [8 , 9 , 11 ], NMIN = [80 , 96 , 128 ], NMAX = [224 , 288 , 384 ]; // 1 波長のマスの数・格子の大きさ(画質 0・1・2) const CELLS = [1.2e6 , 2.5e6 , 5e6 ]; // 1 コマの計算量の上限(マスの数 × 刻みの数) const MAXSTEP = 64 ;const DRIFT = 1.5 ; // ストロボで見た水面の見かけの上下(回/秒) const PW = 256 ; // 断面を読む点の数 const SB = 0.35 ; // 波がくだけ始める傾き const FING = { d: 0.6 , a: 3 }; // 指が押し下げる深さ・広がり mm const WALL = 2.5 ; // 壁ぎわの帯の減衰 ÷ 粘りの減衰 FINGは指で水面を押す強さ、WALLは壁ぎわの摩擦の強さである。どちらも本物の値を測ったものではなく、見え方を決める目安だ。
分散の式から、波長と刻みを決める
solveKは、本文の分散の式ω 0 2 = ( g k + σ k 3 / ρ ) tanh k h ω 0 = π f k k k
分散としきい値(src/sims/faraday-waves.ts の 41〜115 行) // Phys: つまみから作る計算の数。f 周波数・wc 波の角振動数 πf・k と lam 波数と波長・Gc しきい値・N 格子のマス数・dx マスの幅・ // np 波の 1 周期の刻みの数・dt 刻み・hg と hs 重力と表面張力の項の係数・U 1 刻みの減衰・Z 容器の上下の片振幅(mm) interface Phys { f: number; om: number; wc: number; k: number; lam: number; h: number; s: number; nu: number; Gam: number; gam: number; Gc: number; Gf: number; m: number; N: number; dx: number; np: number; dt: number; hg: number; hs: number; ct: number; U: number; Fd: number; rlin: number; rnl: number; loc: number; brk: number; wall: number; wd: number; Z: number; liq: string; shape: string; } function solveK(f: number, h: number, s: number) { const w2 = (Math.PI * f) ** 2 ; let lo = 1e-5 , hi = 60 ; for (let i = 0 ; i < 90 ; i++) { const m = 0.5 * (lo + hi); if ((G * m + s * m * m * m) * Math.tanh(m * h) > w2) hi = m; else lo = m; } return 0.5 * (lo + hi); } /** 9 点のラプラシアンの大きさ × dx²(波数 (tx, ty)/dx の波に対して) */ const L9 = (tx: number, ty: number) => (20 - 8 * Math.cos(tx) - 8 * Math.cos(ty) - 4 * Math.cos(tx) * Math.cos(ty)) / 6 ; /** 1 つの波(k_c)の、この計算の刻みでのしきい値(フロケの方法: 揺れの 1 周期の写像の固有値が 1 を超える Γ) */ function threshold(P: Phys) { const half = P.np / 2 ; const rho = (Gam: number) => { let a1 = 1 , b1 = 0 , a2 = 0 , b2 = 1 ; // (η, v/ω) for (let n = 0 ; n < half; n++) { const W2 = P.ct * P.wc * P.wc * (P.Gf * (1 - Gam * Math.cos(P.om * n * P.dt)) + (1 - P.Gf)); b1 = P.U * (P.Fd * b1 - (P.dt * W2 * a1) / P.wc); a1 += P.dt * P.wc * b1; b2 = P.U * (P.Fd * b2 - (P.dt * W2 * a2) / P.wc); a2 += P.dt * P.wc * b2; } const tr = a1 + b2, det = a1 * b2 - a2 * b1, disc = tr * tr - 4 * det; if (disc >= 0 ) { const q = Math.sqrt(disc); return Math.max(Math.abs(tr + q), Math.abs(tr - q)) / 2 ; } return Math.sqrt(Math.abs(det)); }; if (rho(60 ) <= 1 ) return Infinity; let lo = 0 , hi = 60 ; for (let i = 0 ; i < 36 ; i++) { const m = 0.5 * (lo + hi); if (rho(m) > 1 ) hi = m; else lo = m; } return hi; } function physics(api: Api): Phys { const f = clamp(api.num("freq" ) || 40 , 5 , 400 ), Gam = Math.max(0 , api.num("amp" ) || 0 ), h = clamp(api.num("depth" ) || 10 , 1 , 50 ); const liq = LIQ[api.str("liquid" )] ? api.str("liquid" ) : "water" , L = LIQ[liq]; const s = (L.sigma / L.rho) * 1e9 , nu = L.nu; // σ/ρ mm³/s² const wc = Math.PI * f, om = 2 * wc, k = solveK(f, h, s), lam = TAU / k, q = api.quality; const N = clamp(Math.round(((D / lam) * NLAM[q] + 4 ) / 16 ) * 16 , NMIN[q], NMAX[q]); const dx = D / (N - 4 ); // 格子の上で、波数 k の波に効く Δ と Δ² の大きさ(向きで平均) const th = k * dx; let K1 = 0 , K2 = 0 , q5 = 0 ; const nA = 24 ; for (let i = 0 ; i < nA; i++) { const p = (i / nA) * TAU, tx = th * Math.cos(p), ty = th * Math.sin(p), a = L9(tx, ty) / (dx * dx); K1 += a / nA; K2 += (a * a) / nA; q5 += (4 * Math.sin(tx / 2 ) ** 2 + 4 * Math.sin(ty / 2 ) ** 2 ) / nA; } const he = Math.tanh(k * h) / k, hg = (he * k * k) / K1, hs = (he * k ** 4 ) / K2; const Kmax = 6 / (dx * dx); // 壁のある格子のラプラシアンの大きさの上限(余裕をみた値) const wmax = Math.sqrt(hg * G * (1 + Math.max(Gam, 1 )) * Kmax + hs * s * Kmax * Kmax); let np = Math.max(16 , Math.ceil((TAU * wmax) / (1.8 * wc))); np += np % 2 ; const dt = 2 / f / np; const ct = ((2 / (wc * dt)) * Math.sin((wc * dt) / 2 )) ** 2 ; // 時間の刻みで振動が速くなる分の補正 const x2 = 2 * k * h; const gam = 2 * nu * k * k + (x2 > 40 ? 0 : (k * Math.sqrt((nu * wc) / 2 )) / Math.sinh(x2)); // 粘り + 底の境界層 const rlin = Math.min(0.1 , (4 * nu * dt) / (dx * dx)); const Fd = 1 - rlin * q5; const U = Math.min(1 , Math.exp(-2 * gam * dt) / Fd); // k_c での減衰がちょうど 2γ になるように残りを一様に掛ける const delta = gam / wc, m = clamp((delta - 0.03 ) / 0.12 , 0 , 1 ) * 0.7 ; // 粘る液体ほど縞になりやすい形を混ぜる const bdiv = (470 * gam) / (k * k); // 傾きの 2 乗に比例する減衰の強さ(Γ = 2Γc で傾き 0.15 ほどで飽和) // 壁の境界層の減衰: 壁から λ/4 ほどの帯に、粘りの減衰の WALL 倍を置く。容器全体で平らにならした分をしきい値の計算に足す const wd = lam / 4 , gw = WALL * gam, gwAvg = gw * ((2 * wd) / (D / 2 )); const P: Phys = { f, om, wc, k, lam, h, s, nu, Gam, gam, Gc: 0 , Gf: G / (G + s * k * k), m, N, dx, np, dt, hg, hs, ct, U, Fd, rlin, rnl: (dt * bdiv * (1 - m)) / (dx * dx), loc: dt * 3 * 470 * gam * m, brk: dt * 60 * gam, wall: 2 * gw * dt, wd: wd / dx, Z: (Gam * G) / (om * om), liq, shape: api.str("shape" ), }; P.Gc = threshold({ ...P, U: P.U * Math.exp(-2 * gwAvg * dt) }); return P; } /** 飽和したときの傾きの目安 */ const satSlope = (Gam: number, Gc: number) => Math.sqrt(0.0225 * Math.max(0 , Gam / Gc - 1 ));GPUの上ではtanh k h k tanh k h h e k 2 Δ hg・hsは、マス目の上で計算したときに、ちょうどk ω 0
1刻みの時間dtは、いちばん細かい波でも計算が壊れない長さから決める。表面張力の効く細かい波ほど速く上下するので、波の1周期を数十から数百の刻みに分けることになる。
thresholdは、本文のしきい値Γ c Γ
GPUで水面の式を解く
FS_AとFS_Bは、マスごとに同時に動く小さなプログラム（シェーダー）である。FS_Aが高さη Δ η FS_Bがそれを使って速さv
水面の式(src/sims/faraday-waves.ts の 117〜196 行) // GPU の小さなプログラム(シェーダー)。1 つのマスごとに同時に動く。r = 高さ η(mm)、g = 上下の速さ v(mm/s) const HEAD = ` #version 300 es precision highp float; precision highp int; precision highp sampler2D; uniform float uN; uniform float uR; uniform float uShape; out vec4 o; bool inside(ivec2 c) { int n = int(uN); if (c.x < 0 || c.y < 0 || c.x >= n || c.y >= n) return false; vec2 p = vec2(c) + 0.5 - uN * 0.5; return uShape < 0.5 ? dot(p, p) < uR * uR : max(abs(p.x), abs(p.y)) < uR; } float hash(vec2 p, float s) { vec3 q = fract(vec3(p.xyx) * 0.1031 + s * 0.0137); q += dot(q, q.yzx + 33.33); return fract((q.x + q.y) * q.z); } ` ;// 1. Δη と傾きの 2 乗 const FS_A = ` uniform sampler2D uS; uniform float uIdx2; uniform float uIdx; float at(ivec2 q, float c0) { return inside(q) ? texelFetch(uS, q, 0).r : c0; } void main() { ivec2 c = ivec2(gl_FragCoord.xy); if (!inside(c)) { o = vec4(0.0); return; } float e0 = texelFetch(uS, c, 0).r; float l = at(c + ivec2(-1, 0), e0), r = at(c + ivec2(1, 0), e0), d = at(c + ivec2(0, -1), e0), u = at(c + ivec2(0, 1), e0); float ld = at(c + ivec2(-1, -1), e0), rd = at(c + ivec2(1, -1), e0), lu = at(c + ivec2(-1, 1), e0), ru = at(c + ivec2(1, 1), e0); float lap = (4.0 * (l + r + d + u) + (ld + rd + lu + ru) - 20.0 * e0) / 6.0 * uIdx2; // 傾き: ふつうは両隣の差 ÷ 2。壁ぎわで片方が外なら、内側との差をそのまま使う(半分に小さく見積もらない) float hx = inside(c + ivec2(-1, 0)) && inside(c + ivec2(1, 0)) ? 0.5 : 1.0; float hy = inside(c + ivec2(0, -1)) && inside(c + ivec2(0, 1)) ? 0.5 : 1.0; float gx = hx * (r - l) * uIdx, gy = hy * (u - d) * uIdx; o = vec4(lap, gx * gx + gy * gy, 0.0, 0.0); }` ;// 2. 1 刻み進める const FS_B = ` uniform sampler2D uS; uniform sampler2D uT; uniform float uIdx2, uDt, uGt, uHg, uHs, uSs, uCt, uU, uRlin, uRnl, uLoc, uBrk, uSb2, uNoise, uSeed, uWall, uWd; uniform vec4 uF; float lt(ivec2 q, float c0) { return inside(q) ? texelFetch(uT, q, 0).r : c0; } void main() { ivec2 c = ivec2(gl_FragCoord.xy); vec4 s0 = texelFetch(uS, c, 0); ivec2 nb[4] = ivec2[4](ivec2(1, 0), ivec2(-1, 0), ivec2(0, 1), ivec2(0, -1)); if (!inside(c)) { float sum = 0.0, n = 0.0; for (int i = 0; i < 4; i++) { ivec2 q = c + nb[i]; if (inside(q)) { sum += texelFetch(uS, q, 0).r; n += 1.0; } } o = vec4(n > 0.0 ? sum / n : 0.0, 0.0, 0.0, 0.0); return; } vec4 t0 = texelFetch(uT, c, 0); float a0 = t0.r; float l = lt(c + ivec2(-1, 0), a0), r = lt(c + ivec2(1, 0), a0), d = lt(c + ivec2(0, -1), a0), u = lt(c + ivec2(0, 1), a0); float ld = lt(c + ivec2(-1, -1), a0), rd = lt(c + ivec2(1, -1), a0), lu = lt(c + ivec2(-1, 1), a0), ru = lt(c + ivec2(1, 1), a0); float lap2 = (4.0 * (l + r + d + u) + (ld + rd + lu + ru) - 20.0 * a0) / 6.0 * uIdx2; float force = uCt * (uHg * uGt * a0 - uHs * uSs * lap2); if (uF.w != 0.0) { vec2 q = vec2(c) + 0.5 - uF.xy; float a2 = uF.z * uF.z, r2 = dot(q, q); force += uF.w * (4.0 * r2 / (a2 * a2) - 4.0 / a2) * exp(-r2 / a2) * uIdx2; } float v = s0.g; float vn = v + uDt * force; for (int i = 0; i < 4; i++) { ivec2 q = c + nb[i]; if (!inside(q)) continue; float vq = texelFetch(uS, q, 0).g, s2q = texelFetch(uT, q, 0).g; vn += min(uRlin + uRnl * 0.5 * (t0.g + s2q), 0.12) * (vq - v); } vn += uNoise * (hash(vec2(c), uSeed) - 0.5); float s2 = t0.g; // 壁ぎわの減衰(壁の境界層と接触線): 壁からの距離 dw(マス)で弱まる帯 vec2 pc = vec2(c) + 0.5 - uN * 0.5; float dw = uR - (uShape < 0.5 ? length(pc) : max(abs(pc.x), abs(pc.y))); vn *= uU / (1.0 + uLoc * s2 + uBrk * max(0.0, s2 - uSb2) / uSb2 + uWall * exp(-dw / uWd)); o = vec4(s0.r + uDt * vn, vn, 0.0, 0.0); }` ;FS_Bのforceが、本文の2つの式をつないだ部分だ。uGtに入る揺れる容器の中の重力g ( 1 − Γ cos 2 π f t ) uHg * uGt * a0だけに掛かり、表面張力の項uHs * uSs * lap2には掛からない。しきい値の式のかっこの中の補正は、この違いから生まれる。
速さには、隣のマスとの差をならす減衰をかける。粘りによる分（uRlin）と、傾きが大きい所ほど強くなる分（uRnl）で、後者が模様の高さを頭打ちにする。壁の外のマスは中のマスと同じ値とみなすので、壁の所で水面の傾きは0になる。
実時間とストロボ
水面は1秒に数十回上下するので、そのまま描くとちらついて見えない。planは、計算は実時間で進めながら、波のちょうど1周期分ずつまとめて進めた所だけを画面に出す。1周期後の水面は同じ形に戻るので、少しずつずらして見せると、ゆっくり上下するように見える。ストロボの光で照らすのと同じ仕組みだ。
時間の進め方(src/sims/faraday-waves.ts の 340〜391 行) /** * 1 コマで進める刻みの数を決める。 * 速いとき(実時間)は、波の 1 周期ぶんをまとめて進めて同じ形に戻った所だけを見せる(ストロボ)。見かけの上下は 1 秒に DRIFT 回。 * 指で触った直後は、まとめて進める分をやめて、見かけ 1.5 回/秒のスローで見せる。 * 1 コマの計算量(マスの数 × 刻み)が上限を超えるときは、実時間より遅く進める。 */ function plan(api: Api, s: S) { const P = s.P, dt = api.dt, np = P.np; const budget = Math.max(1 , Math.min(MAXSTEP, Math.floor(CELLS[api.quality] / (P.N * P.N)))); const rate = DRIFT + (P.f / 2 - DRIFT) * s.fast; // 1 秒に進める波の周期の数 if (rate <= 6 || np > budget) { s.strobe = false ; s.acc += rate * dt * np; const n = Math.min(budget, Math.floor(s.acc)); s.acc -= n; if (s.acc > np) s.acc = 0 ; return n; } s.strobe = true ; s.acc += DRIFT * dt * np; const nd = Math.min(budget, Math.floor(s.acc)); s.acc -= nd; if (s.acc > np) s.acc = 0 ; s.accW += (rate - DRIFT) * dt; let nw = Math.floor(s.accW); s.accW -= nw; if (nd + nw * np > budget) { nw = Math.max(0 , Math.floor((budget - nd) / np)); s.accW = Math.min(s.accW, 0.99 ); } return nd + nw * np; } /** steps 回、刻み dt で進める。1 回 = Δη を求める(FS_A)→ 速さと高さを更新(FS_B)の 2 段 */ function run(api: Api, s: S, steps: number) { const g = s.g!, gl = g.gl, P = s.P, N = P.N; const F = s.finger, fOn = !!F && F.on; const toCell = (mm: number) => mm / P.dx + N / 2 ; const pA = g.pA, pB = g.pB; pA.use(); grid(pA, P); pA.set("uIdx2" , 1 / (P.dx * P.dx)); pA.set("uIdx" , 1 / P.dx); pB.use(); grid(pB, P); pB.set("uIdx2" , 1 / (P.dx * P.dx)); pB.set("uDt" , P.dt); pB.set("uHg" , P.hg); pB.set("uHs" , P.hs); pB.set("uSs" , P.s); pB.set("uCt" , P.ct); pB.set("uU" , P.U); pB.set("uRlin" , P.rlin); pB.set("uRnl" , P.rnl); pB.set("uLoc" , P.loc); pB.set("uBrk" , P.brk); pB.set("uSb2" , SB * SB); pB.set("uWall" , P.wall); pB.set("uWd" , P.wd); pB.set("uNoise" , (2e-4 * P.wc * P.wc * P.dt) / P.k); // 床や空気からいつも来ている、ごく小さな揺れ for (let n = 0 ; n < steps; n++) { gl.bindFramebuffer(gl.FRAMEBUFFER, g.Tfb); gl.viewport(0 , 0 , N, N); pA.use(); pA.tex("uS" , 0 , g.S.read); g.quad.draw(); g.S.bindWrite(); gl.viewport(0 , 0 , N, N); pB.use(); pB.tex("uS" , 0 , g.S.read); pB.tex("uT" , 1 , g.T); pB.set("uGt" , G * (1 - P.Gam * Math.cos(P.om * s.t))); // 揺れる容器の中の重力 g(t) pB.set("uSeed" , (s.seedN = (s.seedN + 1 ) % 997 )); if (fOn) { const k = (n + 1 ) / steps, x = F!.px + (F!.x - F!.px) * k, y = F!.py + (F!.y - F!.py) * k; pB.set("uF" , toCell(x), toCell(y), FING.a / P.dx, P.ct * P.hg * G * FING.d); } else pB.set("uF" , 0 , 0 , 1 , 0 ); g.quad.draw(); g.S.swap(); s.t += P.dt; } if (F) { F.px = F.x; F.py = F.y; } } 1コマに進める刻みの数には上限があり、マスの数が多い高い周波数では、実時間より遅く進めることがある。そのときは、画面の下に「実際の何倍」と出る。
指・音・マイク
指は、水面を押す圧力として式に入れる。置いたままならくぼみになり、動かすと引き波が立ち、離すとくぼみが戻って波紋になる。
指で触る(src/sims/faraday-waves.ts の 556〜595 行) /** 画面の点を、容器の中心からの mm に直す(容器の外なら null) */ function fingerAt(api: Api, s: S): [number, number] | null { const L = layout(api), scale = (L.S * 0.96 ) / D; const x = (api.pointer.x - L.cx) / scale, y = -(api.pointer.y - L.cy) / scale, lim = D / 2 - FING.a; const inside = s.P.shape === "square" ? Math.max(Math.abs(x), Math.abs(y)) < D / 2 : Math.hypot(x, y) < D / 2 ; if (!inside) return null ; if (s.P.shape === "square" ) return [clamp(x, -lim, lim), clamp(y, -lim, lim)]; const r = Math.hypot(x, y); return r > lim ? [(x / r) * lim, (y / r) * lim] : [x, y]; } /** * 指は、水面を押す圧力 p(深さ 0.6 mm・広がり 3 mm の山形)として式に入れる(FS_B の uF)。 * 置いたまま = くぼみ、動かす = 引き波(かき混ぜる)、離す = くぼみが戻って波紋になる。触っている間と直後はスローで見せる */ function touch(type : "down" | "move" | "up" | "cancel" | "wheel" , api: Api, s: S) { if (type === "down" ) { const p = fingerAt(api, s); if (!p) return ; s.auto = 0 ; s.finger = { x: p[0 ], y: p[1 ], px: p[0 ], py: p[1 ], on: true , up: false , t0: api.t, auto: false }; s.slowT = 2 ; api.emit("touch" , { kind: "water" }); } else if (type === "move" ) { const F = s.finger; if (!F || !F.on || F.auto || !api.pointer.down) return ; const p = fingerAt(api, s); if (!p) return ; F.x = p[0 ]; F.y = p[1 ]; } else if (type === "up" || type === "cancel" ) { if (s.finger && !s.finger.auto) s.finger.up = true ; } } /** 場面「かき混ぜたあと」: 見えない指で、半径 24 mm の円を 1.3 秒で 1 周なぞる */ function autoStir(api: Api, s: S, dt: number) { if (s.auto <= 0 ) return ; s.auto -= dt; s.autoA += (dt * TAU) / 1.3 ; const x = 24 * Math.cos(s.autoA), y = 24 * Math.sin(s.autoA); if (!s.finger || !s.finger.auto) s.finger = { x, y, px: x, py: y, on: true , up: false , t0: api.t, auto: true }; else { s.finger.x = x; s.finger.y = y; } if (s.auto <= 0 ) s.finger.up = true ; } 音は、揺らす周波数の正弦波を小さな音量で鳴らすだけだ。マイクは、押して許可されたときだけ開く。音を周波数ごとの強さに分け、10〜240 Hzでいちばん強い所を探して、その周波数と大きさでつまみを動かす。音はどこにも送らず、止めるとすぐに閉じる。
音とマイク(src/sims/faraday-waves.ts の 467〜542 行) /** 音の大きさ: 振幅に合わせて少しだけ大きくする(小さな音量のまま) */ const volOf = (api: Api) => clamp(0.02 + 0.025 * api.num("amp" ), 0.02 , 0.07 );function soundOn(api: Api, s: S) { if (s.osc) return ; const ac = api.audio(); if (!ac) { api.note("音を出せませんでした。もう一度「音を鳴らす」を押すと鳴ります" ); api.set("sound" , false ); return ; } const osc = ac.createOscillator(), gain = ac.createGain(); osc.type = "sine" ; osc.frequency.value = s.P.f; gain.gain.value = 0 ; osc.connect(gain); gain.connect(ac.destination); osc.start(); gain.gain.setTargetAtTime(volOf(api), ac.currentTime, 0.05 ); s.osc = osc; s.gain = gain; if (s.P.f < 100 && !s.lowNoted) { s.lowNoted = true ; api.note("100 Hz より低い音は、スマホやノートパソコンのスピーカーでは小さく聞こえます" , 3200 ); } } function soundOff(s: S) { const o = s.osc, gn = s.gain; s.osc = null ; s.gain = null ; if (!o || !gn) return ; try { gn.gain.setTargetAtTime(0 , o.context.currentTime, 0.03 ); } catch { /* 閉じた後 */ } setTimeout(() => { try { o.stop(); o.disconnect(); gn.disconnect(); } catch { /* 止めた後 */ } }, 250 ); } function soundUpdate(api: Api, s: S) { if (!s.osc || !s.gain) return ; const t = s.osc.context.currentTime; s.osc.frequency.setTargetAtTime(s.P.f, t, 0.03 ); s.gain.gain.setTargetAtTime(volOf(api), t, 0.08 ); } /** マイク: 押して許可したときだけ開く。音は端末の中で周波数ごとの強さにするだけで、どこへも送らない。止めたらすぐ閉じる */ function micOn(api: Api, s: S) { if (s.mic || s.micPending) return ; const fail = (why: string) => { s.micPending = false ; s.micErr = why; if (s.alive) { api.note(why, 4200 ); api.set("mic" , false ); } }; const md = (globalThis.navigator as Navigator | undefined )?.mediaDevices; if (!md || typeof md.getUserMedia !== "function" ) { fail("このブラウザではマイクを使えないため、マイクで揺らせません" ); return ; } s.micPending = true ; s.micErr = "" ; md.getUserMedia({ audio: { echoCancellation: false , noiseSuppression: false , autoGainControl: false }, video: false }).then((stream) => { s.micPending = false ; if (!s.alive || !api.bool("mic" )) { stream.getTracks().forEach((t) => t.stop()); return ; } let ac = api.audio(), own: AudioContext | null = null ; if (!ac) { try { own = ac = new AudioContext(); } catch { ac = null ; } } if (!ac) { stream.getTracks().forEach((t) => t.stop()); fail("音の読み取りを始められなかったため、マイクで揺らせません" ); return ; } if (ac.state === "suspended" ) ac.resume().catch (() => {}); const src = ac.createMediaStreamSource(stream), an = ac.createAnalyser(); an.fftSize = 8192 ; an.smoothingTimeConstant = 0.5 ; src.connect(an); s.mic = { stream, src, an, buf: new Float32Array(an.frequencyBinCount), f: s.P.f, db: -120 , tSet: 0 , own }; s.micErr = "" ; api.note("マイクの音で揺らしています。10〜240 Hz でいちばん強い音の高さと大きさを使い、音は端末の外へ送りません" , 4200 ); }).catch ((e: unknown) => { const n = (e as { name?: string })?.name; fail(n === "NotAllowedError" || n === "SecurityError" ? "マイクの使用が許可されなかったため、マイクで揺らせません" : "マイクが見つからないため、マイクで揺らせません" ); }); } function micOff(s: S) { const m = s.mic; s.mic = null ; if (!m) return ; m.stream.getTracks().forEach((t) => t.stop()); try { m.src.disconnect(); m.an.disconnect(); } catch { /* 切った後 */ } m.own?.close().catch (() => {}); } /** マイクの音を周波数ごとの強さに分け、10〜240 Hz でいちばん強い所の周波数と大きさで、つまみを動かす */ function micUpdate(api: Api, s: S, dt: number) { const m = s.mic; if (!m) return ; m.an.getFloatFrequencyData(m.buf); const bw = m.an.context.sampleRate / m.an.fftSize; let best = -Infinity, bi = -1 ; for (let i = Math.max(1 , Math.ceil(10 / bw)); i <= Math.floor(240 / bw) && i < m.buf.length - 1 ; i++) if (m.buf[i] > best) { best = m.buf[i]; bi = i; } if (bi < 0 || !Number.isFinite(best)) return ; const y0 = m.buf[bi - 1 ], y2 = m.buf[bi + 1 ], den = y0 - 2 * best + y2; const fpk = (bi + (Number.isFinite(den) && Math.abs(den) > 1e-6 ? clamp((0.5 * (y0 - y2)) / den, -0.5 , 0.5 ) : 0 )) * bw; m.db = best; const level = clamp((best + 85 ) / 50 , 0 , 1 ); // −85 dB で 0、−35 dB で 1 if (level > 0.05 ) m.f += (fpk - m.f) * Math.min(1 , dt / 0.25 ); m.tSet += dt; if (m.tSet < 0.15 ) return ; m.tSet = 0 ; const fr = Math.round(clamp(m.f, 10 , 240 )), gm = Math.round(Math.min(6 , s.P.Gc * 4 * level) * 100 ) / 100 ; if (fr !== api.num("freq" ) || gm !== api.num("amp" )) { api.set("freq" , fr); api.set("amp" , gm); update(api, s, true ); } } 描く
上から見た水面は、マスの高さを双三次の補間でなめらかにつなぎ、傾きから光の当たり方を求めて網点にしている。断面のグラフはGPUから直径に沿った256点の高さを読み出したもので、読み出しを待たずに次のコマで受け取るので、画面は止まらない。
描画(src/sims/faraday-waves.ts の 597〜803 行) // 上から見た水面: 高さを双三次の補間でなめらかにし、傾きから光の当たり方(陰影と照り返し)を出して、白と黒の網点にする const FS_SHOW = ` uniform sampler2D uS; uniform float uDx, uD, uCell, uBase, uGain, uSpec, uShin; uniform vec3 uL; uniform vec3 uC; ${DITHER} float eta(ivec2 q) { return texelFetch(uS, clamp(q, ivec2(0), ivec2(int(uN) - 1)), 0).r; } void cr(float t, out vec4 w, out vec4 dw) { float t2 = t * t, t3 = t2 * t; w = vec4(-0.5 * t3 + t2 - 0.5 * t, 1.5 * t3 - 2.5 * t2 + 1.0, -1.5 * t3 + 2.0 * t2 + 0.5 * t, 0.5 * t3 - 0.5 * t2); dw = vec4(-1.5 * t2 + 2.0 * t - 0.5, 4.5 * t2 - 5.0 * t, -4.5 * t2 + 4.0 * t + 0.5, 1.5 * t2 - t); } vec4 row(ivec2 i0, int j) { return vec4(eta(i0 + ivec2(-1, j)), eta(i0 + ivec2(0, j)), eta(i0 + ivec2(1, j)), eta(i0 + ivec2(2, j))); } void main() { vec2 fc = gl_FragCoord.xy; vec2 mm = (fc - uC.xy) / uC.z; float rr = uShape < 0.5 ? length(mm) : max(abs(mm.x), abs(mm.y)); if (rr > uD * 0.5) { o = vec4(uPaper, 1.0); return; } vec2 g = mm / uDx + uN * 0.5 - 0.5; ivec2 i0 = ivec2(floor(g)); vec2 f = g - floor(g); vec4 wx, dwx, wy, dwy; cr(f.x, wx, dwx); cr(f.y, wy, dwy); vec4 r0 = row(i0, -1), r1 = row(i0, 0), r2 = row(i0, 1), r3 = row(i0, 2); vec4 a = vec4(dot(wx, r0), dot(wx, r1), dot(wx, r2), dot(wx, r3)); vec4 b = vec4(dot(dwx, r0), dot(dwx, r1), dot(dwx, r2), dot(dwx, r3)); float ex = dot(wy, b) / uDx, ey = dot(dwy, a) / uDx; vec3 n = normalize(vec3(-ex, -ey, 1.0)); float dif = dot(n, uL) - uL.z; vec3 H = normalize(uL + vec3(0.0, 0.0, 1.0)); float sp = max(0.0, pow(max(dot(n, H), 0.0), uShin) - pow(H.z, uShin)); float tone = uBase + uGain * dif + uSpec * sp; o = vec4(monoCell(clamp(tone, 0.0, 1.0), fc, uCell), 1.0); }` ; interface Box { x: number; y: number; w: number; h: number }interface Lay { land: boolean; k: number; cx: number; cy: number; S: number; cap: Box; sec: Box | null ; scope: Box; gauge: Box | null }/** * 配置は描画面の縦横から決める。横長なら左に容器・右にグラフ、縦長や正方形に近ければ上に容器・下にグラフ 2 つ。 * 文字の大きさの単位 k は api.u(1) を、小さな画面でも読める範囲に収めたもの */ function layout(api: Api): Lay { const W = api.w, H = api.h, k = clamp(api.u(1 ), 0.82 , 1.3 ), U = (n: number) => n * k; const pad = U(14 ), cap = U(28 ); if (W >= H * 1.1 ) { const S = Math.max(U(120 ), Math.min(H - 2 * pad - cap, W * 0.6 - pad, W - 4 * pad - U(240 ))); const pw = clamp(W - S - 4 * pad, U(200 ), Math.min(S * 1.05 , U(600 ))); const x0 = Math.max(pad, (W - (S + 2 * pad + pw)) / 2 ), y0 = Math.max(pad, (H - S - cap) / 2 ); const px = x0 + S + 2 * pad, ax = U(44 ), bh = S + cap, plotW = pw - ax; const base = { land: true , k, cx: x0 + S / 2 , cy: y0 + S / 2 , S, cap: { x: x0, y: y0 + S + cap * 0.72 , w: S, h: 0 } }; // 低い描画面(スマホで操作パネルを開いたときなど)では、振動波形だけを大きく出す if (bh < U(340 )) return { ...base, sec: null , gauge: null , scope: { x: px + ax, y: y0 + U(30 ), w: plotW, h: Math.max(U(60 ), bh - U(56 )) } }; const secH = bh * 0.27 , scH = bh * 0.33 , secY = y0 + U(34 ); return { ...base, sec: { x: px + ax, y: secY, w: plotW, h: secH }, scope: { x: px + ax, y: secY + secH + U(62 ), w: plotW, h: scH }, gauge: { x: px, y: y0 + bh - U(16 ), w: pw, h: 0 } }; } const want = U(230 ); const S = Math.max(U(120 ), Math.min(W - 2 * pad, H - 2 * pad - cap - want)); const rest = H - 2 * pad - S - cap, plotsH = Math.min(rest, Math.max(want, S * 0.7 )); const y0 = pad + (rest - plotsH) / 2 ; const top = y0 + S + cap + U(30 ), ax = U(38 ), gap = U(18 ); const bw = (W - 2 * pad - 2 * ax - gap) / 2 , bh = Math.max(U(60 ), plotsH - U(30 ) - U(66 )); const base = { land: false , k, cx: W / 2 , cy: y0 + S / 2 , S, cap: { x: W / 2 - S / 2 , y: y0 + S + cap * 0.72 , w: S, h: 0 } }; if (bh < U(80 )) return { ...base, sec: null , gauge: null , scope: { x: pad + ax, y: top, w: W - 2 * pad - ax, h: Math.max(U(40 ), plotsH - U(56 )) } }; return { ...base, sec: { x: pad + ax, y: top, w: bw, h: bh }, scope: { x: pad + 2 * ax + bw + gap, y: top, w: bw, h: bh }, gauge: { x: pad, y: y0 + S + cap + plotsH - U(18 ), w: W - 2 * pad, h: 0 } }; } const niceUp = (v: number) => { const e = Math.pow(10 , Math.floor(Math.log10(v))), m = v / e; return (m <= 1 ? 1 : m <= 2 ? 2 : m <= 2.5 ? 2.5 : m <= 5 ? 5 : 10 ) * e; };const decs = (v: number) => (v >= 1 ? 1 : v >= 0.1 ? 2 : 3 );/** 文字が幅に収まるように大きさを下げる */ function fitSize(api: Api, str: string, maxW: number, size: number) { api.ctx.font = ` 400 ${size}px sans-serif` ; const w = api.ctx.measureText(str).width; return w > maxW ? Math.max(size * 0.6 , (size * maxW) / w) : size; } function show(api: Api, s: S) { const g = s.g!, gl = g.gl, P = s.P, L = layout(api); toScreen(gl, api.glw, api.glh); const sx = api.glw / api.w, sy = api.glh / api.h, scale = (L.S * 0.96 ) / D; const p = g.pShow; p.use(); grid(p, P); p.tex("uS" , 0 , g.S.read); p.set("uDx" , P.dx); p.set("uD" , D); p.set("uCell" , Math.max(1 , Math.round(sx))); p.set("uC" , L.cx * sx, (api.h - L.cy) * sy, scale * sx); const l = [-0.42 , 0.5 , 0.76 ], ln = Math.hypot(l[0 ], l[1 ], l[2 ]); // 左上から差す光 p.set("uL" , l[0 ] / ln, l[1 ] / ln, l[2 ] / ln); p.set("uBase" , 0.13 ); p.set("uGain" , 6.5 ); p.set("uSpec" , 1.4 ); p.set("uShin" , 90 ); p.set("uInk" , ...api.inkRGB); p.set("uPaper" , ...api.paperRGB); g.quad.draw(); } function overlay(api: Api, s: S) { const ctx = api.ctx, P = s.P, L = layout(api), U = (n: number) => n * L.k; ctx.clearRect(0 , 0 , api.w, api.h); const scale = (L.S * 0.96 ) / D, R = (D / 2 ) * scale; // 容器の縁 if (P.shape === "square" ) api.rect(L.cx - R, L.cy - R, 2 * R, 2 * R, { w: U(2.2 ) }); else api.circle(L.cx, L.cy, R, { w: U(2.2 ) }); // 断面を取る直径と、波形を見る印の所 api.line(L.cx - R, L.cy, L.cx + R, L.cy, { w: U(0.8 ), dash: [U(3 ), U(5 )] }); const px = L.cx + ((s.probeI + 0.5 ) / PW - 0.5 ) * D * 0.99 * scale; if (s.have) { api.circle(px, L.cy, U(6 ), { w: U(1.4 ), fill: "paper" }); api.dot(px, L.cy, U(2 )); } // 指(触っている所)と、マウスの指す所 const F = s.finger, Pt = api.pointer; if (F && F.on) api.circle(L.cx + F.x * scale, L.cy - F.y * scale, FING.a * 1.6 * scale, { w: U(1.4 ), dash: [U(2 ), U(3 )] }); else if (Pt.inside && !Pt.down && Pt.kind === "mouse" && Math.hypot(Pt.x - L.cx, Pt.y - L.cy) < R) api.circle(Pt.x, Pt.y, FING.a * 1.6 * scale, { w: U(1 ), dash: [U(1.5 ), U(4 )] }); // 下の行: 物差しと、時間の流れ(マイクが使えないときは、その知らせ) const cy = L.cap.y, fs = U(11.5 ), lw = 10 * scale, bx = L.cap.x; api.line(bx, cy - fs * 0.35 , bx + lw, cy - fs * 0.35 , { w: U(1.4 ) }); api.line(bx, cy - fs * 0.35 - U(4 ), bx, cy - fs * 0.35 + U(4 ), { w: U(1.4 ) }); api.line(bx + lw, cy - fs * 0.35 - U(4 ), bx + lw, cy - fs * 0.35 + U(4 ), { w: U(1.4 ) }); api.text("1 cm" , bx + lw + U(6 ), cy, { size: fs, font: "mono" }); const msg = s.micErr && !s.mic ? s.micErr : speedText(s), mx = bx + lw + U(48 ), mw = L.cap.x + L.cap.w - mx; if (s.micErr && !s.mic) api.label(msg, mx, cy, { size: fitSize(api, msg, mw - U(8 ), fs) }); else api.text(msg, mx, cy, { size: fitSize(api, msg, mw, fs) }); if (L.sec) section(api, s, L.sec, U); scope(api, s, L.scope, U); if (L.gauge) gauge(api, s, L.gauge, U); } function speedText(s: S) { const r = s.speed; if (s.strobe) return r >= 0.93 ? "実時間・ストロボで見た水面" : ` 実際の ${fmt(r, 2)} 倍・ストロボで見た水面` ; if (r >= 0.93 ) return "実時間" ; return r > 0 ? ` スロー(実際の 1/${Math.max(2, Math.round(1 / r))})` : "止まっている" ; } /** 断面: 直径に沿った高さ。縦は横の何十倍にも拡大し、下の水は点で塗る */ function section(api: Api, s: S, r: Box, U: (n: number) => number) { const P = s.P, narrow = r.w < U(250 ); let mx = 0 ; for (let i = 0 ; i < PW; i++) mx = Math.max(mx, s.env[i]); const Y = niceUp(Math.max(0.02 , mx * 1.15 )); const X = (i: number) => r.x + ((i + 0.5 ) / PW) * r.w, Yp = (y: number) => r.y + r.h / 2 - (clamp(y, -Y * 1.3 , Y * 1.3 ) / Y) * (r.h / 2 ); api.text("直径に沿った断面" , r.x - U(40 ), r.y - U(16 ), { size: U(12.5 ), weight: 600 }); if (!narrow) api.text(` 高さは ${Math.round(r.h / 2 / Y / (r.w / D))} 倍に拡大` , r.x + r.w, r.y - U(16 ), { size: U(10.5 ), align: "right" }); const pts: [number, number][] = []; for (let i = 0 ; i < PW; i++) pts.push([X(i), Yp(s.have ? s.prof[i] : 0 )]); api.dots([...pts, [r.x + r.w, r.y + r.h], [r.x, r.y + r.h]], { gap: U(5 ), r: U(0.75 ) }); api.line(r.x, r.y + r.h, r.x + r.w, r.y + r.h, { w: U(1 ) }); api.line(r.x, Yp(0 ), r.x + r.w, Yp(0 ), { w: U(0.6 ), dash: [U(2 ), U(4 )] }); api.poly(pts, { w: U(1.8 ) }); const d = decs(Y); api.text(` +${fmt(Y, d)}` , r.x - U(5 ), Yp(Y) + U(4 ), { size: U(10 ), align: "right" , font: "mono" }); api.text(` −${fmt(Y, d)}` , r.x - U(5 ), Yp(-Y) + U(4 ), { size: U(10 ), align: "right" , font: "mono" }); api.text("mm" , r.x - U(5 ), Yp(0 ) + U(4 ), { size: U(10 ), align: "right" , font: "mono" }); api.text("−50" , r.x, r.y + r.h + U(14 ), { size: U(10 ), font: "mono" }); api.text("50 mm" , r.x + r.w, r.y + r.h + U(14 ), { size: U(10 ), align: "right" , font: "mono" }); // 波長 λ の物差し const lw = (P.lam / D) * r.w; if (lw < r.w * 0.55 ) { const y = r.y + U(4 ), x0 = r.x + r.w * (narrow ? 0.08 : 0.45 ) - lw / 2 * (narrow ? 0 : 1 ); api.line(x0, y, x0 + lw, y, { w: U(1.2 ) }); api.line(x0, y - U(4 ), x0, y + U(4 ), { w: U(1.2 ) }); api.line(x0 + lw, y - U(4 ), x0 + lw, y + U(4 ), { w: U(1.2 ) }); api.label(` λ ${fmt(P.lam, 1)} mm` , x0 + lw + U(6 ), y + U(4 ), { size: U(10.5 ), font: "mono" }); } if (s.have) { const x = X(s.probeI), y = Yp(s.prof[s.probeI]); api.circle(x, y, U(5 ), { w: U(1.3 ), fill: "paper" }); api.dot(x, y, U(1.8 )); } } /** 振動波形: 上 = 容器の上下(揺らす音の波形)、下 = 印の所の水面の上下。横は波の 1 周期(= 揺れの 2 周期)で折りたたむ */ function scope(api: Api, s: S, r: Box, U: (n: number) => number) { const P = s.P, narrow = r.w < U(250 ), h2 = r.h * 0.42 , gap = r.h - 2 * h2; const c1 = { x: r.x, y: r.y, w: r.w, h: h2 }, c2 = { x: r.x, y: r.y + h2 + gap, w: r.w, h: h2 }; const Xp = (ph: number) => r.x + ph * r.w; api.text("振動波形" , r.x - U(36 ), r.y - U(16 ), { size: U(12.5 ), weight: 600 }); if (!narrow) api.text("波の 1 周期 = 揺れの 2 周期" , r.x + r.w, r.y - U(16 ), { size: U(10.5 ), align: "right" }); const mid1 = c1.y + c1.h / 2 , A1 = c1.h * 0.42 ; api.line(c1.x, mid1, c1.x + c1.w, mid1, { w: U(0.6 ), dash: [U(2 ), U(4 )] }); const n = 120 , pts1: [number, number][] = []; for (let i = 0 ; i <= n; i++) { const ph = i / n; pts1.push([Xp(ph), mid1 - A1 * Math.cos(2 * TAU * ph) * (P.Gam > 0 ? 1 : 0 )]); } api.poly(pts1, { w: U(1.4 ) }); const zt = ` ±${fmt(P.Z * 1000, P.Z < 0.01 ? 1 : 0)} µm` ; api.label(narrow ? ` 容器 ${fmt(P.f, 0)} Hz` : ` 容器の上下 f ${fmt(P.f, 0)} Hz ${zt}` , c1.x, c1.y + U(2 ), { size: U(10.5 ) }); const mid2 = c2.y + c2.h / 2 ; api.line(c2.x, mid2, c2.x + c2.w, mid2, { w: U(0.6 ), dash: [U(2 ), U(4 )] }); const c0 = s.fit.ok ? s.fit.c : 0 ; let mx = 0 ; for (const p of s.samples) mx = Math.max(mx, Math.abs(p.y - c0)); if (s.fit.ok) mx = Math.max(mx, s.fit.amp); const Ym = niceUp(Math.max(0.01 , mx * 1.1 )), k2 = (c2.h * 0.46 ) / Ym, now = api.t; for (const p of s.samples) { const age = clamp((now - p.t) / 2.4 , 0 , 1 ); api.dot(Xp(p.ph), mid2 - (p.y - c0) * k2, U(2.2 - 1.3 * age)); } if (s.fit.ok && s.fit.amp > 2e-4 ) { const pts2: [number, number][] = []; for (let i = 0 ; i <= n; i++) { const ph = i / n; pts2.push([Xp(ph), mid2 - (s.fit.a * Math.cos(TAU * ph) + s.fit.b * Math.sin(TAU * ph)) * k2]); } api.poly(pts2, { w: U(1.8 ) }); } const at = s.fit.ok && s.fit.amp > 2e-4 ? ` ±${fmt(s.fit.amp, decs(s.fit.amp))} mm` : "ほぼ平ら" ; api.label(narrow ? ` 水面 ${fmt(P.f / 2, 1)} Hz` : ` 水面(印の所) f/2 ${fmt(P.f / 2, 1)} Hz ${at}` , c2.x, c2.y + U(2 ), { size: U(10.5 ) }); const Tw = 2000 / P.f; api.line(r.x, r.y + r.h, r.x + r.w, r.y + r.h, { w: U(1 ) }); for (const ph of [0 , 0.5 , 1 ]) { api.line(Xp(ph), r.y + r.h, Xp(ph), r.y + r.h + U(4 ), { w: U(1 ) }); api.text(` ${fmt(ph * Tw, Tw < 20 ? 1 : 0)}${ph === 1 ? " ms" : ""}` , Xp(ph), r.y + r.h + U(15 ), { size: U(10 ), align: ph === 0 ? "left" : ph === 1 ? "right" : "center" , font: "mono" }); } const phNow = (((s.t / (2 / P.f)) % 1 ) + 1 ) % 1 ; api.line(Xp(phNow), r.y, Xp(phNow), r.y + r.h, { w: U(1 ), dash: [U(1.5 ), U(3 )] }); } /** 振幅の目盛り: 今の Γ を斜線で、しきい値 Γc を太い縦線で */ function gauge(api: Api, s: S, r: Box, U: (n: number) => number) { const P = s.P, fin = Number.isFinite(P.Gc), top = fin ? Math.max(P.Gam, P.Gc) : P.Gam; const M = niceUp(Math.max(0.05 , top * 1.5 )), lab = U(118 ), x1 = r.x + r.w - U(44 ); const X = (v: number) => r.x + lab + (clamp(v, 0 , M) / M) * (x1 - r.x - lab); api.text(` 振幅 Γ ${fmt(P.Gam, 2)} g` , r.x, r.y + U(4 ), { size: U(11.5 ), weight: 600 }); api.line(X(0 ), r.y, X(M), r.y, { w: U(1 ) }); if (P.Gam > 0 ) api.hatch([[X(0 ), r.y - U(5 )], [X(P.Gam), r.y - U(5 )], [X(P.Gam), r.y + U(5 )], [X(0 ), r.y + U(5 )]], { gap: U(3 ), angle: 60 , w: U(1 ) }); api.rect(X(0 ), r.y - U(5 ), X(P.Gam) - X(0 ), U(10 ), { w: U(1 ) }); api.text(` ${fmt(M, M < 1 ? 2 : 1)} g` , X(M) + U(5 ), r.y + U(4 ), { size: U(10 ), font: "mono" }); if (fin && P.Gc <= M) { const x = X(P.Gc); api.line(x, r.y - U(10 ), x, r.y + U(10 ), { w: U(2.2 ) }); api.text(` しきい値 Γc ${fmt(P.Gc, 2)} g` , x, r.y + U(22 ), { size: U(10.5 ), align: x > X(M * 0.6 ) ? "right" : "left" }); } else api.text("しきい値 Γc はこの目盛りより上" , X(M), r.y + U(22 ), { size: U(10.5 ), align: "right" }); } 振動波形の下の段の点は、読み出した高さを「波の1周期のうちのどの時点か」で並べ直したものだ。点に正弦波を最小二乗法で当てはめ、太い線で描いている。
ブランコの立ちこぎ
本文の途中の小さなシミュレーションは、長さの変わる振り子の式を4次のルンゲ・クッタ法で1/240秒ずつ解いている。重心までの長さL L ′ − ( 2 L ′ / L ) θ ′
ブランコの式(src/sims/faraday-waves-swing.ts の 8〜35 行) const g = 9.81 ; // m/s² const C = 0.02 ; // 空気の抵抗(1/秒) const HIST = 900 ; // グラフに残す点の数(30 秒) const SEAT = 0.95 ; // 重心(ふだんの位置)から座板までの長さ m interface S { th: number; om: number; t: number; tp: number; hold: boolean; hist: Float32Array; histL: Float32Array; hn: number; hi: number; hT: number; maxNow: number; peak: number; noted45: boolean; step: (dt: number, f: (h: number) => void ) => number } const L0 = (api: Api) => api.num("len" );const pumpF = (api: Api) => (api.str("pump" ) === "none" ? 0 : (Math.sqrt(g / L0(api)) / TAU) * Number(api.str("pump" )));/** 重心までの長さ L と、その変わる速さ L'(こぐ周期の正弦波) */ function lenAt(api: Api, tp: number): [number, number] { const dL = api.str("pump" ) === "none" ? 0 : api.num("dl" ) / 100 , fp = pumpF(api), w = TAU * fp; return [L0(api) - dL * Math.cos(w * tp), dL * w * Math.sin(w * tp)]; } function deriv(api: Api, th: number, om: number, tp: number): [number, number] { const [L, Ld] = lenAt(api, tp); return [om, (-2 * Ld * om) / L - (g / L) * Math.sin(th) - C * om]; } /** 4 次のルンゲ・クッタ法で h 秒進める */ function rk4(api: Api, s: S, h: number) { const { th, om, tp } = s; const [a1, b1] = deriv(api, th, om, tp); const [a2, b2] = deriv(api, th + (h / 2 ) * a1, om + (h / 2 ) * b1, tp + h / 2 ); const [a3, b3] = deriv(api, th + (h / 2 ) * a2, om + (h / 2 ) * b2, tp + h / 2 ); const [a4, b4] = deriv(api, th + h * a3, om + h * b3, tp + h); s.th = th + (h / 6 ) * (a1 + 2 * a2 + 2 * a3 + a4); s.om = om + (h / 6 ) * (b1 + 2 * b2 + 2 * b3 + b4); s.tp = tp + h; } 
