このシミュレーションは、海を一辺96〜192個のますに分け、1コマごとに「水面の傾きで流れを速める → 流れの出入りで水面を上げ下げする」をくり返している。波の速さg h
長さはメートル、時間は秒
計算はすべてメートルと秒で行う。画面の中の海は一辺100の単位で描き、場面ごとに縮尺（W ÷ L）を掛けて直す。640 mの浜も320 kmの海も同じ大きさに見えるのは、このためだ。
定数と単位(src/sims/shallow-waves.ts の 26〜39 行) const G = 9.81 ; // 重力加速度 m/s² const NQ = [96 , 128 , 192 ]; // 一辺のますの数(画質 0・1・2) const STEPMAX = [12 , 8 , 6 ]; // 1 コマに進める刻みの上限(これより速い再生は、計算が追いつく分だけになる) const SPIN = [10e6 , 18e6 , 32e6 ]; // 場面の最初に、波を浜まで届かせておく計算の量の上限(ますの数 × 刻みの数) const SPINMS = 6 ; // その早送りに 1 コマで使ってよい時間(ミリ秒・1/60 秒のコマあたり) const CFL = 0.5 ; // c Δt / Δx(この方法は 2 次元で 0.71 まで安定) const KB = 0.39 ; // 砕ける高さ ÷ 水深(静かな水面から山まで) const W = 100 ; // 画面の中の海の一辺(描画の単位)。縮尺は場面ごとに W ÷ L const SPONGE = 0.06 ; // 端で波を吸い込む帯の幅(一辺に対する割合) const ZONE = 0.14 ; // 造波板の前で波を作る帯の幅(一辺に対する割合) const ZP = 0.5 - SPONGE - 0.01 ; // 造波板の位置(海の中心から沖へ・一辺に対する割合) const UNI = ZP - ZONE - 0.02 ; // ここより沖は水深がそろっている(造波板の作る波が乱れないように) const SPEEDS = ["1" , "3" , "10" , "30" , "100" , "300" ];const NSEG = 24 ; // 造波板の板の数 CFLは、1回の刻みで波が進む距離を、ます1つの幅の何割までにするかを決める数である。この計算方法は2次元で0.71を超えると、刻みごとに誤差がふくらんで計算が発散する。余裕をみて0.5にした。KBは砕ける高さの目安で、静かな水面から山までが水深の0.39倍、つまり波高が水深の約0.78倍になった所で波を崩す。
海底の形を関数で決める
場面ごとに、位置( x , z )
海底の形(src/sims/shallow-waves.ts の 41〜108 行) interface Spec { L: number; // 海の一辺 m h0: number; // 沖の水深 m hcap: number; // 掘れる深さの上限 m(刻み Δt はこの深さの波の速さで決める) hmin: number; // これより浅い所は陸(動かない壁) top: number; // 盛れる高さの上限 m rate: number; // 指でなぞったとき、海底が 1 秒に上下する高さ m fric: number; // 底の摩擦で流れが弱まる割合(1 秒あたり。津波のような長い波では、ほとんど効かない大きさにする) per: boolean; // 左右の端をつなぐ maker: boolean; // 造波板があるか vb: number; vw: number; // 縦の倍率(海底・波) hide: boolean; // 水の下の海底を描かない(波を海底より大きく引き伸ばす場面で、谷が海底を突き抜けて見えないように) levels: number[]; // 等深線を引く深さ m bar: number; barText: string; // 物差しの長さ m と文字 buoys: [number, number][]; // ブイ A・B の位置(一辺に対する割合。z は沖が正) cam: [number, number, number]; // 視点の向き・見下ろす角度(rad)・距離 bed(x: number, z: number, L: number): number; // 海底の高さ m(静かな海面が 0・上が正)。x は右、z は沖が正 } /** 砂浜の断面: 岸から d(m)離れた所の水深。D で沖の深さ h0 に届く(砂浜でよく見られる、d の 2/3 乗に比例する形) */ const dean = (d: number, h0: number, D: number) => h0 * Math.min(1 , Math.pow(Math.max(0 , d) / D, 2 / 3 ));/** 岬の海岸線の z(m)。真ん中が沖へ張り出し、左右の端は湾になる(端どうしがつながる形) */ const capeZ = (x: number, L: number) => -0.36 * L + 0.3 * L * (0.5 + 0.5 * Math.cos((TAU * x) / L)) ** 3 ;const even = (a: number, b: number, n: number) => Array.from ({ length: n }, (_, i) => a + i * b); const SPECS: Record<string, Spec> = { beach: { L: 640 , h0: 14 , hcap: 18 , hmin: 0.3 , top: 6 , rate: 5 , fric: 0.002 , per: true , maker: true , vb: 6 , vw: 6 , hide: false , levels: even(2 , 2 , 6 ), bar: 100 , barText: "100 m" , buoys: [[0.12 , 0.16 ], [0.12 , -0.2 ]], cam: [0.55 , 0.8 , 152 ], bed: (x, z, L) => { const d = z + 0.26 * L; return d > 0 ? -dean(d, 14 , (UNI + 0.26 ) * L) : Math.min(5 , -0.05 * d); }, }, cape: { L: 640 , h0: 14 , hcap: 18 , hmin: 0.3 , top: 8 , rate: 5 , fric: 0.002 , per: true , maker: true , vb: 6 , vw: 6 , hide: false , levels: even(2 , 2 , 6 ), bar: 100 , barText: "100 m" , buoys: [[0 , -0.02 ], [-0.46 , -0.27 ]], cam: [0.15 , 0.95 , 150 ], bed: (x, z, L) => { const zs = capeZ(x, L), d = z - zs; return d > 0 ? -dean(d, 14 , UNI * L - zs) : Math.min(10 , -0.08 * d); }, }, island: { L: 720 , h0: 12 , hcap: 16 , hmin: 0.3 , top: 12 , rate: 5 , fric: 0.002 , per: false , maker: true , vb: 6 , vw: 6 , hide: false , levels: even(2 , 2 , 5 ), bar: 100 , barText: "100 m" , buoys: [[0.12 , 0.16 ], [0 , -0.16 ]], cam: [0.35 , 1.0 , 145 ], bed: (x, z, L) => { const r = Math.hypot(x, z + 0.04 * L), Ri = 0.06 * L, Rs = 0.27 * L; return r < Ri ? 12 * (1 - (r / Ri) ** 2 ) : -12 * Math.min(1 , Math.pow((r - Ri) / (Rs - Ri), 0.8 )); }, }, harbor: { L: 640 , h0: 10 , hcap: 14 , hmin: 0.3 , top: 6 , rate: 5 , fric: 0.002 , per: false , maker: true , vb: 6 , vw: 6 , hide: false , levels: even(2 , 2 , 4 ), bar: 100 , barText: "100 m" , buoys: [[0 , -0.12 ], [-0.27 , -0.12 ]], cam: [0.2 , 1.25 , 140 ], bed: (x, z, L) => { if (Math.abs(z + 0.02 * L) < 0.012 * L && Math.abs(x) > 0.055 * L) return 3 ; // 防波堤(すき間の幅 70 m) const d = z + 0.41 * L; return d > 0 ? -dean(d, 10 , 0.17 * L) : Math.min(4 , -0.04 * d); }, }, tsunami: { L: 320000 , h0: 4000 , hcap: 5200 , hmin: 40 , top: 400 , rate: 1500 , fric: 2e-5 , per: false , maker: false , vb: 10 , vw: 12000 , hide: true , levels: [100 , 150 , 200 , 500 , 1000 , 2000 , 3000 ], bar: 50000 , barText: "50 km" , buoys: [[0 , 0.16 ], [0 , -0.22 ]], cam: [-0.75 , 0.62 , 150 ], bed: (x, z, L) => { const zs = -0.3 * L, zb = -0.03 * L, zd = 0.12 * L; if (z < zs) return Math.min(300 , 0.004 * (zs - z)); // 陸 if (z < zb) return -(60 + (140 * (z - zs)) / (zb - zs)); // 大陸棚(60〜200 m) if (z < zd) { const t = (z - zb) / (zd - zb); return -(200 + 3800 * t * t * (3 - 2 * t)); } // 大陸斜面 return -4000 ; // 深海 }, }, cliff: { L: 640 , h0: 15 , hcap: 19 , hmin: 0.3 , top: 12 , rate: 5 , fric: 0.002 , per: true , maker: true , vb: 6 , vw: 6 , hide: false , levels: even(2 , 2 , 7 ), bar: 100 , barText: "100 m" , buoys: [[0.2 , -0.212 ], [0.02 , -0.17 ]], cam: [0.75 , 0.66 , 148 ], bed: (x, z, L) => (z < -0.22 * L ? 12 : -15 ), }, }; 浜の断面は、岸からの距離の2/3乗で深くなる形（dean）にした。砂浜でよく見られる形で、岸の近くほど急に深くなる。岬は、海岸線そのものを沖へ張り出させ、深さはそこからの距離で決めている。こうすると、等深線も岬に沿って沖へふくらみ、波を岬へ集める形になる。hcapは掘れる深さの上限で、刻みの長さはこの深さの波の速さで決める。指で深く掘っても計算が壊れないようにするためだ。
流れを進め、水面を進める
本文の2本の式を、そのまま格子の上で解いている。水面の高さη P = h u Q
浅水の方程式(src/sims/shallow-waves.ts の 137〜184 行) /** 1 刻み Δt 進める。流量 P・Q を水面の傾きで進め(運動の式)、その流量の出入りで水面 η を進める(水の量の式) */ function advance(s: S, dt: number) { const { N, eta, P, Q, cP, cQ, sP, sQ, sE, h, wet, lim } = s; const N1 = N + 1 , fr = 1 - s.sp.fric * dt; // x の向きの流量: ∂P/∂t = −g h ∂η/∂x(cP = g h Δt ÷ Δx。陸との境目は 0 にしてあり、水が通らない) for (let j = 0 ; j < N; j++) { const r = j * N, f0 = j * N1; for (let i = 1 ; i < N; i++) { const f = f0 + i; P[f] = (P[f] * fr - cP[f] * (eta[r + i] - eta[r + i - 1 ])) * sP[f]; } if (s.sp.per) { P[f0] = (P[f0] * fr - cP[f0] * (eta[r] - eta[r + N - 1 ])) * sP[f0]; P[f0 + N] = P[f0]; } else { // いちばん外の境目は、外へ出ていく波だけを通す(進む波では P = c η。左の端は外が −x なので符号が逆) P[f0] = wet[r] ? -Math.sqrt(G * h[r]) * eta[r] : 0 ; P[f0 + N] = wet[r + N - 1 ] ? Math.sqrt(G * h[r + N - 1 ]) * eta[r + N - 1 ] : 0 ; } } // z の向きの流量(境目 f は、ます f − N と f のあいだ) const nn = N * N; for (let f = N; f < nn; f++) Q[f] = (Q[f] * fr - cQ[f] * (eta[f] - eta[f - N])) * sQ[f]; const last = nn - N; for (let i = 0 ; i < N; i++) { Q[i] = wet[i] ? -Math.sqrt(G * h[i]) * eta[i] : 0 ; Q[nn + i] = wet[last + i] ? Math.sqrt(G * h[last + i]) * eta[last + i] : 0 ; } // 水面: ∂η/∂t = −(∂P/∂x + ∂Q/∂z)。砕ける高さを超えた分は削る const r = dt / s.dx, clip = s.clip, cr = s.cr; let nc = 0 ; for (let j = 0 ; j < N; j++) { const k0 = j * N, f0 = j * N1; for (let i = 0 ; i < N; i++) { const k = k0 + i; if (!wet[k]) continue ; const f = f0 + i; let e = (eta[k] - r * (P[f + 1 ] - P[f] + Q[k + N] - Q[k])) * sE[k]; const m = lim[k]; if (e > m || e < -m) { clip[nc] = k; cr[nc++] = m / Math.abs(e); e = e > 0 ? m : -m; } eta[k] = e; } } // 砕けた所は、流れも同じ割合で弱める(波のエネルギーがそこで失われる)。泡の時刻を記す for (let n = 0 ; n < nc; n++) { const k = clip[n], q = cr[n], f = k + ((k / N) | 0 ); P[f] *= q; P[f + 1 ] *= q; Q[k] *= q; Q[k + N] *= q; s.tb[k] = s.time; } if (s.sp.per) for (let j = 0 ; j < N; j++) P[j * N1 + N] = P[j * N1]; s.brk = nc; s.time += dt; wave(s, dt); } 前半のP[f] = (P[f] * fr - cP[f] * (eta[...] - eta[...]))は、本文の∂ u / ∂ t = − g ∂ η / ∂ x h cPにはg h Δ t /Δ x ∂ η / ∂ t = − ( ∂ P / ∂ x + ∂ Q / ∂ z )
いちばん外の境目では、外へ出ていく波だけを通す。進む波では流量がcη
造波板と、波を吸い込む帯
造波板の前の帯では、水面と流量を、造波板が作るはずの進む波（η = a sin ( ω t − k ⋅ x ) = cη
造波と吸収(src/sims/shallow-waves.ts の 186〜270 行) /** 造波板の波: 周期 T・沖の波長・向き。左右の端をつなぐ場面では、向きを端でつながる値に合わせる */ function waveParams(api: Api, s: S, sync: boolean) { const w = s.wave, sp = s.sp; w.T = clamp(api.num("T" ) || 10 , 2 , 40 ); w.om = TAU / w.T; w.c0 = Math.sqrt(G * sp.h0); w.k0 = w.om / w.c0; let th = (clamp(api.num("theta" ) || 0 , -45 , 45 ) * Math.PI) / 180 ; if (sp.per) { // 一辺 L に入る波の数が、横の向きに整数 n になる向きだけを使う(sin θ = n λ ÷ L) const m = (s.L * w.k0) / TAU, nmax = Math.floor(m * Math.SQRT1_2 + 1e-9 ); const n = clamp(Math.round(m * Math.sin(th)), -nmax, nmax); w.kx = (TAU * n) / s.L; th = Math.asin(clamp(w.kx / w.k0, -1 , 1 )); const deg = Math.round((th * 180 ) / Math.PI); if (sync && deg !== Math.round(api.num("theta" ))) api.set("theta" , deg); } else w.kx = w.k0 * Math.sin(th); w.th = th; w.kz = -Math.sqrt(Math.max(0 , w.k0 * w.k0 - w.kx * w.kx)); w.aim = sp.maker ? Math.max(0 , api.num("H" )) / 2 : 0 ; zone(s); } /** 造波板の前の帯のますと、そこで目指す波の位相(cos と sin)を用意する */ function zone(s: S) { const { N, L, dx } = s, w = s.wave, Z = s.zone, zp = ZP * L, z0 = (ZP - ZONE) * L, ws = Math.max(6 , Math.round(SPONGE * N)); let n = 0 ; if (s.sp.maker) for (let j = 0 ; j < N; j++) { const z = -L / 2 + (j + 0.5 ) * dx; if (z < z0 || z > zp) continue ; const t = (z - z0) / (zp - z0), al = 1 - Math.exp(-3 * w.om * t * t * s.dt); // 造波板に近いほど強く近づける for (let i = 0 ; i < N; i++) { if (!s.sp.per && (i < ws || i > N - 1 - ws)) continue ; const k = j * N + i; if (!s.wet[k] || n >= Z.k.length) continue ; const x = cx(s, i), pe = w.kx * x + w.kz * (z - zp), pp = pe - (w.kx * dx) / 2 , pq = pe - (w.kz * dx) / 2 ; Z.k[n] = k; Z.al[n] = al; Z.cE[n] = Math.cos(pe); Z.sE[n] = Math.sin(pe); Z.cP[n] = Math.cos(pp); Z.sP[n] = Math.sin(pp); Z.cQ[n] = Math.cos(pq); Z.sQ[n] = Math.sin(pq); n++; } } Z.n = n; } /** 造波板: 帯の中の水面と流量を、進む波 η = a sin(ωt − k·x)・流量 = c η(波の向き)に近づける。戻ってきた波もここで吸い込まれる */ function wave(s: S, dt: number) { const w = s.wave, Z = s.zone; w.ph = (w.ph + w.om * dt) % TAU; w.a += (w.aim - w.a) * Math.min(1 , dt / w.T); // 振れ幅は 1 周期ほどかけて変える(急に動かすと乱れる) if (!Z.n) return ; const { eta, P, Q, cP, cQ, N } = s, a = w.a, sn = Math.sin(w.ph), cs = Math.cos(w.ph), ux = (w.c0 * w.kx) / w.k0, uz = (w.c0 * w.kz) / w.k0; for (let n = 0 ; n < Z.n; n++) { const k = Z.k[n], al = Z.al[n], f = k + ((k / N) | 0 ); eta[k] += al * (a * (sn * Z.cE[n] - cs * Z.sE[n]) - eta[k]); // sin(φ − ψ) = sin φ cos ψ − cos φ sin ψ if (cP[f] > 0 ) P[f] += al * (ux * a * (sn * Z.cP[n] - cs * Z.sP[n]) - P[f]); if (cQ[k] > 0 ) Q[k] += al * (uz * a * (sn * Z.cQ[n] - cs * Z.sQ[n]) - Q[k]); } } /** 海底の形が決まったら(または変わったら)、水深・陸・境目の係数・端の吸い込み・表示の形を作り直す */ function depth(s: S) { const { N, dx, dt, sp, bed, h, wet, eta, P, Q, cP, cQ, sE, sP, sQ, lim } = s, N1 = N + 1 , nn = N * N; for (let k = 0 ; k < nn; k++) { const d = -bed[k]; wet[k] = d > sp.hmin ? 1 : 0 ; h[k] = wet[k] ? d : 0 ; lim[k] = KB * h[k]; if (!wet[k]) eta[k] = 0 ; } const gdt = (G * dt) / dx, cf = (a: number, b: number) => (wet[a] && wet[b] ? gdt * 0.5 * (h[a] + h[b]) : 0 ); // 端の帯: 外へ近いほど強く弱める(帯を通り抜ける間に波がおよそ e の 3 乗分の 1 になる強さ) const ws = Math.max(6 , Math.round(SPONGE * N)), c0 = Math.sqrt(G * sp.h0); for (let j = 0 ; j < N; j++) for (let i = 0 ; i < N; i++) { const k = j * N + i, d = Math.min(j, N - 1 - j, sp.per ? 1e9 : i, sp.per ? 1e9 : N - 1 - i); if (d < ws && wet[k]) { const c = Math.max(Math.sqrt(G * h[k]), 0.3 * c0), sig = ((9 * c) / (ws * dx)) * ((ws - d) / ws) ** 2 ; sE[k] = Math.exp(-sig * dt); } else sE[k] = 1 ; } for (let j = 0 ; j < N; j++) { const r = j * N, f0 = j * N1; for (let i = 1 ; i < N; i++) { cP[f0 + i] = cf(r + i - 1 , r + i); sP[f0 + i] = 0.5 * (sE[r + i - 1 ] + sE[r + i]); } if (sp.per) { cP[f0] = cP[f0 + N] = cf(r + N - 1 , r); sP[f0] = sP[f0 + N] = 0.5 * (sE[r + N - 1 ] + sE[r]); } else { cP[f0] = cP[f0 + N] = 0 ; sP[f0] = sP[f0 + N] = 1 ; } for (let i = 0 ; i <= N; i++) if (cP[f0 + i] === 0 && (sp.per || (i > 0 && i < N))) P[f0 + i] = 0 ; } for (let f = N; f < nn; f++) { cQ[f] = cf(f - N, f); sQ[f] = 0.5 * (sE[f - N] + sE[f]); if (cQ[f] === 0 ) Q[f] = 0 ; } for (let i = 0 ; i < N; i++) { cQ[i] = cQ[nn + i] = 0 ; sQ[i] = sQ[nn + i] = 1 ; } shapeGfx(s); } 近づける割合は造波板に近いほど強い。帯の中では、浜や崖で跳ね返って戻ってきた波も目標の波へ引き戻されるので、帯は波を作る所と吸い込む所を兼ねる。波の形のsin ( ω t − ψ ) ψ cos sin L sin θ = nλ / L
海の端の帯では、波を外側ほど強く弱める。帯を抜ける間に波がおよそe − 3
ブイで測る
ブイの所の水面を、刻みごとに記録している。
測る(src/sims/shallow-waves.ts の 272〜297 行) /** ブイの所の水面の上下(波高)・波の向き・いちばん高い水位を測る。向きは「水面の高さ × 流量」の平均(波のエネルギーが運ばれる向き) */ function measure(s: S, dt: number) { const win = s.sp.maker ? 1.5 * s.wave.T : 900 ; for (const b of s.buoys) { const k = b.k; if (k < 0 || !s.wet[k]) continue ; const e = s.eta[k], f = k + ((k / s.N) | 0 ); const pc = 0.5 * (s.P[f] + s.P[f + 1 ]), qc = 0.5 * (s.Q[k] + s.Q[k + s.N]); if (e > b.mx) b.mx = e; if (e < b.mn) b.mn = e; if (e > b.top) b.top = e; b.tw += dt; if (b.tw >= win) { b.pmx = b.mx; b.pmn = b.mn; b.mx = -1e9 ; b.mn = 1e9 ; b.tw = 0 ; } const a = Math.min(1 , dt / win); b.ex += (e * pc - b.ex) * a; b.ez += (e * qc - b.ez) * a; b.e2 += (e * e - b.e2) * a; } } const resetBuoy = (b: Buoy) => { b.mx = b.pmx = -1e9 ; b.mn = b.pmn = 1e9 ; b.tw = 0 ; b.ex = b.ez = b.e2 = 0 ; b.top = 0 ; };/** 波高(山から谷まで)。直近の 2 つの区間(1.5 周期ずつ)の、いちばん高い所と低い所の差 */ const heightOf = (b: Buoy) => { const hi = Math.max(b.mx, b.pmx), lo = Math.min(b.mn, b.pmn); return hi > lo ? hi - lo : 0 ; };/** 波の向き(度・まっすぐ岸へ向かう向きから右回りが正)。進む波でなければ(跳ね返った波と重なっている所など)null */ function dirOf(s: S, b: Buoy) { const c = Math.sqrt(G * s.h[b.k]), m = Math.hypot(b.ex, b.ez); if (!(b.e2 > 1e-8 ) || m < 0.35 * c * b.e2) return null ; return (Math.atan2(b.ex, -b.ez) * 180 ) / Math.PI; } 波高は、1.5周期ずつの2つの区間で、いちばん高い所と低い所の差をとる。波の向きは、水面の高さと流量を掛けた値の平均から求める。進む波では、山の所で水が前へ、谷の所で後ろへ流れるので、この平均は波が進む向き（エネルギーが運ばれる向き）を指す。崖の前のように、来る波と跳ね返った波が重なってその場で上下している所では平均がほとんど0になり、「定まらない」と出る。
指の操作
ブイと矢印の取っ手は、画面の上での距離でつかむ。指の下の海面の位置は、視点から指の方向へのばした線と、高さ0の平面との交点で求める。
操作(src/sims/shallow-waves.ts の 355〜431 行) const SEA = new THREE.Plane(new THREE.Vector3(0 , 1 , 0 ), 0 );/** 画面の点の下にある海面の位置(m)。海の外なら null */ function ground(s: S, x: number, y: number) { const p = s.v.onPlane(x, y, SEA); if (!p) return null ; const X = p.x / s.k, Z = p.z / s.k; return Math.abs(X) < s.L / 2 && Math.abs(Z) < s.L / 2 ? { x: X, z: Z } : null ; } const ARX = -0.28 ; // 波の向きの矢印の根元(一辺に対する割合) const arrowAt = (s: S) => { const x0 = ARX * s.L, zc = (ZP - ZONE - 0.05 ) * s.L, len = 0.16 * s.L, w = s.wave; return [x0, zc, x0 + (len * w.kx) / w.k0, zc + (len * w.kz) / w.k0]; };/** 指の下につかめる物(ブイ A・B・波の向きの取っ手)があるか */ function grab(api: Api, s: S) { const p = api.pointer, R = api.u(28 ); let best: string | null = null , bd = R; s.buoys.forEach((b, n) => { const q = s.v.project(b.m.position); const d = Math.hypot(q.x - p.x, q.y - p.y); if (q.on && d < bd) { bd = d; best = "AB" [n]; } }); if (s.sp.maker) { const a = arrowAt(s), q = s.v.project([a[2 ] * s.k, 0.6 , a[3 ] * s.k]), d = Math.hypot(q.x - p.x, q.y - p.y); if (q.on && d < bd) best = "arrow" ; } return best; } function dragTo(api: Api, s: S, what: string) { const g = ground(s, api.pointer.x, api.pointer.y); if (!g) return ; if (what === "arrow" ) { const zc = (ZP - ZONE - 0.05 ) * s.L, dx = g.x - ARX * s.L, dz = g.z - zc; if (Math.hypot(dx, dz) < 0.04 * s.L) return ; const deg = clamp((Math.atan2(dx, -dz) * 180 ) / Math.PI, -45 , 45 ); api.set("theta" , Math.round(deg)); waveParams(api, s, true ); return ; } const b = s.buoys[what === "A" ? 0 : 1 ], k = cellAt(s, g.x, g.z); if (k < 0 || !s.wet[k]) return ; b.fx = g.x; b.fz = g.z; b.k = k; resetBuoy(b); } /** 指でなぞった所の海底を盛る・掘る(なめらかな山の形に) */ function sculpt(api: Api, s: S) { const c = s.sculpt; if (!c) return ; const sg = api.str("touch" ) === "raise" ? 1 : -1 , R = 0.05 * s.L, amt = sg * s.sp.rate * api.dt, { N } = s; for (let j = 0 ; j < N; j++) { const dz = (cx(s, j) - c.z) / R; if (Math.abs(dz) > 3 ) continue ; for (let i = 0 ; i < N; i++) { const dx = (cx(s, i) - c.x) / R; if (Math.abs(dx) > 3 ) continue ; const k = j * N + i; s.bed[k] = clamp(s.bed[k] + amt * Math.exp(-0.5 * (dx * dx + dz * dz)), -s.sp.hcap, s.sp.top); } } depth(s); zone(s); if (!s.noted.sculpt) { s.noted.sculpt = true ; api.note(sg > 0 ? "盛った所は浅くなり、その上で波が遅くなります" : "掘った所は深くなり、その上で波が速くなります" , 2600 ); } } /** たたいた所に岩を落とす。津波の場面では、その下の海底が持ち上がる(地震) */ function tapAt(api: Api, s: S) { const g = ground(s, api.pointer.x, api.pointer.y); if (!g) return ; const k = cellAt(s, g.x, g.z); if (k < 0 || !s.wet[k]) { if (!s.noted.land) { s.noted.land = true ; api.note("陸の上です。海の上をたたくと岩が落ちます" , 2000 ); } return ; } if (s.key === "tsunami" ) { uplift(s, g.x, g.z, Math.max(0.5 , api.num("H" )), 0.05 * s.L, 0.05 * s.L); api.note("海底が持ち上がり、真上の海面も同じ形に盛り上がりました" , 2400 ); return ; } if (s.rocks.length >= 8 ) { const r = s.rocks.shift(); if (r) s.v.scene.remove(r.m); } const m = new THREE.Mesh(s.rockGeo, s.mats.rock); m.castShadow = true ; s.v.scene.add(m); s.rocks.push({ m, x: g.x, z: g.z, y: 30 , vy: 0 , hit: false }); } /** 岩は本物の重力で落ち、海面に着いた所の水を押し下げる */ function rocks(s: S, dt: number) { for (const r of s.rocks) { r.vy -= G * dt; r.y += r.vy * dt; if (!r.hit && r.y <= 0 ) { r.hit = true ; const w = Math.max(10 , 2 * s.dx); uplift(s, r.x, r.z, -2.5 , w, w); } // 岩が押しのけた水(格子で表せる幅にしてある) r.m.position.set(r.x * s.k, r.y * s.sp.vb * s.k + 1.1 , r.z * s.k); } const gone = s.rocks.filter((r) => r.y < -s.sp.hcap); for (const r of gone) s.v.scene.remove(r.m); if (gone.length) s.rocks = s.rocks.filter((r) => r.y >= -s.sp.hcap); } 海底を盛る・掘るときは、指の下を中心に、なめらかな山の形で海底を上下させる。そのたびに水深と境目の係数、等深線を作り直す。岩は本物の重力加速度で落ち、海面に着いた所の水面を押し下げる。この凹みが輪になって広がる。
three.jsが受け持つ描画
水面と海底は、どちらも格子の点を頂点にした面（heightGrid）で描いている。
描画(src/sims/shallow-waves.ts の 433〜665 行) // 水面は自前の材質(明るさだけを出す)。光の当たり方・照り返し・高さ・泡(砕けた所)・等深線を 1 つの明るさにまとめる。白黒の網点への変換は共通部品がする const WATER_VS = /* glsl */ ` attribute float aC; attribute float aF; attribute float aE; varying vec3 vN; varying vec3 vP; varying float vC; varying float vF; varying float vE; void main() { vec4 w = modelMatrix * vec4(position, 1.0); vP = w.xyz; vN = normalize(mat3(modelMatrix) * normal); vC = aC; vF = aF; vE = aE; gl_Position = projectionMatrix * viewMatrix * w; }` ;const WATER_FS = /* glsl */ ` uniform vec3 uSun; uniform float uLine; varying vec3 vN; varying vec3 vP; varying float vC; varying float vF; varying float vE; void main() { vec3 n = normalize(vN), l = normalize(uSun), e = normalize(cameraPosition - vP); // 平らな水面の明るさを 0.28 にして、傾きによる光の当たり方の差を強めに見せる(小さな波の形も網点で読めるように。大きな波は白と黒に振り切れる) float dif = max(dot(n, l), 0.0) - l.y; float spe = pow(max(dot(n, normalize(l + e)), 0.0), 48.0); float b = 0.28 + 2.2 * dif + 0.45 * spe + 0.2 * clamp(vE, -1.5, 1.5); b = max(b, 0.95 * vF); // 等深線: 深さの目盛り vC が整数になる所に細い線(線が詰まりすぎる急な所では消す) // (深さがそろって目盛りが変わらない所には引かない) float fw = fwidth(vC), d = abs(fract(vC + 0.5) - 0.5), w = max(fw, 1e-4); float ln = (1.0 - smoothstep(0.5 * w, 1.5 * w, d)) * (1.0 - smoothstep(0.2, 0.45, fw)) * step(1e-4, fw) * step(0.5, vC) * uLine; b = mix(b, 0.8, ln); gl_FragColor = vec4(vec3(clamp(b, 0.0, 1.0)), 1.0); }` ;const SUN: [number, number, number] = [-0.35 , 0.55 , 0.76 ]; // 太陽へ向かう向き(沖の左上。カメラの後ろから照らし、平らな水面が光らないように) /** 一辺 N の格子の配列と、水面・海底・断面の形を用意する */ function alloc(s: S, N: number) { const nn = N * N, nf = (N + 1 ) * N, zn = N * (Math.ceil(ZONE * N) + 3 ); s.N = N; s.bed = new Float32Array(nn); s.h = new Float32Array(nn); s.wet = new Uint8Array(nn); s.eta = new Float32Array(nn); s.P = new Float32Array(nf); s.Q = new Float32Array(nf); s.cP = new Float32Array(nf); s.cQ = new Float32Array(nf); s.sP = new Float32Array(nf); s.sQ = new Float32Array(nf); s.sE = new Float32Array(nn); s.lim = new Float32Array(nn); s.tb = new Float32Array(nn); s.clip = new Int32Array(nn); s.cr = new Float32Array(nn); s.zone = { k: new Int32Array(zn), al: new Float32Array(zn), cE: new Float32Array(zn), sE: new Float32Array(zn), cP: new Float32Array(zn), sP: new Float32Array(zn), cQ: new Float32Array(zn), sQ: new Float32Array(zn), n: 0 }; const v = s.v; if (s.g) for (const m of [s.g.water, s.g.land.mesh, s.g.soil, s.g.wall]) { v.scene.remove(m); m.geometry.dispose(); } const side = (W * (N - 1 )) / N; const wg = heightGrid(N - 1 , N - 1 , side, side, s.mats.water); const aC = new Float32Array(nn), aF = new Float32Array(nn), aE = new Float32Array(nn); wg.geo.setAttribute("aC" , new THREE.BufferAttribute(aC, 1 )); wg.geo.setAttribute("aF" , new THREE.BufferAttribute(aF, 1 )); wg.geo.setAttribute("aE" , new THREE.BufferAttribute(aE, 1 )); wg.mesh.frustumCulled = false ; wg.mesh.receiveShadow = false ; const land = heightGrid(N - 1 , N - 1 , side, side, s.mats.land); land.mesh.castShadow = true ; land.mesh.frustumCulled = false ; const soil = new THREE.Mesh(skirtGeo(N), s.mats.soil), wall = new THREE.Mesh(skirtGeo(N), s.mats.wall); soil.frustumCulled = wall.frustumCulled = false ; v.scene.add(wg.mesh, land.mesh, soil, wall); s.g = { water: wg.mesh, wpos: (wg.geo.attributes.position as THREE.BufferAttribute).array as Float32Array, wnor: (wg.geo.attributes.normal as THREE.BufferAttribute).array as Float32Array, wy: new Float32Array(nn), aC, aF, aE, land, landY: new Float32Array(nn), soil, wall, base: 0 , }; } /** 海の四辺の断面(縦の帯)。辺ごとに N 個の点の上と下 */ function skirtGeo(N: number) { const geo = new THREE.BufferGeometry(), idx: number[] = []; geo.setAttribute("position" , new THREE.BufferAttribute(new Float32Array(4 * N * 2 * 3 ), 3 )); for (let e = 0 ; e < 4 ; e++) for (let n = 0 ; n < N - 1 ; n++) { const a = (e * N + n) * 2 ; idx.push(a, a + 1 , a + 2 , a + 1 , a + 3 , a + 2 ); } geo.setIndex(idx); return geo; } /** 辺 e(0 = 岸の側・1 = 沖の側・2 = 左・3 = 右)の n 番目のます */ const edgeCell = (N: number, e: number, n: number) => (e === 0 ? n : e === 1 ? (N - 1 ) * N + n : e === 2 ? n * N : n * N + N - 1 ); /** 海底の形が変わったときの表示: 陸と海底の高さ・等深線の目盛り・断面 */ function shapeGfx(s: S) { const g = s.g; if (!g) return ; const { N, sp, bed, wet, h } = s, nn = N * N, yb = sp.vb * s.k, lv = sp.levels; g.base = -sp.hcap * yb - 2 ; for (let k = 0 ; k < nn; k++) { g.landY[k] = sp.hide && wet[k] ? g.base - 2 : bed[k] * yb; // 等深線の目盛り: 深さが levels[n] のとき n + 1 const d = h[k]; let c = 0 ; if (wet[k]) { if (d <= lv[0 ]) c = d / lv[0 ]; else { let n = 0 ; while (n < lv.length - 1 && d > lv[n + 1 ]) n++; c = n + 1 + (n < lv.length - 1 ? (d - lv[n]) / (lv[n + 1 ] - lv[n]) : 0 ); } } g.aC[k] = c; } // 陸の点は、となりの水の点と同じ目盛りにする(岸の三角形に等深線が引かれないように) for (let k = 0 ; k < nn; k++) if (!wet[k]) { const i = k % N; g.aC[k] = Math.max(i > 0 && wet[k - 1 ] ? g.aC[k - 1 ] : 0 , i < N - 1 && wet[k + 1 ] ? g.aC[k + 1 ] : 0 , k >= N && wet[k - N] ? g.aC[k - N] : 0 , k < nn - N && wet[k + N] ? g.aC[k + N] : 0 ); } g.land.set(g.landY); (g.water.geometry as THREE.BufferGeometry).attributes.aC.needsUpdate = true ; const pos = (g.soil.geometry as THREE.BufferGeometry).attributes.position as THREE.BufferAttribute, a = pos.array as Float32Array; for (let e = 0 ; e < 4 ; e++) for (let n = 0 ; n < N; n++) { const k = edgeCell(N, e, n), i = k % N, j = (k / N) | 0 , o = (e * N + n) * 6 , X = cx(s, i) * s.k, Z = cx(s, j) * s.k; a[o] = X; a[o + 1 ] = bed[k] * yb; a[o + 2 ] = Z; a[o + 3 ] = X; a[o + 4 ] = g.base; a[o + 5 ] = Z; } pos.needsUpdate = true ; g.soil.geometry.computeVertexNormals(); const wp = (g.wall.geometry as THREE.BufferGeometry).attributes.position.array as Float32Array; wp.set(a); } /** 毎コマ: 水面の高さ・面の向き・泡・断面の水の上の辺・造波板・ブイ */ function draw(api: Api, s: S) { const g = s.g; if (!g) return ; const { N, sp, eta, wet, bed } = s, nn = N * N, yw = sp.vw * s.k, yb = sp.vb * s.k, y = g.wy, pos = g.wpos, nor = g.wnor; const aref = sp.maker ? Math.max(0.25 , api.num("H" ) / 2 ) : 0.3 , foamT = sp.maker ? 0.6 * s.wave.T : 300 , t = s.time; for (let k = 0 ; k < nn; k++) { const yy = wet[k] ? eta[k] * yw : g.base - 1 ; // 陸の点では水面を海底の下へ沈め、岸で陸の面と交わらせる y[k] = yy; pos[3 * k + 1 ] = yy; g.aE[k] = wet[k] ? eta[k] / aref : 0 ; const f = 1 - (t - s.tb[k]) / foamT; g.aF[k] = f > 0 ? f : 0 ; } // 面の向き: 水面の傾き(となりのますとの差)から。陸のとなりでは片側の差を使う const d2 = 2 * s.dx * s.k; for (let j = 0 ; j < N; j++) for (let i = 0 ; i < N; i++) { const k = j * N + i, o = 3 * k; if (!wet[k]) { nor[o] = 0 ; nor[o + 1 ] = 1 ; nor[o + 2 ] = 0 ; continue ; } const l = i > 0 && wet[k - 1 ] ? y[k - 1 ] : y[k], r = i < N - 1 && wet[k + 1 ] ? y[k + 1 ] : y[k]; const b = j > 0 && wet[k - N] ? y[k - N] : y[k], f = j < N - 1 && wet[k + N] ? y[k + N] : y[k]; const gx = (r - l) / d2, gz = (f - b) / d2, m = 1 / Math.sqrt(1 + gx * gx + gz * gz); nor[o] = -gx * m; nor[o + 1 ] = m; nor[o + 2 ] = -gz * m; } const at = (g.water.geometry as THREE.BufferGeometry).attributes; at.position.needsUpdate = true ; at.normal.needsUpdate = true ; at.aF.needsUpdate = true ; at.aE.needsUpdate = true ; // 断面の水: 上の辺は水面、下の辺は海底 const wp = (g.wall.geometry as THREE.BufferGeometry).attributes.position as THREE.BufferAttribute, w = wp.array as Float32Array; for (let e = 0 ; e < 4 ; e++) for (let n = 0 ; n < N; n++) { const k = edgeCell(N, e, n), o = (e * N + n) * 6 ; if (wet[k]) { w[o + 1 ] = y[k]; w[o + 4 ] = bed[k] * yb; } else { w[o + 1 ] = w[o + 4 ] = bed[k] * yb; } } wp.needsUpdate = true ; // 造波板: 板が前後に動き、となりの板と少しずつずれて、斜めの波を作る const wv = s.wave, amp = clamp(wv.a * 1.6 , 0 , 1.4 ); for (const m of s.pads) if (m.visible) m.position.z = ZP * W + 0.6 - amp * Math.cos(wv.ph - (wv.kx * m.position.x) / s.k); for (const b of s.buoys) b.m.position.set(b.fx * s.k, (b.k >= 0 && wet[b.k] ? eta[b.k] * yw : 0 ) + 0.2 , b.fz * s.k); } /** つかめる物の印: 点線の輪(地の色で縁取りして、明るい水面の上でも見えるように) */ function ring(api: Api, x: number, y: number, r: number) { api.circle(x, y, r, { w: api.u(4 ), stroke: "paper" }); api.circle(x, y, r, { w: api.u(1.6 ), dash: [api.u(3 ), api.u(3 )] }); } /** 3D の絵の上に重ねる文字と印: ブイの名前・波の向きの矢印と取っ手・物差し・なぞる範囲 */ function overlay(api: Api, s: S) { const v = s.v, U = api.u, k = s.k, sp = s.sp; if (sp.maker) { const a = arrowAt(s), p = v.project([a[0 ] * k, 0.6 , a[1 ] * k]), q = v.project([a[2 ] * k, 0.6 , a[3 ] * k]); if (p.on && q.on) { api.arrow(p.x, p.y, q.x, q.y, { w: U(6 ), head: U(15 ), stroke: "paper" }); api.arrow(p.x, p.y, q.x, q.y, { w: U(2.5 ), head: U(11 ) }); ring(api, q.x, q.y, U(15 )); api.label(` θ ${fmt((s.wave.th * 180) / Math.PI, 0)}°` , q.x + U(20 ), q.y + U(5 ), { size: U(13 ) }); } const span = sp.per ? 1 : 1 - 2 * SPONGE, c = v.project([-W * span * 0.42 , 2.2 , ZP * W + 0.6 ]); if (c.on) api.label("造波板" , c.x, c.y - U(6 ), { size: U(12 ), align: "center" }); } if (s.key === "tsunami" ) { // 岸へ向かう津波の山(左寄りの列で、海面がいちばん高い所) const N = s.N, i = N >> 2 ; let best = -1 , top = 0.05 ; for (let j = 0 ; j < N; j++) { const q = j * N + i; if (s.wet[q] && cx(s, j) < 0.2 * s.L && s.eta[q] > top) { top = s.eta[q]; best = j; } } if (best >= 0 ) { const p = v.project([cx(s, i) * k, top * sp.vw * k + 1 , cx(s, best) * k]); if (p.on) api.label(` 津波の山 ${fmt(top, 2)} m` , p.x, p.y - U(14 ), { size: U(13 ), align: "center" }); } } s.buoys.forEach((b, n) => { const p = v.project(b.m.position); if (!p.on) return ; ring(api, p.x, p.y, U(14 )); api.label("AB" [n], p.x + U(16 ), p.y - U(10 ), { size: U(14 ), weight: 700 }); }); // 物差し(手前の辺に沿って)と縦の倍率 const z = (W / 2 ) * ((s.N - 1 ) / s.N) + 1.5 , x0 = -W * 0.42 , x1 = x0 + sp.bar * k; const a = v.project([x0, 0 , z]), b = v.project([x1, 0 , z]); if (a.on && b.on) { for (const [w, c] of [[U(5 ), "paper" ], [U(2 ), "ink" ]] as const ) { api.line(a.x, a.y, b.x, b.y, { w, stroke: c }); api.line(a.x, a.y - U(5 ), a.x, a.y + U(5 ), { w, stroke: c }); api.line(b.x, b.y - U(5 ), b.x, b.y + U(5 ), { w, stroke: c }); } const vert = sp.vb === sp.vw ? ` 縦は ${sp.vb} 倍` : ` 縦は海底 ${sp.vb} 倍・波 ${sp.vw.toLocaleString("ja-JP")} 倍` ; const tx = clamp(Math.min(a.x, b.x), U(8 ), api.w - U(230 )), ty = Math.min(Math.max(a.y, b.y) + U(22 ), api.h - U(10 )); api.label(` ${sp.barText}(${vert})` , tx, ty, { size: U(12 ) }); } // 盛る・掘るの範囲 if (api.str("touch" ) !== "view" && api.pointer.inside && s.drag !== "A" && s.drag !== "B" && s.drag !== "arrow" ) { const g = ground(s, api.pointer.x, api.pointer.y); if (g) { const R = 0.05 * s.L * k * 1.6 , c = v.project([g.x * k, 0 , g.z * k]), e = v.project([g.x * k + R, 0 , g.z * k]); ring(api, c.x, c.y, Math.max(U(8 ), Math.hypot(e.x - c.x, e.y - c.y))); } } } function readout(api: Api, s: S) { const sp = s.sp, w = s.wave, ts = s.key === "tsunami" , lines: string[] = []; const dep = (h: number) => (h >= 100 ? Math.round(h).toLocaleString("ja-JP" ) : fmt(h, h < 10 ? 1 : 0 )); const vel = (c: number) => (ts ? ` ${Math.round(c * 3.6).toLocaleString("ja-JP")} km/h` : ` ${fmt(c, 1)} m/s` ); if (sp.maker) lines.push(` 沖 水深 ${dep(sp.h0)} m 速さ ${vel(w.c0)} 波長 ${fmt(w.c0 * w.T, 0)} m` ); else lines.push(` 沖 水深 ${dep(sp.h0)} m 速さ ${vel(Math.sqrt(G * sp.h0))}` ); s.buoys.forEach((b, n) => { const name = ` ブイ${"AB"[n]}` ; if (b.k < 0 || !s.wet[b.k]) { lines.push(` ${name} 陸の上` ); return ; } const h = s.h[b.k], c = Math.sqrt(G * h); if (ts) { lines.push(` ${name} 水深 ${dep(h)} m 速さ ${vel(c)} いちばん高い水位 ${fmt(Math.max(0, b.top), 2)} m` ); return ; } const d = dirOf(s, b); lines.push(` ${name} 水深 ${dep(h)} m 速さ ${vel(c)} 波高 ${fmt(heightOf(b), 2)} m 向き ${d === null ? "定まらない" : fmt(d, 0) + "°"}` ); }); const want = Number(api.str("speed" )) || 1 , got = api.playing && s.cap > 0.5 ? clamp(Math.round(s.spd), 1 , want) : want; const sp2 = ` 実際の ${got} 倍の速さ${got < want ? "(計算が追いつく上限)" : ""}` ; if (ts) { const m = Math.floor(s.time / 60 ), sec = Math.floor(s.time % 60 ); lines.push(` 地震から ${m}分${sec}秒(${sp2})` ); } else lines.push(` 時間 ${s.spin > 0 ? "波が岸に届くまで早送り中" : sp2}` ); api.readout(lines); } function notes(api: Api, s: S) { if (s.key !== "tsunami" && s.brk > 0 && !s.noted.brk && api.t > 4 ) { s.noted.brk = true ; api.note("浅い所で波が崩れています。波高が水深の 0.8 倍ほどになると、波は形を保てません" , 3400 ); } const B = s.buoys[1 ]; if (s.key === "tsunami" && !s.noted.shelf && B.top > 0.3 ) { s.noted.shelf = true ; api.note("大陸棚に入った津波は遅くなり、そのぶん高くなります" , 3400 ); } } 水面は毎コマ、頂点の高さと面の向きだけを書き換える。面の向きは、隣の点との高さの差（水面の傾き）から計算する。水面の材質は自前で、光の当たり方と照り返し、水面の高さ、砕けた所の泡、等深線を1つの明るさにまとめて出すだけである。その明るさを白と黒の網点に直す仕事と、輪郭線を引く仕事は、全記事に共通の部品が受け持つ。陸のますでは水面の点を海底より下へ沈めておく。こうすると、水面と陸の面が交わる所が、そのまま波打ち際になる。
縦の長さは引き伸ばして描いている。津波のほかの場面は海底も波も6倍、津波の場面は海底を10倍、波を12,000倍にした。津波の場面で水の下の海底を描かないのは、波だけを大きく引き伸ばすと、谷が海底を突き抜けて見えてしまうからだ。
