このシミュレーションは、1/240秒ごとに「触れそうな粒の組の表を確かめる → 組ごとの力を足す → 床と壁の力を足す → すべての粒の速さと位置を進める」をくり返している。計算はすべて自前の型付き配列の上で動き、three.jsは結果の位置を受け取って球を並べて描くだけだ。粒の物理にthree.jsの機能は使っていない。
長さはメートル、粒は直径20 cm
長さはメートル、時間は秒、質量はキログラム、力はニュートンで数える。粒の直径は平均0.2 mで、±15%の散らばりを持たせた。大きさがそろった球は結晶のように並び、山の形が変わってしまうからである。
定数と単位(src/sims/sand-pile.ts の 20〜44 行) // 長さ m・時間 秒・質量 kg・力 N(1 = 1 m) const G = 9.8 ; // 重力加速度 m/s² const D0 = 0.2 ; // 粒の直径の平均 m const SPREAD = 0.15 ; // 粒の大きさの散らばり ±15 %。そろった球は結晶のように並び、山の形が変わってしまう const RHO = 2650 ; // 粒の密度 kg/m³(石英) const M0 = RHO * Math.PI / 6 * D0 ** 3 ; // 平均の粒の質量 約 11 kg const KN = 40000 ; // 押し合うばね N/m。粒の重さ(約 109 N)で 2.7 mm 縮む const KT = KN * 2 / 7 ; // 接線のばね N/m const LOG_E = -Math.log(0.3 ); // 反発係数 0.3 の対数 const CN = 2 * (LOG_E / Math.hypot(Math.PI, LOG_E)) * Math.sqrt(KN * M0 / 2 ); // 法線の減衰 N·s/m const CT = 0.2 * Math.sqrt(KT * M0 / 2 ); // 接線の減衰 N·s/m const KR = KT * (D0 / 4 ) ** 2 ; // 転がりのばね N·m/rad const CR = 0.2 * Math.sqrt(KR * 0.4 * M0 * (D0 / 2 ) ** 2 ); // 転がりの減衰 N·m·s/rad const CW = 1200 ; // 壁の摩擦の、遅いときの粘り N·s/m const H = 1 / 240 ; // 時間の刻み 秒。1 コマ(1/60 秒)に 4 回 const SKIN = 0.1 ; // 近くの粒の表に入れる余裕 m const MU_FLOOR = 1.0 ; // 床と穴の開いた板の摩擦係数 const MUR_FLOOR = 0.3 ; // 床の転がり抵抗 const VFLOW = 0.15 ; // これより速い粒を「流れている」と数える m/s const VMAX = 30 ; // 速さの上限 m/s(計算が暴れないため) const BASE = [750 , 1500 , 3000 ]; // 画質 0・1・2 の粒の数 const GAP_WET = 0.01 ; // 水の橋が切れるすきま m const RHO_B = 0.6 * RHO; // ベバルー則に入れるかさ密度 kg/m³(粒が占める割合 0.6) const DEG = Math.PI / 180 ;const LAT_A = 0.235 , LAT_G = 0.17 , LAT_V = LAT_A * LAT_A * LAT_G; // 最初に粒を置く格子: 横の間隔・段の間隔(一段おきに半分ずらす)・粒 1 個あたりの体積 ばね定数KNは、ばねの揺れが刻みHで十分追える固さから決めた。粒の質量m k n m / k n CNは、反発係数0.3から求めた。
近くの粒だけを調べる
3,000粒のすべての組を毎回調べると、組の数は約450万になる。そこで「いつか触れるかもしれない組の表」を作っておき、力の計算ではこの表だけを見る。
近くの粒を探す(src/sims/sand-pile.ts の 139〜188 行) // 触れそうな組の表を、粒が表の余裕の半分だけ動くまで使い回す。表を作るときは 1 辺約 0.33 m の格子で近くのますだけを見る private rebuild() { const { P, R, alive, fixed, cellHead, cellNext, cellOf, chainHead, chainNext, cs, nx, ny, nz, ox, oy, oz } = this ; const n = this .n, PA = this .PA, PB = this .PB, HS = this .HS, FC = this .FC, A2 = this .PA2, B2 = this .PB2, H2 = this .HS2, C2 = this .FC2, maxP = this .maxP; cellHead.fill(-1 ); for (let i = 0 ; i < n; i++) { if (!alive[i]) continue ; const k = i * 3 ; let ix = ((P[k] - ox) / cs) | 0 , iy = ((P[k + 1 ] - oy) / cs) | 0 , iz = ((P[k + 2 ] - oz) / cs) | 0 ; ix = ix < 0 ? 0 : ix >= nx ? nx - 1 : ix; iy = iy < 0 ? 0 : iy >= ny ? ny - 1 : iy; iz = iz < 0 ? 0 : iz >= nz ? nz - 1 : iz; const c = ix + nx * (iy + ny * iz); cellOf[i] = c; cellNext[i] = cellHead[c]; cellHead[c] = i; } // 触れていた組の、ばねにたまった伸びを引き継ぐ。小さいほうの番号ごとにつないでおき、新しい表を作るときに探す chainHead.fill(-1 , 0 , n); for (let p = 0 , e = this .np; p < e; p++) if (FC[p] !== 0 ) { const a = PA[p]; chainNext[p] = chainHead[a]; chainHead[a] = p; } let m = 0 ; for (let i = 0 ; i < n; i++) { if (!alive[i]) continue ; const ki = i * 3 , xi = P[ki], yi = P[ki + 1 ], zi = P[ki + 2 ], ri = R[i], c = cellOf[i], fi = fixed[i]; const iz = (c / (nx * ny)) | 0 , rest = c - iz * nx * ny, iy = (rest / nx) | 0 , ix = rest - iy * nx; // 自分のますの残りと、前半分の 13 のますだけを見る(同じ組を二度数えない) for (let s = 0 ; s < 14 ; s++) { let j: number; if (s === 0 ) j = cellNext[i]; else { const jx = ix + STX[s], jy = iy + STY[s], jz = iz + STZ[s]; if (jx < 0 || jx >= nx || jy < 0 || jy >= ny || jz < 0 || jz >= nz) continue ; j = cellHead[jx + nx * (jy + ny * jz)]; } for (; j >= 0 ; j = cellNext[j]) { const kj = j * 3 , dx = P[kj] - xi, dy = P[kj + 1 ] - yi, dz = P[kj + 2 ] - zi, rr = ri + R[j] + SKIN; if (dx * dx + dy * dy + dz * dz >= rr * rr || m >= maxP || (fi && fixed[j])) continue ; const a = i < j ? i : j, b = i < j ? j : i, q = m * 6 ; A2[m] = a; B2[m] = b; C2[m] = 0 ; H2[q] = 0 ; H2[q + 1 ] = 0 ; H2[q + 2 ] = 0 ; H2[q + 3 ] = 0 ; H2[q + 4 ] = 0 ; H2[q + 5 ] = 0 ; for (let o = chainHead[a]; o >= 0 ; o = chainNext[o]) if (PB[o] === b) { const u = o * 6 ; H2[q] = HS[u]; H2[q + 1 ] = HS[u + 1 ]; H2[q + 2 ] = HS[u + 2 ]; H2[q + 3 ] = HS[u + 3 ]; H2[q + 4 ] = HS[u + 4 ]; H2[q + 5 ] = HS[u + 5 ]; C2[m] = FC[o]; break ; } m++; } } } this .np = m; this .PA = A2; this .PB = B2; this .HS = H2; this .FC = C2; this .PA2 = PA; this .PB2 = PB; this .HS2 = HS; this .FC2 = FC; this .disp = 0 ; this .dirty = false ; } 表を作るときは、空間を1辺約0.33 mのますに分け、各粒が自分のますと隣のますの粒だけを調べる。隣のますを半分の13個に絞っているのは、同じ組を二度数えないためだ。表には、中心の距離が半径の和より0.1 m長い組まで入れる。この余裕の8割だけ粒が動くまで、同じ表を使い回す。止まった山なら表の作り直しはほとんど起きない。
表を作り直すときに忘れてはならないのが、ばねにたまった伸びである。接線と転がりのばねは組ごとの記憶なので、触れていた組の伸びを新しい表へ写している。
粒どうしの力
本文の3つの式が、この部分にそのまま入っている。
粒どうしの力(src/sims/sand-pile.ts の 190〜262 行) private contacts() { const { P, V, W, F, T, R, PA, PB, HS, FC } = this ; const mu = this .mu, mur = this .mur, coh = this .coh, gap = this .gapR; let omax = 0 ; for (let p = 0 , np = this .np; p < np; p++) { const a = PA[p], b = PB[p], ia = a * 3 , ib = b * 3 ; const dx = P[ib] - P[ia], dy = P[ib + 1 ] - P[ia + 1 ], dz = P[ib + 2 ] - P[ia + 2 ]; const ra = R[a], rb = R[b], sum = ra + rb, reach = sum + gap, d2 = dx * dx + dy * dy + dz * dz; if (d2 >= reach * reach) { if (FC[p] !== 0 ) { const q = p * 6 ; FC[p] = 0 ; HS[q] = 0 ; HS[q + 1 ] = 0 ; HS[q + 2 ] = 0 ; HS[q + 3 ] = 0 ; HS[q + 4 ] = 0 ; HS[q + 5 ] = 0 ; } continue ; } const dist = Math.sqrt(d2) + 1e-9 , inv = 1 / dist, nx = dx * inv, ny = dy * inv, nz = dz * inv, over = sum - dist; if (over <= 0 ) { // 触れる手前: 水の橋の引力だけが働く(離れるほど弱まり、すきま 1 cm で切れる)。乾いた砂(gap = 0)では 0 で割らないよう何もしない if (FC[p] !== 0 ) { const q = p * 6 ; FC[p] = 0 ; HS[q] = 0 ; HS[q + 1 ] = 0 ; HS[q + 2 ] = 0 ; HS[q + 3 ] = 0 ; HS[q + 4 ] = 0 ; HS[q + 5 ] = 0 ; } if (gap > 0 ) { const f = coh * (1 + over / gap); F[ia] += f * nx; F[ia + 1 ] += f * ny; F[ia + 2 ] += f * nz; F[ib] -= f * nx; F[ib + 1 ] -= f * ny; F[ib + 2 ] -= f * nz; } continue ; } if (over > omax) omax = over; // 接点での相対速度(b から a を引いた向き) const rvx = V[ib] - V[ia], rvy = V[ib + 1 ] - V[ia + 1 ], rvz = V[ib + 2 ] - V[ia + 2 ]; const vn = rvx * nx + rvy * ny + rvz * nz; // 押し合う力 N = k_n δ − c_n v_n(引っ張りにはならない) let fn = KN * over - CN * vn; if (fn < 0 ) fn = 0 ; FC[p] = fn > 1e-6 ? fn : 1e-6 ; const wax = W[ia], way = W[ia + 1 ], waz = W[ia + 2 ], wbx = W[ib], wby = W[ib + 1 ], wbz = W[ib + 2 ]; const sx = ra * wax + rb * wbx, sy = ra * way + rb * wby, sz = ra * waz + rb * wbz; const vtx = rvx - vn * nx - (sy * nz - sz * ny), vty = rvy - vn * ny - (sz * nx - sx * nz), vtz = rvz - vn * nz - (sx * ny - sy * nx); // 接線のばね: 接点が横にずれた長さをためる。ためた向きは今の接平面の上に入れ直す const q = p * 6 ; let ux = HS[q], uy = HS[q + 1 ], uz = HS[q + 2 ]; const un = ux * nx + uy * ny + uz * nz; ux += vtx * H - un * nx; uy += vty * H - un * ny; uz += vtz * H - un * nz; let ftx = -KT * ux, fty = -KT * uy, ftz = -KT * uz; const lim = mu * fn, ft2 = ftx * ftx + fty * fty + ftz * ftz; if (ft2 > lim * lim) { // クーロンの上限 μN を超えたらすべる。ばねの伸びは上限の分まで縮める const c = lim / Math.sqrt(ft2); ftx *= c; fty *= c; ftz *= c; ux = -ftx / KT; uy = -fty / KT; uz = -ftz / KT; } else { ftx -= CT * vtx; fty -= CT * vty; ftz -= CT * vtz; } HS[q] = ux; HS[q + 1 ] = uy; HS[q + 2 ] = uz; // 粒 b が受ける力 = (押し合い − 水の引力)の法線の力 + 接線の力。粒 a はその逆(作用と反作用) const fnn = fn - coh, fx = fnn * nx + ftx, fy = fnn * ny + fty, fz = fnn * nz + ftz; F[ib] += fx; F[ib + 1 ] += fy; F[ib + 2 ] += fz; F[ia] -= fx; F[ia + 1 ] -= fy; F[ia + 2 ] -= fz; const cx = ny * ftz - nz * fty, cy = nz * ftx - nx * ftz, cz = nx * fty - ny * ftx; T[ia] -= ra * cx; T[ia + 1 ] -= ra * cy; T[ia + 2 ] -= ra * cz; T[ib] -= rb * cx; T[ib + 1 ] -= rb * cy; T[ib + 2 ] -= rb * cz; // 転がりの抵抗: 転がり合った角度をためる。上限 μr R* N を超えたら転がる if (mur > 0 ) { let rx = wax - wbx, ry = way - wby, rz = waz - wbz; const rn = rx * nx + ry * ny + rz * nz; rx -= rn * nx; ry -= rn * ny; rz -= rn * nz; let tx = HS[q + 3 ], ty = HS[q + 4 ], tz = HS[q + 5 ]; const tn = tx * nx + ty * ny + tz * nz; tx += rx * H - tn * nx; ty += ry * H - tn * ny; tz += rz * H - tn * nz; let mx = -KR * tx, my = -KR * ty, mz = -KR * tz; const mlim = mur * (ra * rb / sum) * fn, m2 = mx * mx + my * my + mz * mz; if (m2 > mlim * mlim) { const c = mlim / Math.sqrt(m2); mx *= c; my *= c; mz *= c; tx = -mx / KR; ty = -my / KR; tz = -mz / KR; } else { mx -= CR * rx; my -= CR * ry; mz -= CR * rz; } HS[q + 3 ] = tx; HS[q + 4 ] = ty; HS[q + 5 ] = tz; T[ia] += mx; T[ia + 1 ] += my; T[ia + 2 ] += mz; T[ib] -= mx; T[ib + 1 ] -= my; T[ib + 2 ] -= mz; } } this .overMax = omax; } まず2つの粒の距離からめり込みδ N = k n δ − c n v n sx、sy、szの項で足している。
接線のばねは、横にずれた長さux、uy、uzをためていく。接点の向きは粒の動きとともに変わるので、ためた向きを今の接平面に入れ直してから足す。力がμ N μ N ∣ F t ∣ ≤ μ N
力は粒bに足し、粒aには逆向きに足す。作用と反作用である。接線の力は粒を回す力にもなるので、半径を掛けて回転の力Tにも足す。濡れた砂では、押し合う力から水の引力を引くだけで、接点が離れても1 cmまでは引力が残る。摩擦の上限は押し合う力で決まるので、水が粒を押し付ける分だけ、摩擦も強くなる。
床と壁
床と壁(src/sims/sand-pile.ts の 264〜340 行) /** 壁との接触: n は壁から粒へ向かう向き、over はめり込み。摩擦はばねを持たず、遅いときは粘りで止める */ private wall(i: number, nx: number, ny: number, nz: number, over: number, mu: number, wvx = 0 , wvy = 0 , wvz = 0 ) { const { V, W, F, T, R } = this ; const k = i * 3 , r = R[i], wx = W[k], wy = W[k + 1 ], wz = W[k + 2 ]; const vx = V[k] - wvx - r * (wy * nz - wz * ny), vy = V[k + 1 ] - wvy - r * (wz * nx - wx * nz), vz = V[k + 2 ] - wvz - r * (wx * ny - wy * nx); const vn = vx * nx + vy * ny + vz * nz; const fn = KN * over - CN * vn; if (fn <= 0 ) return ; let fx = fn * nx, fy = fn * ny, fz = fn * nz; const tx = vx - vn * nx, ty = vy - vn * ny, tz = vz - vn * nz, ts = Math.sqrt(tx * tx + ty * ty + tz * tz); if (ts > 1e-9 && mu > 0 ) { const lim = mu * fn, ft = CW * ts < lim ? CW * ts : lim, c = -ft / ts, ftx = c * tx, fty = c * ty, ftz = c * tz; fx += ftx; fy += fty; fz += ftz; T[k] -= r * (ny * ftz - nz * fty); T[k + 1 ] -= r * (nz * ftx - nx * ftz); T[k + 2 ] -= r * (nx * fty - ny * ftx); } F[k] += fx; F[k + 1 ] += fy; F[k + 2 ] += fz; } private boundaries() { const { P, V, W, F, T, R, alive, fixed, onFloor, HF } = this ; const n = this .n, krf = KR * 4 , crf = CR * 2 ; const xMin = this .xMin, xMax = this .xMax, zMin = this .zMin, zMax = this .zMax, muW = this .muW, silo = this .siloR, plate = this .hasPlate, py = this .plateY, hr = this .holeR, damOn = this .damOn, damX = this .damX, damEdge = this .damEdge; const fy0 = this .floorY, fx0 = this .fx0, fx1 = this .fx1, fz0 = this .fz0, fz1 = this .fz1, tool = this .toolOn, tr = this .tr, tx = this .tx, ty = this .ty, tz = this .tz; for (let i = 0 ; i < n; i++) { if (!alive[i] || fixed[i]) continue ; const k = i * 3 , x = P[k], y = P[k + 1 ], z = P[k + 2 ], r = R[i]; // 床(穴の開いた板も同じ)。平らな面なので、接線と転がりのばねの伸びを粒ごとに持つ let ground = fy0, flat = x >= fx0 && x <= fx1 && z >= fz0 && z <= fz1; if (plate && y > py - r * 0.5 ) { const rho2 = x * x + z * z; if (rho2 >= hr * hr) { ground = py; flat = true ; } else { flat = false ; // 穴の縁(輪)に当たる const rho = Math.sqrt(rho2) + 1e-9 , e = hr - rho, dy = y - py, dd = e * e + dy * dy; if (dd < r * r) { const dist = Math.sqrt(dd) + 1e-9 ; this .wall(i, -x / rho * e / dist, dy / dist, -z / rho * e / dist, r - dist, muW); } } } const over = ground + r - y; if (flat && over > 0 ) { const fn0 = KN * over - CN * V[k + 1 ], fn = fn0 > 0 ? fn0 : 0 , q = i * 4 , wx = W[k], wz = W[k + 2 ]; const vcx = V[k] + r * wz, vcz = V[k + 2 ] - r * wx; let ux = HF[q] + vcx * H, uz = HF[q + 1 ] + vcz * H; let fx = -KT * ux, fz = -KT * uz; const lim = MU_FLOOR * fn, f2 = fx * fx + fz * fz; if (f2 > lim * lim) { const c = lim / Math.sqrt(f2); fx *= c; fz *= c; ux = -fx / KT; uz = -fz / KT; } else { fx -= CT * vcx; fz -= CT * vcz; } HF[q] = ux; HF[q + 1 ] = uz; F[k] += fx; F[k + 1 ] += fn; F[k + 2 ] += fz; T[k] -= r * fz; T[k + 2 ] += r * fx; let ax = HF[q + 2 ] + wx * H, az = HF[q + 3 ] + wz * H; let mx = -krf * ax, mz = -krf * az; const ml = MUR_FLOOR * r * fn, m2 = mx * mx + mz * mz; if (m2 > ml * ml) { const c = ml / Math.sqrt(m2); mx *= c; mz *= c; ax = -mx / krf; az = -mz / krf; } else { mx -= crf * wx; mz -= crf * wz; } HF[q + 2 ] = ax; HF[q + 3 ] = az; T[k] += mx; T[k + 2 ] += mz; onFloor[i] = 1 ; } else if (onFloor[i]) { onFloor[i] = 0 ; const q = i * 4 ; HF[q] = 0 ; HF[q + 1 ] = 0 ; HF[q + 2 ] = 0 ; HF[q + 3 ] = 0 ; } // 容器の壁 if (x - r < xMin) this .wall(i, 1 , 0 , 0 , xMin - x + r, muW); if (x + r > xMax) this .wall(i, -1 , 0 , 0 , x + r - xMax, muW); if (z - r < zMin) this .wall(i, 0 , 0 , 1 , zMin - z + r, muW); if (z + r > zMax) this .wall(i, 0 , 0 , -1 , z + r - zMax, muW); // 引き抜く壁: 引き上げた縁より上の粒だけを止める if (damOn && x + r > damX && y > damEdge) this .wall(i, -1 , 0 , 0 , x + r - damX, muW); if (silo > 0 && y > py - r) { const rho2 = x * x + z * z, lim = silo - r; if (rho2 > lim * lim) { const rho = Math.sqrt(rho2); this .wall(i, -x / rho, 0 , -z / rho, rho - lim, muW); } } // 指の玉に押される if (tool) { const dx = x - tx, dy = y - ty, dz = z - tz, rr = tr + r, d2 = dx * dx + dy * dy + dz * dz; if (d2 < rr * rr && d2 > 1e-12 ) { const d = Math.sqrt(d2); this .wall(i, dx / d, dy / d, dz / d, rr - d, 0.3 , this .tvx, this .tvy, this .tvz); } } } } 床は平らなので、粒どうしと同じばねの伸びを粒ごとに持つ。穴の開いた板では、中心が穴の中にある粒だけ、穴の縁（輪）との距離で当たりを調べる。壁の摩擦はばねを持たない簡単な形で、遅いときは粘りで止める。
時間を進める
時間を進める(src/sims/sand-pile.ts の 342〜391 行) // 1 刻みの流れ: 表を確かめる → 力を 0 にする → 粒どうしの力 → 床と壁の力 → 速さと位置を進める private integrate() { const { P, V, W, F, T, IM, II, alive, fixed } = this ; const n = this .n, h = H, gx = this .gx, gy = this .gy, killY = this .killY, vf2 = VFLOW * VFLOW; // 格子に置いた粒を詰まらせる最初の間だけ、空気の抵抗のように速さを落とす(勢いよく詰まって、ばねが縮みすぎないように) const drag = this .t < this .settle ? 0.97 : 1 ; let vm2 = 0 , mov = 0 , cnt = 0 ; for (let i = 0 ; i < n; i++) { if (!alive[i] || fixed[i]) continue ; const k = i * 3 , im = IM[i], ii = II[i]; // 半陰的オイラー法: 新しい速さで位置を進める let vx = (V[k] + (F[k] * im + gx) * h) * drag, vy = (V[k + 1 ] + (F[k + 1 ] * im + gy) * h) * drag, vz = (V[k + 2 ] + F[k + 2 ] * im * h) * drag; let v2 = vx * vx + vy * vy + vz * vz; if (!(v2 < VMAX * VMAX)) { if (v2 !== v2) { this .kill(i); continue ; } const c = VMAX / Math.sqrt(v2); vx *= c; vy *= c; vz *= c; v2 = VMAX * VMAX; } V[k] = vx; V[k + 1 ] = vy; V[k + 2 ] = vz; W[k] += T[k] * ii * h; W[k + 1 ] += T[k + 1 ] * ii * h; W[k + 2 ] += T[k + 2 ] * ii * h; P[k] += vx * h; P[k + 1 ] += vy * h; P[k + 2 ] += vz * h; if (v2 > vm2) vm2 = v2; if (v2 > vf2) mov++; cnt++; if (P[k + 1 ] < killY) this .kill(i); } this .moving = cnt > 0 ? mov / cnt : 0 ; this .vmax = Math.sqrt(vm2); this .disp += this .vmax * h; } step() { if (!this .sleeping) { // 粒が表の余裕の 8 割だけ動いたら表を作り直す(落ちてくる粒が止まっている粒に近づく速さは、落ちる速さと同じ) if (this .dirty || this .disp > (SKIN - this .gapR) * 0.8 ) this .rebuild(); const e = this .n * 3 ; this .F.fill(0 , 0 , e); this .T.fill(0 , 0 , e); this .contacts(); this .boundaries(); this .integrate(); // すべての粒が止まって 1 秒たったら、計算を休む(何かが変わるまで結果は同じなので)。 // 壁の摩擦はばねを持たないので、壁のそばの粒はわずかに回り続ける。水路と箱では少し大きな速さまで止まったとみなす const calm = this .key === "slump" || this .key === "tilt" ? 0.08 : 0.02 ; this .still = this .vmax < calm && !this .pouring && !this .toolOn ? this .still + H : 0 ; if (this .still > 1 ) { this .sleeping = true ; this .moving = 0 ; } } this .stepN++; this .t += H; this .sceneStep(); } /** 止まって休んでいる計算を起こす(材料や容器を変えたとき) */ wake() { this .sleeping = false ; this .still = 0 ; } 半陰的オイラー法で進める。先に力から速さを更新し、その新しい速さで位置を動かす。ばねの揺れのエネルギーが増えにくい、単純で扱いやすい方法である。回転の刻みも同じだが、小さな粒は慣性モーメントが半径の5乗で小さくなるので、接線と転がりの減衰を弱めてある。強すぎると、減衰の計算そのものが揺れを増やしてしまう。
すべての粒が止まって1秒たったら、力の計算を休む。何かが変わるまで結果は同じだからだ。つまみを動かしたり砂を注いだりすると、wakeで計算を起こす。水路と箱では、壁のそばの粒がわずかに回り続けるので、止まったとみなす速さを少し大きくしてある。
格子に置いた粒を詰まらせる最初の間だけ、空気の抵抗のように速さを落としている。勢いよく詰まると、柔らかいばねが縮みすぎて、砂が固まってしまうためである。
山の形から角度を測る
山を測る(src/sims/sand-pile.ts の 528〜650 行) /** 高さの表: 1 辺 20 cm のますごとに、止まっている粒のてっぺんの高さ(落ちている粒は入れない) */ surface(y0 = this .floorY, onlyAbove = -Infinity) { const { P, V, R, alive, top, tcs, tnx, tnz, tx0, tz0 } = this ; top.fill(y0); for (let i = 0 ; i < this .n; i++) { if (!alive[i]) continue ; const k = i * 3 ; if (P[k + 1 ] < onlyAbove || V[k] * V[k] + V[k + 1 ] * V[k + 1 ] + V[k + 2 ] * V[k + 2 ] > 1 ) continue ; const a = ((P[k] - tx0) / tcs) | 0 , b = ((P[k + 2 ] - tz0) / tcs) | 0 ; if (a < 0 || b < 0 || a >= tnx || b >= tnz) continue ; const c = a + tnx * b, h = P[k + 1 ] + R[i]; if (h > top[c]) top[c] = h; } } /** 山の形: 中心からの距離ごとに高さをならし、高さが 25〜80 % の所に直線を当てはめた傾きを、斜面の角度とする */ pileAngle() { this .surface(); const { top, tcs, tnx, tnz, tx0, tz0, bin, binN, m } = this , y0 = this .floorY; let hm = 0 ; for (let i = 0 ; i < top.length; i++) if (top[i] - y0 > hm) hm = top[i] - y0; m.top = hm; m.angle = NaN; if (hm < 0.5 ) return ; let sw = 0 , cx = 0 , cz = 0 ; for (let b = 0 ; b < tnz; b++) for (let a = 0 ; a < tnx; a++) { const h = top[a + tnx * b] - y0; if (h > 0.8 * hm) { sw += h; cx += h * (tx0 + (a + 0.5 ) * tcs); cz += h * (tz0 + (b + 0.5 ) * tcs); } } cx /= sw; cz /= sw; m.cx = cx; m.cz = cz; bin.fill(0 ); binN.fill(0 ); for (let b = 0 ; b < tnz; b++) for (let a = 0 ; a < tnx; a++) { const r = Math.hypot(tx0 + (a + 0.5 ) * tcs - cx, tz0 + (b + 0.5 ) * tcs - cz), k = (r / tcs) | 0 ; if (k < NBIN) { bin[k] += top[a + tnx * b] - y0; binN[k]++; } } let peak = 0 ; for (let k = 0 ; k < 6 ; k++) if (binN[k]) { const v = bin[k] / binN[k]; if (v > peak) peak = v; } this .fit(0 , NBIN, 0.25 * peak, 0.8 * peak, (k) => (k + 0.5 ) * tcs, -1 ); } /** 区切りごとの高さに直線を当てはめる。m.x1〜x2 に両端の位置、m.y1・y2 に高さを入れ、傾きを返す */ private fit(k0: number, k1: number, lo: number, hi: number, pos: (k: number) => number, sign: number) { const { bin, binN, m } = this ; let sx = 0 , sy = 0 , sxx = 0 , sxy = 0 , c = 0 , a = NaN, b = NaN; for (let k = k0; k < k1; k++) { if (!binN[k]) continue ; const v = bin[k] / binN[k], x = pos(k); if (v < lo) break ; if (v > hi) continue ; if (c === 0 ) a = x; b = x; sx += x; sy += v; sxx += x * x; sxy += x * v; c++; } if (c < 3 ) { m.angle = NaN; return 0 ; } const sl = (c * sxy - sx * sy) / (c * sxx - sx * sx), y0 = (sy - sl * sx) / c; m.angle = sign * Math.atan(sl) / DEG; m.x1 = a; m.y1 = y0 + sl * a; m.x2 = b; m.y2 = y0 + sl * b; return sl; } /** 水路の形: 長さの向きに高さをならす。場面「壁を抜いて崩す」で使う */ private channelProfile() { this .surface(); const { top, tcs, tnx, tnz, tx0, tz0, bin, binN, m } = this , y0 = this .floorY, W = this .WZ; bin.fill(0 ); binN.fill(0 ); for (let b = 0 ; b < tnz; b++) { const z = tz0 + (b + 0.5 ) * tcs; if (Math.abs(z) > W) continue ; for (let a = 0 ; a < tnx && a < NBIN; a++) { bin[a] += top[a + tnx * b] - y0; binN[a]++; } } let peak = 0 , kp = 0 , front = -this .L0; for (let k = 0 ; k < tnx && k < NBIN; k++) if (binN[k]) { const v = bin[k] / binN[k]; if (v > peak) { peak = v; kp = k; } if (v > 0.1 ) front = tx0 + (k + 1 ) * tcs; } m.top = peak; m.front = front; this .fit(kp, Math.min(tnx, NBIN), 0.2 * peak, 0.8 * peak, (k) => tx0 + (k + 0.5 ) * tcs, -1 ); } /** 容器の中の形: 容器の座標で直線を当てはめ、容器の傾きを足して水平からの角度に直す */ trayProfile() { this .surface(); const { top, tcs, tnx, tnz, tx0, tz0, bin, binN, m } = this , y0 = this .floorY, W = this .WZ, L = this .LX; bin.fill(0 ); binN.fill(0 ); let hm = 0 ; for (let b = 0 ; b < tnz; b++) { const z = tz0 + (b + 0.5 ) * tcs; if (Math.abs(z) > W) continue ; for (let a = 0 ; a < tnx && a < NBIN; a++) { const v = top[a + tnx * b] - y0; bin[a] += v; binN[a]++; if (v > hm) hm = v; } } m.top = hm; // 両端の壁のそばを除いた中ほどだけを見る const k0 = Math.ceil((-L + 0.15 * 2 * L - tx0) / tcs), k1 = Math.floor((L - 0.15 * 2 * L - tx0) / tcs); const sl = this .fit(k0, k1, -1 , 1e9 , (k) => tx0 + (k + 0.5 ) * tcs, 1 ); if (Number.isFinite(m.angle)) m.angle = this .phi / DEG - Math.atan(sl) / DEG; } private siloMeasure() { const { P, IM, alive, m } = this ; let mass = 0 ; for (let i = 0 ; i < this .n; i++) if (alive[i] && P[i * 3 + 1 ] > -0.1 ) mass += 1 / IM[i]; this .above = mass; // 2 秒前からの減り方で流量を求める const L = this .aboveLog, N = L.length; L[this .logN % N] = mass; this .logN++; const old = this .logN >= N ? L[this .logN % N] : L[0 ], span = (Math.min(this .logN, N) - 1 ) * 0.25 ; m.flow = !this .plugged && span > 0 ? Math.max(0 , old - mass) / span : 0 ; this .surface(this .plateY, this .plateY); const { top, tcs, tnx, tnz, tx0, tz0 } = this , Rs = this .siloR * 0.8 ; let sum = 0 , c = 0 ; for (let b = 0 ; b < tnz; b++) for (let a = 0 ; a < tnx; a++) { const x = tx0 + (a + 0.5 ) * tcs, z = tz0 + (b + 0.5 ) * tcs; if (x * x + z * z < Rs * Rs) { sum += top[a + tnx * b] - this .plateY; c++; } } m.fill = c ? sum / c : 0 ; if (!this .plugged && !this .jam && !this .done && this .opened >= 0 && this .t > this .opened + 3 && m.flow < 5 && this .moving < 0.02 ) { // 砂の高さが、穴の縁から筒の壁まで斜面で積もった形より高ければ、粒が橋を組んで詰まった。低ければ、流れ終わって穴のまわりに山が残っただけ if (m.fill > 0.8 * (this .siloR - this .holeR)) { this .jam = true ; this .msg = "穴の上で粒が橋を組み、流れが止まりました" ; } else { this .done = true ; this .msg = "流れ終わりました。穴のまわりの砂は、斜面の形で残ります" ; } } } /** 0.25 秒ごとに呼ぶ */ measure() { const m = this .m; if (this .sleeping && this .measured) return ; this .measured = this .sleeping; if (this .key === "pour" || this .key === "wet" || this .key === "chains" ) { this .pileAngle(); if (this .pouring) { const h = this .topAt(this .sx, this .sz); this .sy = Math.max(0.9 , h + 0.6 ); } } else if (this .key === "slump" ) this .channelProfile(); else if (this .key === "tilt" ) this .trayProfile(); else this .siloMeasure(); m.over = this .overMax; } /** ベバルー則: 穴から流れる質量 W = C ρ_b √g (D − k d)^(5/2) kg/s。C = 0.58・k = 1.5 */ beverloo(D = this .holeD) { const x = D - 1.5 * D0; return x > 0 ? 0.58 * RHO_B * Math.sqrt(G) * x ** 2.5 : 0 ; } 20 cm四方のますごとに、止まっている粒のてっぺんの高さを持つ。山なら、いちばん高い所からの距離ごとに高さをならし、高さが25〜80%の範囲に直線を当てはめ、その傾きを角度とする。てっぺんの丸みと、ふもとに散らばった粒を外すためである。画面に引く線も、この直線そのものだ。
箱を傾ける場面では、重力の向きのほうを回している。箱の中の座標で計算すれば壁も床も軸に沿ったままなので、当たりの計算を変えずに済む。測った傾きに箱の傾きを足し引きして、水平からの角度に直している。
three.jsで描く
three.jsで描く(src/sims/sand-pile.ts の 660〜823 行) const SP_V2 = 0.5 * 0.5 ; // この速さ(m/s)で粒をいちばん濃く描く const MAX_RODS = 1400 ;interface S { v: View3; sim: Sand; step: ReturnType<typeof stepper>; key: SceneKey; group: THREE.Group; deco: THREE.Group; world: THREE.Group; grains: THREE.InstancedMesh; rods: THREE.InstancedMesh; spout: THREE.Mesh; tool: THREE.Mesh; dam: THREE.Mesh | null ; plateMesh: THREE.Mesh | null ; grab: "" | "spout" | "dig" | "view" | "tilt" | "dam" ; px0: number; py0: number; phi0: number; moved: boolean; tNext: number; hist: Float32Array; histN: number; } const CUT = new THREE.Vector3(), CAM = new THREE.Vector3();const isPile = (k: SceneKey) => k === "pour" || k === "wet" || k === "chains" ; function makeGrains(s: { v: View3; group: THREE.Group }, cap: number) { // 三角形 80 個の球(各面を 4 つに分けた正二十面体)。濃さは粒ごとに変えるので、材質は濃さ 1 にして、粒の色で掛ける const mesh = new THREE.InstancedMesh(new THREE.IcosahedronGeometry(1 , 1 ), s.v.mat(1 ), cap); mesh.instanceColor = new THREE.InstancedBufferAttribute(new Float32Array(cap * 3 ), 3 ); mesh.instanceColor.setUsage(THREE.DynamicDrawUsage); mesh.instanceMatrix.setUsage(THREE.DynamicDrawUsage); const A = mesh.instanceMatrix.array as Float32Array; A.fill(0 ); for (let i = 0 ; i < cap; i++) A[i * 16 + 15 ] = 1 ; mesh.frustumCulled = false ; mesh.castShadow = true ; mesh.receiveShadow = true ; mesh.count = 0 ; s.group.add(mesh); return mesh; } function makeRods(s: { v: View3; group: THREE.Group }) { // 力の鎖の棒。粒の奥に隠れないよう、奥行きを無視して最後に描く const m = s.v.flat(1 ); m.depthTest = false ; m.depthWrite = false ; const mesh = new THREE.InstancedMesh(new THREE.CylinderGeometry(1 , 1 , 1 , 5 , 1 , true ), m, MAX_RODS); mesh.frustumCulled = false ; mesh.renderOrder = 5 ; mesh.count = 0 ; mesh.instanceMatrix.setUsage(THREE.DynamicDrawUsage); s.group.add(mesh); return mesh; } /** 粒の位置と濃さを、描画用の配列へ写す(毎コマ。部品は作り直さない) */ function syncGrains(api: Api, s: S) { const sim = s.sim, A = s.grains.instanceMatrix.array as Float32Array, C = s.grains.instanceColor!.array as Float32Array; const { P, V, R, alive, fixed } = sim, chain = api.bool("chain" ); // 断面の面: 注目点を通り、カメラの水平の向きを向いた面。手前の粒を消す let cx = 0 , cz = 0 , nx = 0 , nz = 0 ; if (chain) { s.group.updateWorldMatrix(true , false ); CAM.copy(s.v.camera.position); s.group.worldToLocal(CAM); CUT.set(isPile(s.key) ? sim.m.cx : s.key === "slump" ? 2 : 0 , 0 , isPile(s.key) ? sim.m.cz : 0 ); nx = CAM.x - CUT.x; nz = CAM.z - CUT.z; const l = Math.hypot(nx, nz) || 1 ; nx /= l; nz /= l; cx = CUT.x; cz = CUT.z; } const rest = chain ? 0.3 : 0.8 ; for (let i = 0 ; i < sim.n; i++) { const o = i * 16 , k = i * 3 ; if (!alive[i]) { A[o] = 0 ; A[o + 5 ] = 0 ; A[o + 10 ] = 0 ; continue ; } const x = P[k], y = P[k + 1 ], z = P[k + 2 ]; let r = R[i]; if (chain && (x - cx) * nx + (z - cz) * nz > 0 ) r = 0 ; A[o] = r; A[o + 5 ] = r; A[o + 10 ] = r; A[o + 12 ] = x; A[o + 13 ] = y; A[o + 14 ] = z; const v2 = V[k] * V[k] + V[k + 1 ] * V[k + 1 ] + V[k + 2 ] * V[k + 2 ]; // 止まっている粒は薄く、流れている粒は濃く。床に固定した粒はさらに薄く const t = fixed[i] ? rest * 0.7 : v2 >= SP_V2 ? 1.6 : rest + (1.6 - rest) * v2 / SP_V2; C[k] = t; C[k + 1 ] = t; C[k + 2 ] = t; } s.grains.count = sim.n; s.grains.instanceMatrix.needsUpdate = true ; s.grains.instanceColor!.needsUpdate = true ; // 力の鎖: 断面のすぐ奥(粒 3 個分)の、平均の 1.6 倍より大きな力の接触を棒で描く let nr = 0 ; if (chain) { const R2 = s.rods.instanceMatrix.array as Float32Array, { PA, PB, FC, np } = sim, slab = 3 * D0; let sum = 0 , c = 0 ; for (let p = 0 ; p < np; p++) if (FC[p] > 0 ) { sum += FC[p]; c++; } const mean = c ? sum / c : 1 ; for (let p = 0 ; p < np && nr < MAX_RODS; p++) { const f = FC[p]; if (f < mean * 1.6 ) continue ; const a = PA[p] * 3 , b = PB[p] * 3 , ax = P[a], ay = P[a + 1 ], az = P[a + 2 ], bx = P[b], by = P[b + 1 ], bz = P[b + 2 ]; const dep = ((ax + bx) / 2 - cx) * nx + ((az + bz) / 2 - cz) * nz; if (dep > 0 || dep < -slab) continue ; let dx = bx - ax, dy = by - ay, dz = bz - az; const len = Math.sqrt(dx * dx + dy * dy + dz * dz) || 1 ; dx /= len; dy /= len; dz /= len; // 棒の向きの行列: y 軸を接触の向きに、x・z 軸をそれと直角に。太さは力に比例させる let ux = dy, uy = -dx, uz = 0 ; const ul = Math.hypot(ux, uy) || 1 ; if (Math.hypot(ux, uy) < 1e-3 ) { ux = 1 ; uy = 0 ; } else { ux /= ul; uy /= ul; } const vx = dy * uz - dz * uy, vy = dz * ux - dx * uz, vz = dx * uy - dy * ux; const w = Math.min(0.05 , 0.0065 * f / mean), o = nr * 16 ; R2[o] = ux * w; R2[o + 1 ] = uy * w; R2[o + 2 ] = uz * w; R2[o + 3 ] = 0 ; R2[o + 4 ] = dx * len; R2[o + 5 ] = dy * len; R2[o + 6 ] = dz * len; R2[o + 7 ] = 0 ; R2[o + 8 ] = vx * w; R2[o + 9 ] = vy * w; R2[o + 10 ] = vz * w; R2[o + 11 ] = 0 ; R2[o + 12 ] = (ax + bx) / 2 ; R2[o + 13 ] = (ay + by) / 2 ; R2[o + 14 ] = (az + bz) / 2 ; R2[o + 15 ] = 1 ; nr++; } } s.rods.count = nr; s.rods.instanceMatrix.needsUpdate = true ; } const lines = (v: View3, pts: number[]) => { const g = new THREE.BufferGeometry(); g.setAttribute("position" , new THREE.Float32BufferAttribute(pts, 3 )); return new THREE.LineSegments(g, v.lineMat(1 )); };/** 長方形の枠の線 */ const rect = (x0: number, x1: number, y0: number, y1: number, z0: number, z1: number) => x0 === x1 ? [x0, y0, z0, x0, y1, z0, x0, y1, z0, x0, y1, z1, x0, y1, z1, x0, y0, z1, x0, y0, z1, x0, y0, z0] : z0 === z1 ? [x0, y0, z0, x1, y0, z0, x1, y0, z0, x1, y1, z0, x1, y1, z0, x0, y1, z0, x0, y1, z0, x0, y0, z0] : [x0, y0, z0, x1, y0, z0, x1, y0, z0, x1, y0, z1, x1, y0, z1, x0, y0, z1, x0, y0, z1, x0, y0, z0]; const circle = (R: number, y: number, n = 48 ) => { const p: number[] = []; for (let i = 0 ; i < n; i++) { const a = i / n * Math.PI * 2 , b = (i + 1 ) / n * Math.PI * 2 ; p.push(R * Math.cos(a), y, R * Math.sin(a), R * Math.cos(b), y, R * Math.sin(b)); } return p; };function plane(v: View3, x0: number, x1: number, z0: number, z1: number, y: number, tone: number) { const g = new THREE.PlaneGeometry(x1 - x0, z1 - z0); g.rotateX(-Math.PI / 2 ); const m = new THREE.Mesh(g, v.mat(tone, { double: true })); m.position.set((x0 + x1) / 2 , y, (z0 + z1) / 2 ); m.receiveShadow = true ; return m; } function clearDeco(s: S) { for (const gr of [s.deco, s.world]) for (const o of [...gr.children]) { gr.remove(o); const m = o as THREE.Mesh; m.geometry?.dispose?.(); (Array.isArray(m.material) ? m.material : [m.material]).forEach((x) => x?.dispose?.()); } for (const o of [s.dam, s.plateMesh]) if (o) { o.parent?.remove(o); o.geometry.dispose(); (o.material as THREE.Material).dispose(); } s.dam = null ; s.plateMesh = null ; } function plateRing(s: S) { if (s.plateMesh) { s.plateMesh.parent?.remove(s.plateMesh); s.plateMesh.geometry.dispose(); (s.plateMesh.material as THREE.Material).dispose(); } const sim = s.sim, g = new THREE.RingGeometry(sim.holeD / 2 , sim.siloR + 0.04 , 48 , 1 ); g.rotateX(-Math.PI / 2 ); const m = new THREE.Mesh(g, s.v.mat(0.45 , { double: true })); m.receiveShadow = true ; m.castShadow = true ; s.plateMesh = m; s.group.add(m); } /** 場面の容器と視点を作る */ function build(api: Api, s: S, key: SceneKey) { const v = s.v, sim = s.sim, N = Math.round(BASE[api.quality] * api.num("amount" ) / 100 ), sc = Math.cbrt(N / 3000 ); clearDeco(s); s.key = key; s.tNext = 0 ; s.histN = 0 ; s.hist.fill(0 ); s.grab = "" ; sim.mu = api.num("mu" ); sim.mur = api.num("mur" ); sim.setCoh(api.num("coh" )); sim.setup(key, N, api.num("hole" ), api.reduced); s.group.rotation.set(0 , 0 , 0 ); s.spout.visible = isPile(key); s.tool.visible = false ; v.orbit.auto = 0 ; if (isPile(key)) { s.deco.add(plane(v, -6 , 6 , -6 , 6 , 0 , 0.13 ), lines(v, rect(-6 , 6 , 0.002 , 0 , -6 , 6 ))); if (key === "chains" ) v.orbit.set({ yaw: 0.7 , pitch: 0.1 , dist: 7 * sc + 1 , target: [0 , 0.75 * sc, 0 ] }); else v.orbit.set({ yaw: 0.7 , pitch: 0.36 , dist: 8.2 * sc + 1.2 , target: [0 , 0.75 * sc, 0 ] }); v.shadowSize(8 ); v.sunAt([0.55 , 0.85 , 0.35 ], new THREE.Vector3(0 , 0 , 0 )); } else if (key === "slump" ) { const L0 = sim.L0, W = sim.WZ, X1 = sim.fx1, Hh = sim.damTop; s.deco.add(plane(v, -L0, X1, -W, W, 0 , 0.13 )); s.deco.add(lines(v, [...rect(-L0, X1, 0.5 , 0.5 , W, W), ...rect(-L0, X1, 0.5 , 0.5 , -W, -W), ...rect(-L0, 0 , 0 , Hh, W, W), ...rect(-L0, 0 , 0 , Hh, -W, -W), ...rect(-L0, -L0, 0 , Hh, -W, W)])); const dam = new THREE.Mesh(new THREE.BoxGeometry(0.08 , Hh, 2 * W), v.mat(0.8 )); dam.position.set(0.04 , Hh / 2 , 0 ); dam.castShadow = true ; s.dam = dam; s.group.add(dam); v.orbit.set({ yaw: 0.3 , pitch: 0.24 , dist: 12.5 * sc + 1 , target: [2.0 * sc, 1.2 * sc, 0 ] }); v.shadowSize(11 ); v.sunAt([0.3 , 0.85 , 0.45 ], new THREE.Vector3(2.5 , 0 , 0 )); } else if (key === "tilt" ) { const L = sim.LX, W = sim.WZ, Hw = 3.1 * sc; s.deco.add(plane(v, -L, L, -W, W, -0.005 , 0.18 )); s.deco.add(lines(v, [...rect(-L, L, 0 , Hw, W, W), ...rect(-L, L, 0 , Hw, -W, -W), ...rect(-L, -L, 0 , Hw, -W, W), ...rect(L, L, 0 , Hw, -W, W)])); s.world.add(plane(v, -9 , 9 , -7 , 7 , -2.6 * sc, 0.07 )); // 箱を支える台(回る軸は原点) s.world.add(lines(v, [0 , 0 , W + 0.1 , -1.1 * sc, -2.6 * sc, W + 0.1 , 0 , 0 , W + 0.1 , 1.1 * sc, -2.6 * sc, W + 0.1 , 0 , 0 , -W - 0.1 , -1.1 * sc, -2.6 * sc, -W - 0.1 , 0 , 0 , -W - 0.1 , 1.1 * sc, -2.6 * sc, -W - 0.1 , 0 , 0 , W + 0.1 , 0 , 0 , -W - 0.1 ])); v.orbit.set({ yaw: 0.22 , pitch: 0.18 , dist: 12.5 * sc + 1 , target: [0 , 0.1 * sc, 0 ] }); v.shadowSize(8 ); v.sunAt([0.3 , 0.85 , 0.45 ], new THREE.Vector3(0 , 0 , 0 )); } else { const Rs = sim.siloR, y = sim.floorY, top = N * LAT_V / (Math.PI * Rs * Rs) * 0.75 + 0.3 ; s.deco.add(plane(v, -6 , 6 , -6 , 6 , y, 0.13 )); const pts = [...circle(Rs, 0 ), ...circle(Rs, top)]; for (let k = 0 ; k < 10 ; k++) { const a = k / 10 * Math.PI * 2 ; pts.push(Rs * Math.cos(a), 0 , Rs * Math.sin(a), Rs * Math.cos(a), top, Rs * Math.sin(a)); } // 板を支える脚 for (const [a, b] of [[1 , 1 ], [-1 , 1 ], [1 , -1 ], [-1 , -1 ]]) pts.push(a * Rs * 0.7 , 0 , b * Rs * 0.7 , a * Rs * 0.7 , y, b * Rs * 0.7 ); s.deco.add(lines(v, pts)); plateRing(s); v.orbit.set({ yaw: 0.45 , pitch: 0.16 , dist: 13.5 * sc + 1 , target: [0 , 0.3 * sc, 0 ] }); v.shadowSize(8 ); v.sunAt([0.55 , 0.85 , 0.35 ], new THREE.Vector3(0 , 0 , 0 )); } } 球は三角形80個で、InstancedMeshという、同じ形をまとめて描く部品で並べる。毎コマ、位置と半径を行列の配列へ、濃さを色の配列へ写すだけで、部品は作り直さない。影と視点と光の当たり方は、three.jsと共通部品に任せている。
力の鎖は、平均の1.6倍を超える力の接点を細い円柱で描く。粒の奥に隠れないよう、奥行きを無視して最後に描き、断面より手前の粒は半径を0にして消している。
画質を下げると粒の数が減り、同じ場面を作り直す。「高」の3,000粒は、山が育つと1コマ（刻み4回）の計算が5〜8 msほどになる。「中」の1,500粒なら半分ほどで済む。砂時計だけは、粒が減っても筒を細くしすぎないよう、半径を粒の数の6乗根に合わせて変える。筒が穴の直径の2.5倍より細いと、流量が砂の高さで変わってしまうからである。
