今回は「波動方程式」という方程式を使って波の表現にチャレンジしました。
ちなみに実装にはこの動画（波の伝播シミュレーションの実装解説 | ニコニコ動画）が非常に参考になりました。
波動方程式
まずはざっくり方程式。
ちなみに色々調べてみたところ、電波や弦の振動、バネなど「波」に関する公式があるようですが、そのどれも突き詰めていくとこの波動方程式という形に落ち着くらしいです。
\frac{1}{s^2} \frac{\partial^2 u}{\partial t^2} = \Delta u
波の高さの変化の加速度に定数を掛けたもの = 現在の波の高さに「ラプラシアンフィルタ」を掛けたもの
ということのようです。
ちなみに「ラプラシアンフィルタ」は以下の形です。
| 0 | 1 | 0 |
|---|---|---|
| 1 | -4 | 1 |
| 0 | 1 | 0 |
これをどう使うかは、以前書いた畳込み処理の記事を見てもらうと分かると思います。
ラプラスの演算子
\Delta = \frac{\partial^2}{\partial x_1^2} + \frac{\partial^2}{\partial x_2^2} + ... \frac{\partial^2}{\partial x_n^2}
\Delta u は、 u = u(x, y) において、
\Delta u = \frac{\partial^2}{\partial x^2}u + \frac{\partial^2}{\partial y^2}u
二階微分の計算は、
// 地点 x = n における x 方向の傾きの変化
(u[n+1] - u[n]) - (u[n] - u[n-1]) = u[n+1] + u[n-1] - 2 * u[n]
これを2次元平面に拡張し、地点(x, y)において傾きの変化（X方向の傾きの変化＋Y方向の傾きの変化）を求める計算は以下になります。
(u[x+1][y] + u[x-1][y]) + (u[x][y+1] + u[x][y-1]) - (4 * u[x][y]);
上記の計算が各波の点に「ラプラシアンフィルタ」をかけたものになります。
すごくざっくり書くと、各頂点に対して毎フレームごとにこの「ラプラシアンフィルタ」を掛けてやるとあら不思議、波の状況が再現できる、というわけです。
実際のサンプルをjsdo.itにあげているのでそれを見てもらうとだいぶ短い計算で実現しているのが分かるかと思います。
サンプルコード解説
値の初期化
これは単純に各頂点の配列を作り、0で初期化しています。
（波を作るために1点にだけ値を設定しています）
initialize
// 各頂点の高さ（波の高さ）
var U = [];
// 波の高さの速度
var Vel = [];
// 初期化
for (var x = 0; x < MAX_X; x++) {
U[x] = [];
Vel[x] = [];
for (var y = 0; y < MAX_Y; y++) {
U[x][y] = 0;
Vel[x][y] = 0;
}
}
// 中央の頂点に値を設定
U[MAX_X/2][MAX_Y/2] = 10;
アップデート処理
update
// 伝播スピード
var k = 0.5;
// 減衰率
var attenuation = 0.999;
// 波のアップデート処理
function updateWave() {
for (var x = 1; x < MAX_X - 1; x++) {
for (var y = 1; y < MAX_Y - 1; y++) {
// 現在の波の高さにラプラシアンフィルタをかけ、
// それを加速度とする。
var accel = U[x ][y-1]
+ U[x ][y+1]
+ U[x-1][y ]
+ U[x+1][y ]
- U[x ][y ] * 4;
// 伝播速度を掛ける
accel *= k;
// 現在の速度に加速度を足し、さらに減衰率を掛ける
Vel[x][y] = (Vel[x][y] + accel) * attenuation;
}
}
// 現在の波の高さに速度を足し、波の高さを更新する
for (var x = 1; x < MAX_X - 1; x++) {
for (var y = 1; y < MAX_Y - 1; y++) {
U[x][y] += Vel[x][y];
}
}
// 計算した波の高さをメッシュの頂点に適用する
for (var x = 0; x < MAX_X - 1; x++) {
for (var y = 0; y < MAX_Y - 1; y++) {
var idx = x + (MAX_X * y);
geometry.vertices[idx].z = U[x][y];
}
}
geometry.computeFaceNormals();
geometry.verticesNeedUpdate = true;
geometry.normalsNeedUpdate = true;
geometry.uvsNeedUpdate = true;
}
頂点とインデックスから平面を作成
今回の波を表現するにあたって、単純に平面を作ってしまうとそれぞれの頂点にアクセスする処理がめんどくさくなるので自前で生成、操作するようにしました。
以下が平面の頂点作成とUV座標配列の生成処理です。
var uvs = [];
var w = MAX_X - 1;
var h = MAX_Y - 1;
for (var i = 0; i < MAX_X; i++) {
for (var j = 0; j < MAX_Y; j++) {
geometry.vertices.push(new THREE.Vector3(i, j, 0));
uvs.push(new THREE.Vector2(i / w, j / h));
}
}
for (var j = 0; j < MAX_Y - 1; j++) {
for (var i = 0; i < MAX_X - 1; i++) {
var idx0 = (j * MAX_X) + i;
var idx1 = idx0 + 1;
var idx2 = idx0 + MAX_X;
var idx3 = idx1;
var idx4 = idx2;
var idx5 = idx2 + 1;
geometry.faces.push(new THREE.Face3(idx0, idx1, idx2));
geometry.faces.push(new THREE.Face3(idx3, idx4, idx5));
geometry.faceVertexUvs[0].push([
uvs[idx0],
uvs[idx1],
uvs[idx2]
]);
geometry.faceVertexUvs[0].push([
uvs[idx3],
uvs[idx4],
uvs[idx5]
]);
}
}