へぇー。剛体球＋重力だけで水面っぽく見えるのか。物理屋だったら「引力なきゃ表面張力ないから駄目じゃん」で試しもしないところだろうけど、もしかしたら深い物理的意味があるのかもしれない。剛体球の斥力だけで固体液体相転移が出ると示してそれまでの常識を覆したアルダー転移みたいに。
http://homepage1.nifty.com/kaneko/water2.htm
論文誌に載るようなちゃんとした計算だと剛体球の計算って非常にテクニカルで難しいんだけど（タイムドリブンでなくイベントドリブンになる）、上のコードは見た目が水面に見えればいいからけっこういいかげんな計算で高速化してる。しかしその方向にいくなら、まだまだ高速化できるぜ。
多分一番重いのは球同士の距離を求める時の sqrt()の計算。しかし dt が常識的な値なら距離が 2R より小さくなってもそれほど深く埋まることはないはず。玉を離す処理は中心の相対位置を(2rx,2ry,2rz), 距離 2r=|2rx,2ry,2rz|, 半径を Rとすると (R-r)*(rx,ry,rz)/r だけ移動することになる。rx,ry,rz,R は既知なんで、 (R-r)/r = R/r -1 を計算すればよい。このとき r は R とあまり違わなくて、 r^2 = R^2 -δ と書けるなら、R/r は (1-δ/R^2)^(-1/2) 〜 1+(R^2-r^2)/2R^2 で近似できて sqrt() 計算は必要ない。
あるいはsqrt をテーブル引きにしてもいいかも。
まあしかし一番重いのは可視化の部分だろうから、あまり効果ないかもな。
追記：げげ！自分は彼と同じオフィスにいた時期があったらしい。話したことないから気がつかなかった。顔も覚えてないし。あんまりめったなこと書けないなぁ。
しかしＭＤの話を聞く機会もあったろうに、粒子系のシミュレーションがN^2 アルゴリズムってのはユルいと思うね。
シミュレーション領域を一辺 2R のセルで区切る。セル構造体には隣接する26個のセルへのポインタと、セルにいる粒子の代表へのポインタを持たせる。粒子の構造体にはそれが属するセルへのポインタと、同じセルに属する粒子にアクセスするための prev, next というポインタストリングのためのメンバを持たせる。衝突する可能性のある粒子は ThisParticle->MyCell->NeighborCell->FirstParticle ->next->next->next ... とすればスキャンできる。計算時間は粒子数に比例。
トリビア：メトロポリスによるモンテカルロ法の原論文は剛体球の問題を扱っているが、これは核爆弾の爆縮のための研究であったというもっぱらの噂 via T大O先生。google:implosion "metropolis algorithm" 実はこないだのオンラインカジノの乱数の話も。
昨日の自分のコメント補足。
格子ボルツマン方は位置と速度の同時分布関数という位相空間での確率分布の時間発展を数値的にシミュレーションする方法。空間は細かいメッシュに切るけど、速度は二次元なら上下左右斜め＋動かないの９方向だけにする。これでたとえばある格子で上向き速度をもってる振幅が上の格子点に移る、という感じで時間発展し、同じ格子に左右から来た成分が衝突して速度が上下に変わるとかの衝突の効果を入れてシミュレーションする。こんなんでナビエ・ストークス方程式が再現できるらしい。
