雪の結晶の6本の枝は、同じ空を一緒に落ちるからそろう

6本の枝は互いに連絡を取り合っていない。同じ時に同じ温度と湿り気の中を通るからそろい、通ってきた道が結晶ごとに違うから、複雑に育った結晶で同じ形はまず見つからない。

作成日 2026年11月4日更新日 2026年11月4日執筆・監修佐々木航

記事

作成日 2026年11月4日・更新日 2026年11月4日・執筆・監修佐々木航

現象

枝どうしは、連絡を取り合っていない

手袋に落ちた雪の結晶を見ると、6本の枝が同じ形をしている。向かい合う枝は数 mmも離れているのに、脇に出た小枝の位置までよく似ているのはなぜだろう。一方で「同じ形の雪は2つとない」ともいう。

答えは、結晶が通ってきた道にある。雪の結晶を研究するカリフォルニア工科大学のリブレクトは、枝がそろう理由を「枝のあいだの連絡ではなく、共通の履歴だ」と説明する¹。6本は1つの結晶として一緒に雲の中を落ちるので、同じ時に同じ温度と同じ湿り気を受ける。だから同じ形に育つ。そして、乱れた空気の中で2つの結晶がまったく同じ道をたどることはない。

画面の結晶も、六角形の小さな種から育っていく。画面の下の図は中谷ダイヤグラムで、点線の輪の中の点が今の温度と過飽和度（氷から見た湿り気）を指す。点をドラッグすると形が変わり、図の上をなぞると、その道すじのとおりに結晶が育つ。結晶のまわりをドラッグすれば、視点が回る。

仕組み

角が先に水蒸気を受け取り、枝になる

雪の結晶は、雲の中の水蒸気が氷の表面に直接付いて育つ。水の分子が六角形に並ぶ氷は、小さなうちは平らな六角形の板や柱になる¹。枝の始まりは、この六角形の角だ。

水蒸気は濃い所から薄い所へ広がる。結晶の表面で水蒸気が氷になると、そのまわりは薄くなり、遠くの空気から水蒸気が流れ込んでくる。空気中の水蒸気の濃さ $c$ （m⁻³）は拡散の式に従う。

\frac{\partial c}{\partial t} = D \nabla^{2} c

$D$ は水蒸気が空気の中を広がる速さを表す拡散係数で、1気圧でおよそ $2 \times 1 0^{- 5}$ m²/秒である²。距離の2乗を $D$ で割って見積もると、1 mmを広がるのにかかる時間はおよそ0.05秒にすぎない。何十分もかけて育つ結晶のまわりでは、濃さの分布はすぐに落ち着く。

六角形の角は、まわりの空気へ少し突き出している。濃さの分布を解くと、角のそばでは濃さの傾きが急になり、辺の真ん中より多くの水蒸気が届く。角が速く伸びるとさらに突き出して、もっと受け取る。リブレクトはこの循環を枝分かれの不安定性と呼び、6本の枝の出発点だと説明している¹。

角は辺の真ん中より多くの水蒸気を受け取り、伸びるとさらに多く受け取る

ただし、水蒸気が届くだけでは結晶にならない。分子が平らな面にうまく収まるかどうかは面の種類で違い、ここで形が決まる。表面の成長の速さ $v_{n}$ は次の形に書かれる³。

v_{n} = α v_{kin} σ_{surf}, σ = \frac{c - c _{sat}}{c _{sat}}

$σ$ は過飽和度で、氷にとっての飽和の濃さ $c_{sat}$ を超えた分の割合だ。 $σ_{surf}$ は表面のすぐ上での値、 $v_{kin}$ は水蒸気の分子が飛んでくる速さで決まる速度である。 $α$ は当たった分子が氷に収まる割合で、乱れた面では1に近く、平らな面では小さい。この $α$ が六角形の上下の面と横の面で違い、どちらが育ちやすいかは温度で入れ替わる。低い過飽和度では面が平らなまま育ち、高くなると枝が出る。

雲の中の水滴に接した空気は、水にとってちょうど飽和している。氷の飽和は水の飽和より低いので、その空気は氷から見ると過飽和で、−5 °Cで5%、−15 °Cで15.7%、−30 °Cで34%になる⁴。図の太い曲線がこの「水の飽和」である。

仕組み

上に積もるか、横へ広がるか

画面の計算は、グラブナーとグリフィースが2009年に示した3次元の格子モデルの型を使う⁵。空間を六角柱のますで埋め、ますごとに3つの量を持たせる。結晶の一部かどうか、結晶の縁に付いたまだ固まっていない水の量 $b$ 、水蒸気の量 $d$ だ。1刻みの計算は、混ぜる、凍る、付く、融けるの4段階で進む。

まず水蒸気を隣のますと混ぜる。横の隣6つに重み $a$ 、上下の隣2つに重み $1.5 a$ を掛けて足すと、どの向きにも同じ速さで広がる。

d^{'} = (1 - 9 a) d + a 横 \sum d_{隣} + 1.5 a 縦 \sum d_{隣}

結晶の縁のますでは、水蒸気の90%が $b$ に移る。 $b$ がしきい値 $β$ を超えたますが結晶になり、残った $b$ の一部は水蒸気に戻る。

b \geq β (n_{T}, n_{Z})

ここが形を決める。 $n_{T}$ は横の隣のうち結晶のますの数、 $n_{Z}$ は上下の隣のうち結晶のますの数である。隣が多い段の角ほど $β$ は小さく、平らな面の上に新しい層を始めるときほど大きい。分子が段差には収まりやすく、平らな面には収まりにくいことを、この差が表す。

横から見た格子。上に積もるますは真下の1つ、横へ広がるますは横の1つとだけ接している

とくに効くのが、上に積もるときの $β$ と横へ広がるときの $β$ の比である。グラブナーらが示した値を比べると、このことがよく分かる。

形	$β$ （上）	$β$ （横）	水蒸気 $ρ$
羊歯状	3	1.4	0.14
樹枝状	1.6	1.5	0.1
針状	2	4	0.1
中空の柱	1	2	0.1

$β$ （上）のほうが大きければ板に、 $β$ （横）のほうが大きければ柱になる。樹枝状の行では差が1割もないのに、それでも板になる。

温度と形をつなぐ部分は、計算で導いたものではない。温度ごとの $β$ の表は、この4行と論文のほかの例をもとに、中谷ダイヤグラムに合うように決めた想定の値で、過飽和度が水の飽和より低いほど大きくしてある。遠くの水蒸気の量 $ρ$ だけは、余分な水蒸気の密度 $σ \times ρ_{sat} (T)$ に比例させてある⁴。なぜ温度で板と柱が入れ替わるのかは、今も完全には分かっていない。リブレクトの2019年の論文のように、面の上を動く分子のふるまいから説明しようとする研究が続いている⁶。

試す

温度と湿り気を変えて育てる

操作パネルには、経過時間と今の条件、結晶の差し渡しと厚さ、形の名前が出る。数字は画質「高」のときの値である。画質を下げると格子が粗くなって形も少し変わり、同じ時刻の差し渡しは画質「低」で1割あまり大きくなる。

−15 °Cで水の飽和の空気では、種はすぐ六角形を離れ、11分で差し渡し1.9 mm、22分で3.6 mmの樹枝状になる。枝先は1時間に約5 mm伸びるが、これは時間の目盛りを合わせて出した値である。「見せ方」が水蒸気の濃さのときは、結晶のまわりの点がまばらになる。水蒸気が使われて薄くなった所で、枝のあいだほど点が少なく、枝先のそばにはまだ残っている。「拡散の強さ」を2倍にして育て直すと、水蒸気がすき間まで届くので枝が太くなり、0.5倍では31分でも2.7 mmにとどまる。
同じ−15 °Cのまま、過飽和度を5%（水の飽和の約3分の1）に下げる。角は先へ伸びず、差し渡し0.84 mm（32分）の六角形の板になる。水蒸気が少ないと、 $β$ を超えるまでの時間が長くなり、そのあいだに段が面の上を走って六角形を保つからだ。枝を出すには、角が辺より十分に速く伸びなければならない。
−7 °Cで水の飽和の空気では、23分で長さ2.0 mm、太さ0.54 mmの柱になる。「断面」で切り口を見ると、両端の中心がくぼんでいる。角に水蒸気が先に届き、柱の縁だけが伸びて、真ん中に穴が残るからだ。図の点を−5 °Cまでドラッグして過飽和度を5%にすると、31分で長さ2.0 mm、太さ0.41 mmの針状になる。
−7 °Cで22分育てたあと、2分かけて−15 °Cへ移す。柱の伸びはほとんど止まり、両端の縁から板が広がって、鼓の形になる。図の上を自分でなぞれば、別の履歴を試せる。なぞった時間の全体が36分に引き伸ばされるので、指を止めた所に長くとどまる。−18 °Cで止めてから−15 °Cへ動かせば、角板の先に枝がつく。
AとBは同じ履歴をたどるが、Bの温度はAの2%（−15 °Cなら0.3 °C）だけ低い。どちらかが計算の端に届く30分で比べると、結晶のますの約38%が食い違う。ゆらぎを0.5%（−15 °Cなら0.08 °C）にしても26%が違い、0にすると0%に戻る。わずかな温度の差で、枝の幅や小枝の位置が変わるのだ。
ここだけは対称性を使わず、全面を解く。ますごとに凍る割合をゆらぎの分だけ乱すと、26分の結晶で6本の枝の長さは2.34〜2.38 mmとほとんどそろうのに、60°回して重ねると約9%が食い違う。小枝の出方が枝ごとに変わるからだ。ゆらぎを1%に下げてもまだ8%ほど食い違い、0にして育て直したときだけ0%になって6本がぴったり重なる。空いた所をたたくと種を足せて、向かい合う枝は水蒸気を取り合って伸びなくなる。

現実では

兎の毛の先に、天からの手紙を読む

人工雪を初めて作ったのは北海道大学の中谷宇吉郎の研究室で、1936年3月のことだった。北海道大学の記事によると、3月12日に佐藤礒之介が羅紗の毛の上で作り、確実には再現できなかった。3月14日に助手の関戸弥太郎が、防寒服の襟から抜いた兎の毛で作ることに成功し、これが実験の中心になった⁷。中谷宇吉郎雪の科学館は、1936年3月12日を世界初の成功の日としている⁸。

中谷らは、温度と水蒸気の量を変えて結晶を作り、条件と形の関係を1枚の表にした。後に中谷ダイヤグラムと呼ばれる図である⁸。リブレクトは「−2 °C近くで板、−5 °C近くで細い柱、−15 °C近くでとくに薄い板、−30 °C近くでまた柱」とまとめている⁶。大きな星形の結晶ができるのは−15 °C近くの狭い範囲だけだという¹。形から逆に上空の条件が読めるわけで、中谷は1938年の『雪』に「雪の結晶は、天から送られた手紙であるということが出来る」と書いた⁹。

測った数字も残っている。中谷が人工雪で測った平均の成長速度は、羊歯状が1時間に4.6 mm、ふつうの樹枝と広幅の枝が1.3 mm、扇形と角板が0.66 mmだった⁹。本には「平均生長速度」とあるだけで、半径か差し渡しかは書かれていない。画面の時間の目盛りは、この4.6 mmを枝先の速さと読み、−15 °Cで水の飽和にある結晶の枝先が近い速さになるように決めた仮定である。差し渡しの速さだったなら、本物の結晶は画面の倍の時間をかけて育つことになる。

この目盛りで測ると、−15 °C、過飽和度5%の角板の縁は1時間に0.3〜0.6 mm伸びる。中谷の角板の0.66 mmに近いが、一致とは呼べない。中谷の値はこの温度と過飽和度を決めて測ったものではなく、画面の $β$ と $ρ$ も形に合わせて決めた値だからだ。言えるのは、水蒸気が少ないほど遅くなるという向きが同じことまでである。

ただ、本物の空から落ちてくる雪は、画面ほど整っていない。1931年に2,000枚を超える結晶の写真集を出したベントレーは、対称な美しいものを選んで撮っていた¹。中谷はこれを指摘し、実際には「見た眼には汚い形のものが非常に多い」と書いている。枝のそろい方も「量的の問題」で、整った結晶でも詳しく調べると枝ごとの細かな構造は対称になっていないという⁹。分類も増え続け、菊地らが2012年に提案したグローバル分類には121の小分類がある¹⁰。

もっとも「2つとない」は程度の話である。リブレクトによれば、単純な六角形の小さな結晶なら、顕微鏡で見分けのつかない2つが見つかってもおかしくない。一方、複雑な結晶では事情が変わる。本100冊の並べ方は $1 0^{158}$ 通りに近く、顕微鏡で100か所の特徴を拾える結晶でも、厳密な数え方は決まらないものの計算はこれに似てくるという¹。

同じ履歴なら本当に同じ形になるのか。リブレクトは実験室で、透明な板の上に2つの結晶を並べて同時に育て、ほとんど同じ形の「一卵性の双子」を作ってみせた。温度を大きく急に変えるのがこつで、小さくゆるやかな変化では、そろい方が悪くなるという¹¹。

省いていること

画面の結晶が省いていること

1ますは約22 μmで、元の論文が想定する1 μmほどのますよりずっと粗い⁵。そのため板は本物より厚くなり、差し渡し4.9 mmの樹枝状で厚さが約0.24 mmある。格子の端（半径2.5 mm）に届くか、描く面の数が上限に達すると計算を止める
温度と形をつなぐ表と時間の目盛りは、測定に合わせて決めた想定の値である。雲の中の風や、凝結で出る熱は計算しない。水蒸気の量が変わると、場全体がすぐに新しい濃さへ移るとした
格子の上の水蒸気は、本物と広がり方が違う。時間の目盛りを枝先の速さに合わせると、格子の拡散係数は画質「高」で本物の約8万分の1になり、代わりに水蒸気の量が氷との比で本物の数十万倍ある。画面で薄くなった所が結晶のすぐまわりに集まるのはこのためで、本物の薄い所はもっと遠くまで広がる。「拡散の強さ」は、この格子の値に掛ける倍率である
分子が面の上を動く様子や、結晶の中の乱れは、 $b$ と $β$ にまとめた。2つの結晶がくっつく雪片や、複数の結晶が集まる形は作れない
「ゆらぎ」の場面のほかは全体の24分の1だけを解いて写しており、6本がそろうのは計算の作りでもある。ゆらぎの場面も上下は鏡で写している
落ちていく向きや回転は入れていない。実際の板状の結晶は水平に近い姿勢で落ちるとされ、グラブナーらは、それが上と下で模様が違う一因だとみている⁵

参考文献

Libbrecht, K. G.「Snow Crystals」arXiv:1910.06389, 2019. 書籍版は『Snow Crystals: A Case Study in Spontaneous Structure Formation』Princeton University Press, 2021. https://arxiv.org/abs/1910.06389
Libbrecht, K. G.「A Critical Look at Ice Crystal Growth Data」arXiv:cond-mat/0411662, 2004. https://arxiv.org/abs/cond-mat/0411662
Libbrecht, K. G.「A Quantitative Physical Model of the Snow Crystal Morphology Diagram」arXiv:1910.09067, 2019. https://arxiv.org/abs/1910.09067
Murphy, D. M., Koop, T.「Review of the vapour pressures of ice and supercooled water for atmospheric applications」『Quarterly Journal of the Royal Meteorological Society』131(608), 1539–1565, 2005.
Gravner, J., Griffeath, D.「Modeling snow-crystal growth: A three-dimensional mesoscopic approach」『Physical Review E』79(1), 011601, 2009. 著者の公開版（題はModeling snow crystal growth III、2007年の版） https://www.math.ucdavis.edu/~gravner/papers/h3l.pdf
中谷宇吉郎『雪』岩波新書, 岩波書店, 1938. 青空文庫 https://www.aozora.gr.jp/cards/001569/card52468.html
北海道大学 CoSTEP「#3 天からの手紙、中谷からの手紙」. https://costep.open-ed.hokudai.ac.jp/like_hokudai/article/6708
中谷宇吉郎雪の科学館「中谷宇吉郎について」. https://yukinokagakukan.kagashi-ss.com/ukichiro/

参考文献1の1章。枝がそろう理由は「not by any communication between them, but rather by their common history」。枝分かれの不安定性、−15 °Cの大きな星形、写真集の枚数とベントレーの選び方も同じ章による。 ↩ ↩² ↩³ ↩⁴ ↩⁵ ↩⁶
参考文献2。 ↩
参考文献3の式(1)。 ↩
参考文献4の水と氷の飽和水蒸気圧の式から計算した。 $ρ_{sat}$ は氷の飽和水蒸気密度で、−15 °Cで1.39 g/m³。 ↩ ↩²
参考文献5。数値と章の番号は著者の公開版による。表の値は8章（樹枝状）、10章（羊歯状）、11章（針状と中空の柱）、ますの大きさは7章、落ちる向きの話は10章。式の重みは、論文が2段階に分けて混ぜているのを1段階にまとめ、等方に広がる比は同じにした。論文が先行の試みに挙げるのは、Packard, N. H.「Lattice models for solidification and aggregation」（1984年のプレプリント）とReiter, C. A.「A local cellular model for snow crystal growth」『Chaos, Solitons & Fractals』23(4), 1111–1119, 2005である。 ↩ ↩² ↩³
参考文献3。 ↩ ↩²
参考文献7。 ↩
参考文献8。 ↩ ↩²
参考文献6。手紙の一文と成長速度は第四の章「雪を作る話」、「見た眼には汚い形」は第二の章「「雪の結晶」雑話」、「量的の問題」は第三の章「北海道における雪の研究」による。 ↩ ↩² ↩³
菊地勝弘, 亀田貴雄, 樋口敬二, 山下晃「中緯度と極域での観測に基づいた新しい雪結晶の分類－グローバル分類－」『雪氷』74(3), 223–241, 2012. 元のMagono and Leeの分類は81小分類。 ↩
参考文献1の9.4節。 ↩

プログラム

このシミュレーションのプログラム

このシミュレーションは、1刻みごとに「水蒸気を隣のますと混ぜる → 結晶の縁で凍らせる → しきい値を超えたますを結晶にする → 少し融かして戻す」をくり返している。立体の描画と影と視点はthree.jsが受け持ち、このファイルが計算するのは、どのますが結晶になったかである。描くときは、表面に出ている面だけを三角形にしてthree.jsに渡す。

長さはmm、時間は分

計算する範囲は半径2.5 mmで、格子の半径（RQ）を画質で変える。1ますの幅は $2.5 \div R$ mmで、画質「高」では約22 μmになる。

定数と単位(src/sims/snow-crystal.ts の 31〜49 行)

const SQ = Math.sqrt(3) / 2; // 隣の列までの距離(ますの幅を 1 として)const SIZE = 2.5; // 計算する範囲の半径 mmconst RQ = [56, 84, 116], ZQ = [46, 70, 96]; // ひとつの結晶の格子: 横の半径と縦の半分の高さ(ます。画質 0・1・2)const RT = [40, 58, 76], ZT = [34, 48, 64]; // 2 つ並べるときの格子const RF = [36, 48, 64], ZF = [14, 18, 24]; // 横の全面を解くときの格子const RP = 32, ZP = 28; // 図の点を動かしている間の下書きの格子const PLAN_RS = 0.72; // 履歴をたどるときの格子の倍率(柱の両端に板が 2 枚広がると面の数が倍になるので、格子を少し粗くする)const SEC_120 = 0.27; // 半径 120 ますの格子の 1 刻み(秒)。ますが粗くなると枝先がほぼ同じ速さになるよう、その分だけ長くするconst STEP_RATE = [240, 340, 480]; // 「早送り 1 倍」で画面の 1 秒に進める刻みの数(画質 0・1・2)const BUDGET_MS = 6; // 1 コマで格子の計算に使ってよい時間const TRI_MAX = 300000; // 1 コマに描く三角形の上限。結晶がここまで細かくなったら、育てるのを止めるconst STOP_MIN = 90; // 条件を変えないとき、ここまで育てる(分)const PLAN_MIN = 36; // 図の上をなぞったときの履歴の長さ(分)const PREVIEW_MIN = 35; // 下書きを育てる時間(分)const TMIN = -30, TMAX = -1, SMIN = 1, SMAX = 40; // 温度 °C と過飽和度 % の範囲const EXPOSURE = 2.05; // 明るさの倍率(面の明るさを網点の濃さへ)const TONE0 = 0.66; // いちばん低い上の面の濃さconst FIELD = 0.115; // 遠くの水蒸気の濃さを描くときの明るさconst A0 = 0.1; // 拡散の重み(横の隣 1 つあたり)。縦の隣は 1.5 倍にすると、横と縦で同じ速さに広がる

1刻みの時間はSEC_120から決める。枝先の速さは格子の粗さによらず1刻みあたり約0.02ます（画質「高」で0.017、「低」で0.020）なので、ますが粗いほど1刻みを長くすると、同じ条件の枝先がほぼ同じ速さになる。粗い格子のほうが1割あまり速い。SEC_120は、中谷の測った羊歯状の結晶の速さ4.6 mm/時を枝先の速さと読み、それに近くなるように決めた想定である。BUDGET_MSは1コマで計算に使ってよい時間で、これを超えた刻みはそのコマでは進めない。遅い端末では結晶がゆっくり育つだけで、形は変わらない。

温度と過飽和度を係数にする

coefは、温度 $T$ と過飽和度 $σ$ から、しきい値 $β$ と融ける割合 $μ$ と遠くの水蒸気の量 $ρ$ を決める。

温度と過飽和度から係数へ(src/sims/snow-crystal.ts の 51〜91 行)

/** 氷と水の飽和水蒸気圧(Pa)。Murphy と Koop の式(2005 年) */const eIce = (T: number) => { const K = T + 273.15; return Math.exp(9.550426 - 5723.265 / K + 3.53068 * Math.log(K) - 0.00728332 * K); };const eLiq = (T: number) => { const K = T + 273.15; return Math.exp(54.842763 - 6763.22 / K - 4.21 * Math.log(K) + 0.000367 * K + Math.tanh(0.0415 * (K - 218.8)) * (53.878 - 1331.22 / K - 9.44523 * Math.log(K) + 0.014025 * K)); };/** 水滴のある雲の中の過飽和度(氷に対して・%)。この値を中谷ダイヤグラムに「水の飽和」の線で引く */const sigmaWater = (T: number) => 100 * (eLiq(T) / eIce(T) - 1);/** 氷の飽和水蒸気密度(g/m³) */const rhoSat = (T: number) => (eIce(T) / (461.5 * (T + 273.15))) * 1000; /** 温度ごとの係数の表(想定の値)。β は 1 + 値 × g(σ) で決める。行の順: T °C・β(上)・β(横)・β(横 1 縦 1)・β(横 2 縦 1)・μ(横)・μ(ほか) */const TABLE = [  [-1, 1.5, 0.8, 0.8, 0, 0.002, 0.001],  [-3, 1.0, 1.0, 1.0, 0, 0.002, 0.001],  [-5, 0.4, 3.0, 3.0, 0, 0.001, 0.001],  [-7, 0.0, 1.0, -0.5, -0.5, 0.01, 0.01],  [-10, 0.5, 0.6, 0, 0, 0.015, 0.015],  [-12, 1.5, 1.0, 1.0, 0, 0.001, 0.001],  [-15, 2.0, 0.4, 0.4, 0, 0.005, 0.001],  [-18, 1.8, 1.2, 0.6, 0.6, 0.001, 0.0005],  [-22, 1.5, 1.0, 1.0, 0, 0.001, 0.001],  [-26, 0.5, 0.6, 0, 0, 0.015, 0.015],  [-30, 0.7, 1.0, 0, 0, 0.01, 0.01],];interface Coef { beta: Float32Array; mu: Float32Array; kappa: number; rho: number; A: number; B: number; C0: number; passes: number }const newCoef = (): Coef => ({ beta: new Float32Array(8), mu: new Float32Array(8), kappa: 0.1, rho: 0.1, A: A0, B: 1.5 * A0, C0: 1 - 9 * A0, passes: 1 }); /** T(°C)・σ(%)・拡散の強さ D から係数を決める。表の添え字は nT(3 まで)+ 4 ×(nZ が 1 以上か) */function coef(C: Coef, T: number, sig: number, D: number) {  const t = clamp(T, TMIN, TMAX);  let i = 0;  while (i < TABLE.length - 2 && t < TABLE[i + 1][0]) i++;  const f = (TABLE[i][0] - t) / (TABLE[i][0] - TABLE[i + 1][0]), at = (n: number) => TABLE[i][n] + (TABLE[i + 1][n] - TABLE[i][n]) * f;  // σ が水の飽和より低いほど、結晶の角から新しい層を始めにくくする(g が大きいほど β が高い)  const x = sig / Math.max(sigmaWater(t), 3), g = Math.min(3, Math.sqrt(1 / Math.max(x, 0.05)));  const up = at(1), side = at(2), b11 = at(3), b21 = at(4), mt = at(5), mo = at(6), B = C.beta, M = C.mu, th = (v: number) => Math.max(0.3, 1 + v * g);  B[1] = B[2] = th(side); B[3] = th(Math.min(0, b21)); B[4] = th(up); B[5] = th(b11); B[6] = th(b21); B[7] = 1;  M[1] = M[2] = mt; M[3] = M[4] = M[5] = M[6] = M[7] = mo;  // 遠くの水蒸気の量 ρ: 余分な水蒸気の密度(g/m³)に比例。0.64 は −15 °C の水の飽和が手本の ρ = 0.14 になる倍率  C.rho = clamp(0.64 * (sig / 100) * rhoSat(t), 0.004, 0.22);  // 拡散の強さ: 1 回の混ぜる量を変える。1 より強いときは 1 刻みに 2 回混ぜる  C.passes = D > 1 ? 2 : 1; C.A = (A0 * D) / C.passes; C.B = 1.5 * C.A; C.C0 = 1 - 9 * C.A;}

TABLEが温度ごとの表で、値はグラブナーとグリフィースの論文の例をもとに、中谷ダイヤグラムに合うように決めた。表のあいだは直線でつなぐ。 $β$ は $1 + 値 \times g$ の形で、 $g$ は水の飽和の過飽和度 $σ_{w}$ と今の $σ$ から $g = σ_{w} / σ$ （3まで）と決めた想定である。 $σ$ が低いと新しい層を始めにくくなり、面が平らなまま育つ。eIceとeLiqはマーフィーとコープの飽和水蒸気圧の式で、sigmaWaterが図の「水の飽和」の曲線を与える。 $ρ$ は余分な水蒸気の密度に比例させた。比例の係数0.64は、−15 °Cの水の飽和が論文の羊歯状の $ρ = 0.14$ になるように決めた想定である。

全体の24分の1だけを解く

結晶が上下と6回の対称を保つ場面では、全体の24分の1を解けば足りる。makeLatは柱（横の1か所に縦に並ぶます）を中心からの距離の順に並べ、隣の表を作る。

格子を作る(src/sims/snow-crystal.ts の 93〜168 行)

const DI = [1, 0, -1, 0, 1, -1], DJ = [0, 1, 0, -1, -1, 1]; // 横の隣 6 つ(軸の座標 i・j)const hexDist = (i: number, j: number) => Math.max(Math.abs(i), Math.abs(j), Math.abs(i + j));interface Lat {  R: number; Nz: number; S: number; wedge: boolean; nc: number; ci: Int16Array; cj: Int16Array; sh: Int16Array; wt: Uint8Array; nb: Int32Array; shellEnd: Int32Array; index: Int32Array;  n: number; d: Float32Array; dn: Float32Array; b: Float32Array; a: Uint8Array; nt: Uint8Array; nz: Uint8Array; isb: Uint8Array; flat: Uint8Array; top: Int16Array;  bl: Int32Array; nbl: number; srf: Int32Array; nsrf: number; att: Int32Array; rho: number; rmax: number; zmax: number; t: number; mx: number; mz: number; area: number; cells: number; dirty: boolean;}/** 回転と鏡で、全体の 1/12 の扇(0 ≦ j ≦ i)へ写す */function canon(i: number, j: number): [number, number] {  for (let r = 0; r < 6; r++) { if (i >= 0 && j >= 0) return j > i ? [j, i] : [i, j]; const t = i; i = -j; j = t + j; }  return [0, 0];}function colAt(L: Lat, i: number, j: number) {  if (L.wedge) [i, j] = canon(i, j);  if (hexDist(i, j) > L.R) return L.nc;  return L.index[(i + L.R) * (2 * L.R + 1) + j + L.R];}/** 柱(横の 1 か所に縦に並ぶます)を、中心からの距離の順に並べる。隣の表には、扇の外へ出るところに鏡で写した先を入れる */function makeLat(R: number, Nz: number, wedge: boolean): Lat {  const ci: number[] = [], cj: number[] = [], sh: number[] = [];  if (wedge) { for (let s = 0; s <= R; s++) for (let j = 0; j <= s >> 1; j++) { ci.push(s - j); cj.push(j); sh.push(s); } }  else { const a: number[][] = []; for (let i = -R; i <= R; i++) for (let j = -R; j <= R; j++) { const h = hexDist(i, j); if (h <= R) a.push([h, i, j]); } a.sort((p, q) => p[0] - q[0] || p[1] - q[1] || p[2] - q[2]); for (const [h, i, j] of a) { ci.push(i); cj.push(j); sh.push(h); } }  const nc = ci.length, S = Nz + 1, n = (nc + 1) * S, index = new Int32Array((2 * R + 1) ** 2).fill(nc);  const wt = new Uint8Array(nc), shellEnd = new Int32Array(R + 2);  for (let c = 0; c < nc; c++) { index[(ci[c] + R) * (2 * R + 1) + cj[c] + R] = c; wt[c] = !wedge ? 1 : ci[c] === 0 && cj[c] === 0 ? 1 : cj[c] === 0 || cj[c] === ci[c] ? 6 : 12; shellEnd[sh[c]] = c + 1; }  const L: Lat = {    R, Nz, S, wedge, nc, ci: Int16Array.from(ci), cj: Int16Array.from(cj), sh: Int16Array.from(sh), wt, nb: new Int32Array(nc * 6), shellEnd, index,    n, d: new Float32Array(n), dn: new Float32Array(n), b: new Float32Array(n), a: new Uint8Array(n), nt: new Uint8Array(n), nz: new Uint8Array(n), isb: new Uint8Array(n), flat: new Uint8Array(n), top: new Int16Array(nc + 1),    bl: new Int32Array(n >> 1), nbl: 0, srf: new Int32Array(n), nsrf: 0, att: new Int32Array(1 << 15), rho: 0.1, rmax: 0, zmax: 0, t: 0, mx: Math.round(R / 4), mz: Math.round(R / 5), area: 0, cells: 0, dirty: true,  };  for (let c = 0; c < nc; c++) for (let d = 0; d < 6; d++) L.nb[c * 6 + d] = colAt(L, ci[c] + DI[d], cj[c] + DJ[d]);  return L;}/** ます k の、結晶の隣の数を数え直す。鏡で写した隣は 2 回数えるので、全体の格子で数えるのと同じになる */function recount(L: Lat, k: number) {  const { S, nb, a } = L, w = (k / S) | 0, z = k - w * S;  let t = 0;  for (let e = w * 6, f = e + 6; e < f; e++) t += a[nb[e] * S + z];  L.nt[k] = t; L.nz[k] = z === 0 ? 2 * a[k + 1] : a[k - 1] + a[k + 1];}function touch(L: Lat, q: number) {  if (L.a[q]) return;  if (!L.isb[q]) { L.isb[q] = 1; L.bl[L.nbl++] = q; }  recount(L, q);}/** ます k を結晶にする。まわりのますが結晶の外の縁になる */function attach(L: Lat, k: number) {  const { S, nb, nc, sh, top } = L, w = (k / S) | 0, z = k - w * S;  L.a[k] = 1; L.isb[k] = 0; L.srf[L.nsrf++] = k; L.d[k] = 0; L.dn[k] = 0; L.dirty = true;  if (sh[w] > L.rmax) L.rmax = sh[w];  if (z > L.zmax) L.zmax = z;  if (top[w] < 0) L.area += L.wt[w];  if (z > top[w]) top[w] = z;  L.cells += L.wt[w] * (z === 0 ? 1 : 2);  for (let e = w * 6, f = e + 6; e < f; e++) if (nb[e] !== nc) touch(L, nb[e] * S + z);  if (z + 1 < L.Nz) touch(L, k + 1);  if (z > 0) touch(L, k - 1);}/** 種: 半径 r の薄い六角形。初めから水蒸気の量 ρ で満たす */function clearLat(L: Lat, rho: number) {  L.d.fill(rho); L.dn.fill(rho); L.b.fill(0); L.a.fill(0); L.nt.fill(0); L.nz.fill(0); L.isb.fill(0); L.flat.fill(0); L.top.fill(-1);  L.nbl = 0; L.nsrf = 0; L.rho = rho; L.rmax = 0; L.zmax = 0; L.t = 0; L.area = 0; L.cells = 0; L.dirty = true;}function addSeed(L: Lat, ci0: number, cj0: number, r = 2) {  for (let c = 0; c < L.nc; c++) if (hexDist(L.ci[c] - ci0, L.cj[c] - cj0) <= r) { const k = c * L.S; if (!L.a[k]) { L.b[k] = 1; attach(L, k); } }  // 種を順に置く途中で、次に種になるますも縁の表に入る。結晶になったますは表から外す(残すと、種から水蒸気が漏れ、置いた順番で枝がずれる)  let w = 0;  for (let i = 0; i < L.nbl; i++) { const k = L.bl[i]; if (!L.a[k]) L.bl[w++] = k; }  L.nbl = w;}/** 遠くの水蒸気の量が変わったら、場全体を同じ倍率にする(本物の空気では、拡散は成長よりずっと速く、場はすぐに新しい形に落ち着く) */function setRho(L: Lat, rho: number) {  const f = rho / L.rho, { d, dn, n } = L;  for (let i = 0; i < n; i++) { d[i] *= f; dn[i] *= f; }  L.rho = rho;}

扇の外へ出る隣は、canonで鏡に写した先のますを指す。高さ0の面も鏡で、下の隣は上の隣と同じになる。こうすると計算量は24分の1になり、半径116ますの格子でも、1刻みはおよそ1 ms以下で済む。「ゆらぎ」の場面ではwedgeをfalseにして、横の全面を解く。ゆらぎが0なら、全面を解いても24分の1を写したのと同じ結晶になる。recountは隣の結晶のますを数えるところで、鏡で写した隣は2回数えるので、全体を解いたときと同じ数になる。

混ぜて、凍らせて、付ける

本文の4段階をそのまま書いた部分である。

格子の1刻み(src/sims/snow-crystal.ts の 170〜209 行)

/** 拡散: 横 6 つと縦 2 つの隣と混ぜる。結晶のますは d = 0 として足し、結晶の隣のますは跳ね返った分をあとで足す */function diffuse(L: Lat, C: Coef) {  const { d, dn, nb, S } = L, A = C.A, B = C.B, C0 = C.C0;  // 結晶のまわりだけを解く。外側のますは遠くの値 ρ のまま(結晶から少し離れた所は、ほぼ一様とみなす)  const wa = L.shellEnd[Math.min(L.R, L.rmax + L.mx)], za = Math.min(L.Nz, L.zmax + L.mz);  for (let w = 0; w < wa; w++) {    const base = w * S, e = w * 6, n0 = nb[e] * S, n1 = nb[e + 1] * S, n2 = nb[e + 2] * S, n3 = nb[e + 3] * S, n4 = nb[e + 4] * S, n5 = nb[e + 5] * S;    dn[base] = C0 * d[base] + A * (d[n0] + d[n1] + d[n2] + d[n3] + d[n4] + d[n5]) + 2 * B * d[base + 1]; // 高さ 0 の面は上下の鏡    for (let z = 1; z < za; z++) {      const k = base + z;      dn[k] = C0 * d[k] + A * (d[n0 + z] + d[n1 + z] + d[n2 + z] + d[n3 + z] + d[n4 + z] + d[n5 + z]) + B * (d[k - 1] + d[k + 1]);    }  }  const { srf, bl, nt, nz } = L;  for (let i = 0, m = L.nsrf; i < m; i++) dn[srf[i]] = 0;  for (let i = 0, m = L.nbl; i < m; i++) { const k = bl[i]; dn[k] += d[k] * (A * nt[k] + B * nz[k]); }  L.d = dn; L.dn = d;}/** 凍る・付く・融ける。結晶の縁のますだけを調べる。付くますはいったん集めてから、まとめて結晶にする(どの枝も同じ順序になる) */function settle(L: Lat, C: Coef, noise: number, rng: (() => number) | null) {  const { d, b, bl, nt, nz, att } = L, beta = C.beta, mu = C.mu, keep = C.kappa;  let na = 0, w = 0;  for (let i = 0, m = L.nbl; i < m; i++) {    const k = bl[i], nT = nt[k], nZ = nz[k], t = (nT > 3 ? 3 : nT) + (nZ > 0 ? 4 : 0);    let f = 1 - keep;    if (rng) f *= 1 + noise * (2 * rng() - 1); // ゆらぎ: 凍る割合を、ますごと刻みごとに少しだけ変える    const v = d[k], fr = f * v, bb = b[k] + fr;    // 横に 4 つ以上・上下に 1 つ以上の隣があるますは、穴を埋めるために β によらず結晶にする(論文と同じ)    if ((bb >= beta[t] || (nT >= 4 && nZ >= 1)) && na < att.length) { b[k] = bb + v - fr; d[k] = 0; att[na++] = k; continue; }    const me = mu[t] * bb;    b[k] = bb - me; d[k] = v - fr + me; bl[w++] = k;  }  L.nbl = w;  for (let i = 0; i < na; i++) attach(L, att[i]);  L.t++;}function stepLat(L: Lat, C: Coef, noise: number, rng: (() => number) | null) {  for (let p = 0; p < C.passes; p++) diffuse(L, C);  settle(L, C, noise, rng);}

diffuseは本文の $d^{'}$ の式で、横6つと縦2つの隣を足す。論文は横と縦を2段階に分けて混ぜるが、ここでは1段階にまとめ、横の重み $a$ と縦の重み $1.5 a$ の比で等方に広げている。 $a = 0.1$ は、中心の重み $1 - 9 a$ が負にならない範囲で大きい値だ。

この混ぜ方の拡散係数は、ますの幅 $Δ x$ と1刻みの時間 $Δ t$ から $1.5 a Δ x^{2} /Δ t$ と決まる。画質「高」では約 $2.5 \times 1 0^{- 10}$ m²/秒で、本物の $2 \times 1 0^{- 5}$ m²/秒の約8万分の1だ。水蒸気の量 $ρ = 0.14$ も、氷との比で本物の約60万倍ある（本物は−15 °Cの水の飽和で、余分な水蒸気0.22 g/m³と氷0.92 g/cm³の比の約 $2.4 \times 1 0^{- 7}$ ）。論文のモデルの時間は物理の単位を持たないので、時間の目盛りは枝先の速さだけで決めた。

結晶のますは $d = 0$ として足し、あとで結晶の隣のますに、跳ね返った分A * nt + B * nzを足す。こうすると、ループの中で結晶かどうかを調べずに済む。

settleは結晶の縁のますだけを調べる。bb >= beta[t]が本文のしきい値の式で、横に4つ以上、上下に1つ以上の隣があるますは、穴を埋めるために $β$ によらず結晶にする（論文と同じ規則）。付くますはいったんattに集めてから、まとめて結晶にする。調べる順番で結果が変わらないようにするためで、順番の影響があると、対称な条件でも枝がずれてしまう。外側のますはwaとzaの範囲の外にあり、計算せず遠くの値 $ρ$ のままにしてある。結晶が小さいうちは、これで計算がさらに軽くなる。

表面の面だけをthree.jsに渡す

ここから先は描画である。計算の結果のうち、three.jsに渡すのは三角形の座標と向きと濃さだけで、照明と影、視点は共通部品とthree.jsが受け持つ。白黒の網点への変換は全記事に共通の部品が行う。

形を描く(src/sims/snow-crystal.ts の 246〜430 行)

// three.js が受け持つのは、ここで作った三角形を照らし、影を落として画面に描くところ。計算の結果は「面がどこに出ているか」だけを渡す/** 1 つの結晶の描画の部品を作る。全体の 1/24 の形を 1 つの geometry に入れ、回転と鏡で並べて描く */function makeCry(s: S, R: number, Nz: number, wedge: boolean, offX: number): Cry {  const L = makeLat(R, Nz, wedge), dx = SIZE / R, capTri = (wedge ? 3 : 16) * R * R;  const pos = new Float32Array(capTri * 9), nor = new Float32Array(capTri * 9), tone = new Float32Array(capTri * 9), idx = new Uint32Array(capTri * 3);  const geo = new THREE.BufferGeometry();  geo.setAttribute("position", new THREE.BufferAttribute(pos, 3).setUsage(THREE.DynamicDrawUsage));  geo.setAttribute("normal", new THREE.BufferAttribute(nor, 3).setUsage(THREE.DynamicDrawUsage));  geo.setAttribute("color", new THREE.BufferAttribute(tone, 3).setUsage(THREE.DynamicDrawUsage));  geo.setIndex(new THREE.BufferAttribute(idx, 1).setUsage(THREE.DynamicDrawUsage));  const count = wedge ? 12 : 1;  const mA = new THREE.InstancedMesh(geo, s.mats.ice, count), mB = new THREE.InstancedMesh(geo, s.mats.ice, count);  mA.frustumCulled = mB.frustumCulled = false; mA.castShadow = mB.castShadow = true; mA.receiveShadow = mB.receiveShadow = true;  const flip = new THREE.Group(); flip.scale.y = -1; flip.add(mB); // 上下を返した側。three.js が面の表裏を自動で入れ替える  const group = new THREE.Group(); group.position.x = offX; group.add(mA, flip);  // 断面の切り口  const capPos = new Float32Array((wedge ? 2 : 1) * 2 * (2 * R + 2) * 4 * 18), capGeo = new THREE.BufferGeometry();  capGeo.setAttribute("position", new THREE.BufferAttribute(capPos, 3).setUsage(THREE.DynamicDrawUsage));  const cap = new THREE.Mesh(capGeo, s.mats.cap); cap.frustumCulled = false; cap.visible = false; group.add(cap);  // 水蒸気の濃さ: 高さ 0 の面。粗い三角形の頂点に濃さを持たせる  const st = R >= 100 ? 3 : 2, verts: number[] = [], cols: number[] = [], tris: number[] = [], ids = new Map<number, number>();  const vid = (i: number, j: number) => {    const key = (i + 4096) * 8192 + j + 4096;    let v = ids.get(key);    if (v === undefined) { v = cols.length; ids.set(key, v); cols.push(colAt(L, i, j)); verts.push((i + j * 0.5) * dx, 0, j * SQ * dx); }    return v;  };  const inside = (i: number, j: number) => (wedge ? i >= 0 && j >= 0 : true) && hexDist(i, j) <= R;  const lo = wedge ? 0 : -R;  for (let i = lo; i <= R; i += st) for (let j = lo; j <= R; j += st) {    if (!inside(i, j)) continue;    if (inside(i + st, j) && inside(i, j + st)) tris.push(vid(i, j), vid(i, j + st), vid(i + st, j));    if (inside(i + st, j) && inside(i, j + st) && inside(i + st, j + st)) tris.push(vid(i + st, j), vid(i, j + st), vid(i + st, j + st));  }  const fGeo = new THREE.BufferGeometry(), fCol = new Float32Array(cols.length * 3);  fGeo.setAttribute("position", new THREE.BufferAttribute(new Float32Array(verts), 3));  fGeo.setAttribute("color", new THREE.BufferAttribute(fCol, 3).setUsage(THREE.DynamicDrawUsage));  fGeo.setIndex(tris);  const field = new THREE.InstancedMesh(fGeo, s.mats.field, wedge ? 6 : 1); field.frustumCulled = false; group.add(field);  const m = new THREE.Matrix4();  for (let k = 0; k < (wedge ? 6 : 1); k++) field.setMatrixAt(k, m.makeRotationY((k * Math.PI) / 3));  s.v.scene.add(group);  const c: Cry = { L, dx, h: stepMin(R), C: newCoef(), T: s.T, s: s.sig, done: false, group, geo, pos, nor, tone, idx, mA, mB, flip, capGeo, capPos, cap, field, fGeo, fCol, fIdx: Int32Array.from(cols), tris: 0, lastR: 0, lastM: 0, rate: 0, rot: null, phase: 0, bn: 0, bm: 0, nv: 0, ni: 0 };  // 全面を解くときは、60° 回した先のますを控えておく(6 本の枝を比べるため)  if (!wedge) { const rot = new Int32Array(L.nc); for (let w = 0; w < L.nc; w++) rot[w] = colAt(L, -L.cj[w], L.ci[w] + L.cj[w]); c.rot = rot; }  setInstances(c, false);  return c;}function dropCry(s: S, c: Cry) { s.v.scene.remove(c.group); c.mA.dispose(); c.mB.dispose(); c.field.dispose(); c.geo.dispose(); c.capGeo.dispose(); c.fGeo.dispose(); } /** 並べ方の行列。断面を見るときは、手前(z が正)の半分を外す */function setInstances(c: Cry, section: boolean) {  if (!c.L.wedge) { c.mA.setMatrixAt(0, new THREE.Matrix4()); c.mB.setMatrixAt(0, new THREE.Matrix4()); c.mA.count = c.mB.count = 1; c.mA.instanceMatrix.needsUpdate = c.mB.instanceMatrix.needsUpdate = true; return; }  const all: THREE.Matrix4[] = [], mirror = new THREE.Matrix4().makeScale(1, -1, -1);  for (let k = 0; k < 6; k++) { const r = new THREE.Matrix4().makeRotationY((k * Math.PI) / 3); all.push(r.clone(), r.clone().multiply(mirror)); }  const P = new THREE.Vector3(), keep = (m: THREE.Matrix4) => P.set(Math.cos(0.26), 0, Math.sin(0.26)).applyMatrix4(m).z < 0;  const sorted = section ? [...all.filter(keep), ...all.filter((m) => !keep(m))] : all;  sorted.forEach((m, n) => { c.mA.setMatrixAt(n, m); c.mB.setMatrixAt(n, m); });  c.mA.count = c.mB.count = section ? sorted.filter(keep).length : 12;  c.mA.instanceMatrix.needsUpdate = c.mB.instanceMatrix.needsUpdate = true;} /** 結晶の表面に出ている面を三角形にする。上の面が平らに続く所は大きな三角形でまとめ、縦に続く横の面は 1 枚にする。 *  ms ミリ秒で区切り、数コマに分けて進める。書き終えたら true を返し、three.js に渡す */function buildMesh(c: Cry, section: boolean, ms: number) {  const L = c.L, { a, nb, S, Nz, ci, cj, nc, flat, srf, wedge } = L, dx = c.dx, P = c.pos, N = c.nor, G = c.tone, I = c.idx;  const maxV = P.length / 3 - 16, maxI = I.length - 48, t0 = performance.now();  if (c.phase === 0) { c.phase = 1; c.bn = 0; c.bm = 0; L.dirty = false; }  if (c.phase === 1) {    // 全部のますが結晶に囲まれたら、もう描かない。上の面が平らかどうかを調べる    let n = c.bn, m = c.bm;    for (; n < L.nsrf; n++) {      if ((n & 255) === 255 && performance.now() - t0 > ms) { c.bn = n; c.bm = m; return false; }      const k = srf[n], w = (k / S) | 0, z = k - w * S, e = w * 6, open0 = !a[k + 1];      // flat の印: 1 = 上の面が隣と同じ高さで続く・2 = 下の面が続く      let f = (open0 ? 1 : 0) | (z > 0 && !a[k - 1] ? 2 : 0), covered = !open0 && (z === 0 || !!a[k - 1]);      for (let d = 0; d < 6; d++) { const q = nb[e + d] * S + z; if (!a[q]) { covered = false; f = 0; } else { if (a[q + 1]) f &= 2; if (z === 0 || a[q - 1]) f &= 1; } }      flat[k] = f;      if (!covered) srf[m++] = k;    }    L.nsrf = m; c.phase = 2; c.bn = 0; c.nv = 0; c.ni = 0;  }  let nv = c.nv, ni = c.ni, g = 1; // g は面の濃さ(高い所ほど明るくして、段の重なりが読めるように)  const vert = (x: number, y: number, z: number, nx: number, ny: number, nz: number) => { const o = nv * 3; P[o] = x; P[o + 1] = y; P[o + 2] = z; N[o] = nx; N[o + 1] = ny; N[o + 2] = nz; G[o] = G[o + 1] = G[o + 2] = g; return nv++; };  const gz = 1 / (L.zmax + 1);  /** 平らに続く面: このますと、右と奥の隣を結ぶ三角形を書く(3 つとも平らなときだけ。扇の外へ写った隣は使わない)。周りが全部平らなら六角形は要らない */  const flatPart = (e: number, i: number, j: number, z: number, X: number, Z: number, bit: number, y: number, ny: number) => {    let deep = !(wedge && i - j <= 1);    for (let d = 0; d < 6 && deep; d++) if (!(flat[nb[e + d] * S + z] & bit)) deep = false;    const direct = !wedge || j + 1 <= i, w0 = nb[e], w1 = nb[e + 1];    if (direct && w0 !== nc && w1 !== nc && flat[w0 * S + z] & bit && flat[w1 * S + z] & bit) {      const v0 = vert(X, y, Z, 0, ny, 0), v1 = vert(X + dx, y, Z, 0, ny, 0), v2 = vert(X + dx * 0.5, y, Z + SQ * dx, 0, ny, 0);      I[ni++] = v0; I[ni++] = ny > 0 ? v2 : v1; I[ni++] = ny > 0 ? v1 : v2;      const w3 = nb[w0 * 6 + 1];      if (w3 !== nc && flat[w3 * S + z] & bit) { const v3 = vert(X + dx * 1.5, y, Z + SQ * dx, 0, ny, 0); I[ni++] = v1; I[ni++] = ny > 0 ? v2 : v3; I[ni++] = ny > 0 ? v3 : v2; }    }    return !deep;  };  for (let n = c.bn; n < L.nsrf && nv < maxV && ni < maxI; n++) {    if ((n & 127) === 127 && performance.now() - t0 > ms) { c.bn = n; c.nv = nv; c.ni = ni; return false; }    const k = srf[n], w = (k / S) | 0, z = k - w * S, e = w * 6, i = ci[w], j = cj[w];    if (section && !wedge && j > 0) continue;    const X = (i + j * 0.5) * dx, Z = j * SQ * dx, y0 = z === 0 ? 0 : (z - 0.5) * dx, y1 = (z + 0.5) * dx;    // 横の面: 下のますから続いている面は、下のますのところでまとめて描く    g = 0.5;    for (let d = 0; d < 6; d++) {      const q = nb[e + d] * S;      if (a[q + z] || (z > 0 && a[k - 1] && !a[q + z - 1])) continue;      let e2 = z;      while (e2 + 1 < Nz && a[w * S + e2 + 1] && !a[q + e2 + 1]) e2++;      const f = SIDE[d], ca = f[0], sa = f[1], yt = (e2 + 0.5) * dx, c1x = X + f[2] * dx, c1z = Z + f[3] * dx, c2x = X + f[4] * dx, c2z = Z + f[5] * dx;      const b0 = vert(c2x, y0, c2z, ca, 0, sa), b1 = vert(c1x, y0, c1z, ca, 0, sa), b2 = vert(c1x, yt, c1z, ca, 0, sa), b3 = vert(c2x, yt, c2z, ca, 0, sa);      I[ni++] = b0; I[ni++] = b1; I[ni++] = b2; I[ni++] = b0; I[ni++] = b2; I[ni++] = b3;    }    // 上の面    g = TONE0 + (1 - TONE0) * Math.min(1, (z + 1) * gz);    if (!a[k + 1]) {      let fan = true, nx = 0, ny = 1, nz = 0;      if (flat[k] & 1) {        fan = flatPart(e, i, j, z, X, Z, 1, y1, 1);      } else {        // 段差のそばの面は、段差の向きへ少し傾けて照らす(尾根と溝が明暗で読める)        let gx = 0, gy = 0;        for (let d = 0; d < 6; d++) { const q = nb[e + d] * S + z, h = a[q + 1] ? 1 : a[q] ? 0 : -1; gx += h * DX[d]; gy += h * DZ[d]; }        nx = -0.3 * gx; nz = -0.3 * gy; const l = 1 / Math.hypot(nx, 1, nz); nx *= l; ny = l; nz *= l;      }      if (fan) {        const b = nv;        for (let q = 0; q < 6; q++) vert(X + CORN[q][0] * dx, y1, Z + CORN[q][1] * dx, nx, ny, nz);        I[ni++] = b; I[ni++] = b + 2; I[ni++] = b + 1; I[ni++] = b; I[ni++] = b + 3; I[ni++] = b + 2; I[ni++] = b; I[ni++] = b + 4; I[ni++] = b + 3; I[ni++] = b; I[ni++] = b + 5; I[ni++] = b + 4;      }    }    // 下の面(柱の両端に広がった板の裏など)    if (z > 0 && !a[k - 1]) {      g = 0.42;      if (!(flat[k] & 2) || flatPart(e, i, j, z, X, Z, 2, y0, -1)) {        const b = nv;        for (let q = 0; q < 6; q++) vert(X + CORN[q][0] * dx, y0, Z + CORN[q][1] * dx, 0, -1, 0);        I[ni++] = b; I[ni++] = b + 1; I[ni++] = b + 2; I[ni++] = b; I[ni++] = b + 2; I[ni++] = b + 3; I[ni++] = b; I[ni++] = b + 3; I[ni++] = b + 4; I[ni++] = b; I[ni++] = b + 4; I[ni++] = b + 5;      }    }  }  const ia = c.geo.index as THREE.BufferAttribute;  for (const key of ["position", "normal", "color"]) { const at = c.geo.attributes[key] as THREE.BufferAttribute; at.clearUpdateRanges(); at.addUpdateRange(0, nv * 3); at.needsUpdate = true; }  ia.clearUpdateRanges(); ia.addUpdateRange(0, ni); ia.needsUpdate = true;  c.geo.setDrawRange(0, ni); c.tris = ni / 3; c.phase = 0;  c.cap.visible = section;  if (section) buildCap(c);  return true;}/** 断面の切り口: 枝の中心線(j = 0)を通る面で切った長方形 */function buildCap(c: Cry) {  const L = c.L, dx = c.dx, CP = c.capPos, face = 0.59 * dx;  let o = 0;  const rect = (x0: number, x1: number, y0: number, y1: number) => {    if (o + 18 > CP.length) return;    CP[o++] = x0; CP[o++] = y0; CP[o++] = face; CP[o++] = x1; CP[o++] = y0; CP[o++] = face; CP[o++] = x1; CP[o++] = y1; CP[o++] = face;    CP[o++] = x0; CP[o++] = y0; CP[o++] = face; CP[o++] = x1; CP[o++] = y1; CP[o++] = face; CP[o++] = x0; CP[o++] = y1; CP[o++] = face;  };  for (let w = 0; w < L.nc; w++) {    if (L.cj[w] !== 0 || L.top[w] < 0) continue;    const i = L.ci[w], cx = i * dx;    let z = 0;    while (z <= L.top[w]) {      if (!L.a[w * L.S + z]) { z++; continue; }      let e = z;      while (e + 1 <= L.top[w] && L.a[w * L.S + e + 1]) e++;      const y0 = (z - 0.5) * dx, y1 = (e + 0.5) * dx;      for (const sx of L.wedge && i !== 0 ? [1, -1] : [1]) {        rect(sx * cx - dx / 2, sx * cx + dx / 2, z === 0 ? 0 : y0, y1);        rect(sx * cx - dx / 2, sx * cx + dx / 2, -y1, z === 0 ? 0 : -y0);      }      z = e + 1;    }  }  const pa = c.capGeo.attributes.position as THREE.BufferAttribute;  pa.clearUpdateRanges(); pa.addUpdateRange(0, o); pa.needsUpdate = true;  c.capGeo.setDrawRange(0, o / 3);}/** 水蒸気の濃さを頂点の明るさにする。遠くの値 ρ を基準に、少なくなった所ほど点がまばらになる */function paintField(c: Cry) {  const L = c.L, F = c.fCol, I = c.fIdx, inv = FIELD / Math.max(1e-6, L.rho);  for (let v = 0; v < I.length; v++) { const t = Math.min(0.23, L.d[I[v] * L.S] * inv); F[v * 3] = F[v * 3 + 1] = F[v * 3 + 2] = t; }  (c.fGeo.attributes.color as THREE.BufferAttribute).needsUpdate = true;}

buildMeshは、結晶の表面のますを順に見て、外に出ている面を三角形にする。縦に続く横の面は1枚にまとめ、平らに続く上の面は六角形でなく大きな三角形で埋める。こうして作った24分の1の形を、setInstancesの行列で12回ずつ並べる（InstancedMesh）。上下を返した側はflipの中に入れると、three.jsが面の表裏を入れ替えてくれる。断面を見るときは、手前の半分の行列を外し、切り口の長方形をbuildCapで足す。三角形を書く作業は1コマ1.2 msで区切り、書き終えたときだけGPUへ送る。描く三角形がTRI_MAX（30万個）を超えたら、育てるのを止める。水蒸気の濃さはpaintFieldが高さ0の面の頂点の明るさにする。

図と履歴

中谷ダイヤグラムは3Dの絵の上に2Dで描く。結晶は図を除いた広さの真ん中に写すので、layoutでカメラの中心をずらしている。

図と履歴(src/sims/snow-crystal.ts の 500〜585 行)

/** 画面の配置を決める。横長なら図を左下、縦長なら下に置き、結晶は残りの広さの真ん中に写す */function layout(api: Api, s: S) {  const W = api.w, H = api.h, U = api.u, land = W >= H * 1.05, m = U(10);  let pw: number, ph: number;  if (land) { pw = Math.min(W * 0.38, H * 0.66, U(300)); ph = pw * 0.8; } else { pw = W - 2 * m; ph = Math.min(H * 0.36, pw * 0.62); }  pw = Math.max(pw, Math.min(W - 2 * m, 200)); ph = Math.max(ph, 150);  const x = m, y = H - ph - Math.max(m, U(34)), b = s.box, al = U(30), ab = U(26), at = U(24), ar = U(16);  b.x = x; b.y = y; b.w = pw; b.h = ph; b.x0 = x + al; b.x1 = x + pw - ar; b.y0 = y + at; b.y1 = y + ph - ab;  // 上の端は空けておく(画面の知らせが結晶に重ならないように)  const top = U(land ? 44 : 54);  if (land) { s.freeX = x + pw; s.freeY = top; s.freeW = W - s.freeX; s.freeH = H - top; } else { s.freeX = 0; s.freeY = top; s.freeW = W; s.freeH = y - top; }  s.v.camera.setViewOffset(W, H, W / 2 - (s.freeX + s.freeW / 2), H / 2 - (s.freeY + s.freeH / 2), W, H);}const tx = (b: Box, T: number) => b.x0 + (-T / 30) * (b.x1 - b.x0);const sy = (b: Box, sg: number) => b.y1 - Math.sqrt(clamp(sg, 0, SMAX) / SMAX) * (b.y1 - b.y0);const txi = (b: Box, x: number) => clamp((-(x - b.x0) / (b.x1 - b.x0)) * 30, TMIN, TMAX);const syi = (b: Box, y: number) => { const v = (b.y1 - y) / (b.y1 - b.y0); return clamp(SMAX * v * Math.max(v, 0), SMIN, SMAX); }; /** 中谷ダイヤグラムの大まかな区切りと、今の条件・履歴を描く */function drawDiagram(api: Api, s: S) {  const b = s.box, U = api.u, fs = Math.max(U(10), 9.5), small = fs * 0.92;  api.rect(b.x, b.y, b.w, b.h, { fill: "paper", w: U(1) });  api.text("中谷ダイヤグラム", b.x + U(8), b.y + U(14), { size: fs * 1.05, weight: 600 });  // 板と柱の帯  for (const T of [-4, -10, -22]) { const x = tx(b, T); api.line(x, b.y0, x, b.y1, { w: U(0.7), dash: [U(2), U(3)] }); }  const band = (T: number, sg: number, t: string, o?: { weight?: number }) => api.text(t, tx(b, T), sy(b, sg), { size: small, align: "center", baseline: "middle", weight: o?.weight ?? 400 });  band(-2, 32, "板", { weight: 700 }); band(-7, 32, "柱", { weight: 700 }); band(-16, 36, "板", { weight: 700 }); band(-26, 36, "柱", { weight: 700 });  band(-15.5, 24, "樹枝"); band(-15, 2.2, "角板"); band(-7, 1.2, "角柱"); band(-5.5, 11, "針");  // 軸  api.line(b.x0, b.y1, b.x1, b.y1, { w: U(1) }); api.line(b.x0, b.y0, b.x0, b.y1, { w: U(1) });  for (const T of [0, -10, -20, -30]) { const x = tx(b, T); api.line(x, b.y1, x, b.y1 + U(3), { w: U(1) }); api.text(fmt(T, 0), x, b.y1 + U(13), { size: small, align: "center", font: "mono" }); }  api.text("温度 °C", b.x1, b.y1 + U(23), { size: small, align: "right" });  for (const v of [1, 5, 10, 20, 40]) { const y = sy(b, v); api.line(b.x0 - U(3), y, b.x0, y, { w: U(1) }); api.text(String(v), b.x0 - U(5), y + small * 0.34, { size: small, align: "right", font: "mono" }); }  api.text("過飽和度 %", b.x + U(4), b.y1 + U(23), { size: small });  // 水の飽和(雲の水滴に接している空気)  const pts: number[] = [];  for (let T = -0.5; T >= -30; T -= 1) pts.push(tx(b, T), sy(b, sigmaWater(T)));  api.poly(pts, { w: U(1.6) });  { const q = -25.5; api.text("水の飽和", tx(b, q), sy(b, sigmaWater(q)) + small * 1.5, { size: small, align: "center" }); }  // 履歴  if (s.grab === "draw" && s.draft.length > 1) {    const d: number[] = [];    for (const p of s.draft) d.push(tx(b, p[1]), sy(b, p[2]));    api.poly(d, { w: U(2.6) });  } else if (s.plan) {    const P = s.plan, d: number[] = [];    for (let i = 0; i < P.n; i++) d.push(tx(b, P.T[i]), sy(b, P.s[i]));    api.poly(d, { w: U(2.2) });    const o = { T: 0, s: 0 };    for (let m = 0; m <= P.end + 1e-3; m += 10) { planAt(P, m, o); api.dot(tx(b, o.T), sy(b, o.s), U(2.6)); }  } else if (s.trail.length > 4) {    const d: number[] = [];    for (let i = 0; i + 1 < s.trail.length; i += 2) d.push(tx(b, s.trail[i]), sy(b, s.trail[i + 1]));    api.poly(d, { w: U(1.6), dash: [U(3), U(3)] });  }  // 今の点: 地の色で縁どった丸と、つかめる印の点線の輪  const c0 = s.cry[0], T = s.grab === "handle" || !c0 ? s.T : c0.T, sg = s.grab === "handle" || !c0 ? s.sig : c0.s, X = tx(b, T), Y = sy(b, sg);  if (s.mode === "twin" && s.cry[1]) { const c1 = s.cry[1]; api.circle(tx(b, c1.T), sy(b, c1.s), U(4), { fill: "paper", w: U(1.4) }); }  api.circle(X, Y, U(5), { fill: "ink", stroke: "paper", w: U(1.5) });  api.circle(X, Y, U(11), { stroke: "paper", w: U(3) });  api.circle(X, Y, U(11), { w: U(1.4), dash: [U(3), U(3)] });} /** 結晶の絵に重ねる文字: 物差しと名札 */function drawLabels(api: Api, s: S) {  const U = api.u, v = s.v, d = v.orbit.dist, px = api.h / (2 * Math.tan((v.camera.fov * Math.PI) / 360) * d);  let L = 1;  for (const q of [0.05, 0.1, 0.2, 0.5, 1, 2, 5]) if (q * px <= s.freeW * 0.28) L = q;  const x0 = s.freeX + s.freeW / 2 - (L * px) / 2, y = s.freeY + s.freeH - U(18), x1 = x0 + L * px;  for (const [w, c] of [[U(5), "paper"], [U(2), "ink"]] as const) { api.line(x0, y, x1, y, { w, stroke: c }); api.line(x0, y - U(5), x0, y + U(5), { w, stroke: c }); api.line(x1, y - U(5), x1, y + U(5), { w, stroke: c }); }  api.label(`${L} mm`, (x0 + x1) / 2, y - U(9), { size: U(12), align: "center" });  if (s.mode === "twin" && s.cry.length === 2) ["A", "B"].forEach((t, n) => { const c = s.cry[n], p = v.project([c.group.position.x, 0, -(c.L.rmax + 4) * c.dx]); if (p.on) api.label(t, p.x, Math.max(U(18), p.y - U(8)), { size: U(15), weight: 700, align: "center" }); });} function readouts(api: Api, s: S) {  const c = s.prev && s.grab === "handle" ? s.prev : s.cry[0];  if (!c) return;  const L = c.L, d = (2 * L.rmax + 1) * c.dx, th = (2 * L.zmax + 1) * c.dx, wr = sigmaWater(c.T), lines: string[] = [];  lines.push(`経過  ${fmt(s.grab === "handle" ? L.t * c.h : s.minute, 1)} 分${s.grab === "handle" ? "(粗い下書き)" : ""}`);  lines.push(`温度と過飽和度  ${fmt(c.T, 1)} °C  ${fmt(c.s, 1)} %(水の飽和の ${fmt(c.s / Math.max(wr, 0.1), 2)} 倍)`);  lines.push(`大きさ  差し渡し ${fmt(d, 2)} mm  厚さ ${fmt(th, 2)} mm`);  lines.push(`形  ${s.shape}  先の伸びる速さ ${fmt(c.rate, 1)} mm/時`);  if (s.mode === "twin") lines.push(`A と B の違い  結晶のますの ${fmt(s.diff, 1)} %`);  if (s.mode === "full") lines.push(`6 本の枝  長さ ${fmt(s.arms[0], 2)}〜${fmt(s.arms[1], 2)} mm  回して重ねると ${fmt(s.sym, 1)} % が違う`);  api.readout(lines);}

横軸は温度で、縦軸は過飽和度の平方根にした。水の飽和の曲線は−5 °Cの5%から−30 °Cの34%まで広がるので、平方根の目盛りなら低い所も読める。板と柱の区切りの点線は−4 °C、−10 °C、−22 °Cに引いた。資料にある−2、−5、−15、−30 °C付近の記述から決めた目安で、計算には使っていない。

点を動かす、なぞる、種を置く

図の点をつかむとrunPreviewが半径32ますの粗い格子で結晶を育て、形をすぐに見せる。離すと、本来の格子で育て直す。

操作(src/sims/snow-crystal.ts の 675〜710 行)

function startDraw(api: Api, s: S) { s.grab = "draw"; s.draft = []; s.draftT0 = performance.now(); addDraft(api, s); }function addDraft(api: Api, s: S) {  const b = s.box, p = api.pointer;  s.draft.push([(performance.now() - s.draftT0) / 1000, txi(b, p.x), syi(b, p.y)]);}/** なぞり終えたら、指の動いた時間をそのまま履歴の時間にする(ゆっくりなぞった所に長くとどまる) */function finishDraw(api: Api, s: S) {  const D = s.draft, span = D.length ? D[D.length - 1][0] : 0;  const moved = D.length > 1 && Math.hypot(tx(s.box, D[0][1]) - tx(s.box, D[D.length - 1][1]), sy(s.box, D[0][2]) - sy(s.box, D[D.length - 1][2])) > api.u(10);  s.grab = null;  if (!moved || span <= 0) {    // たたいただけなら、その条件で育て直す    s.plan = null; s.T = D[0]?.[1] ?? s.T; s.sig = D[0]?.[2] ?? s.sig; s.rs = 1; syncKnobs(api, s); if (s.curRs === s.rs) restart(api, s);    return;  }  // 点の数を 120 までに間引く  const st = Math.max(1, Math.floor(D.length / 120)), pts: [number, number, number][] = [];  for (let i = 0; i < D.length; i += st) pts.push([(D[i][0] / span) * PLAN_MIN, D[i][1], D[i][2]]);  const last = D[D.length - 1];  if (pts[pts.length - 1][0] < PLAN_MIN) pts.push([PLAN_MIN, last[1], last[2]]);  s.plan = makePlan(pts);  if (s.mode === "one" && s.curRs !== PLAN_RS) s.rs = PLAN_RS; else restart(api, s);  api.note("なぞった道すじのとおりに育てます", 2200);}function syncKnobs(api: Api, s: S) { api.set("T", Math.round(s.T * 10) / 10); api.set("sig", Math.round(s.sig * 10) / 10); }function placeSeed(api: Api, s: S) {  if (s.mode !== "full") { if (!s.noted.seed) { s.noted.seed = true; api.note("種を足せるのは「ゆらぎ」の場面です(ほかの場面は 6 本がそろう計算です)", 2600); } return; }  const c = s.cry[0], L = c.L, q = s.v.onPlane(api.pointer.x, api.pointer.y, new THREE.Plane(new THREE.Vector3(0, 1, 0), 0));  if (!q) return;  const j = Math.round(q.z / (SQ * c.dx)), i = Math.round(q.x / c.dx - j / 2);  if (hexDist(i, j) > L.R - 6) return;  addSeed(L, i, j, 2);  L.rmax = Math.max(L.rmax, hexDist(i, j) + 2);  c.done = false;  api.note("新しい種を置きました。近くの枝と水蒸気を取り合います", 2200);}

点のない所から図をなぞるとaddDraftが時刻と条件を記録し、finishDrawが全体を36分に引き伸ばして履歴（Plan）にする。動かした時間を残しているので、指を止めた所に結晶が長くとどまる。履歴をたどる間はupdateCoefが16刻みごとに係数を更新し、 $ρ$ が変わればsetRhoで場全体を同じ倍率にする。本物の空気では拡散が成長よりはるかに速いので、遠くから新しい濃さがしみ込むのを待たない。「ゆらぎ」の場面で空いた所をたたくとplaceSeedが種を足し、向かい合う枝が水蒸気を取り合う。

画面の部品(描画・場面の切り替え・画質の自動調整・GPU の網点)は全記事で共通です。シミュレーションの共通の仕組みで解説しています。

プログラム全体を見る(src/sims/snow-crystal.ts・844 行)

読み込んでいます

この記事は、人が監修しながらAIを活用して執筆しています。執筆・監修: 佐々木航(株式会社Game)

自然のほかの動き

公転する惑星は、太陽へ落ち続けながら横へ進んでいる

地球は1秒に約30 km進むあいだに、太陽の方へ約3 mm落ちている。

ブラックホールに近づいた星は、なぜ回り続けられないのか

太陽のまわりなら、どれほど近くても円を描いて回れる。

脳も設計図もないのに、粘菌の網や動物の模様が形になるわけ

迷路を満たした粘菌は、入口と出口の餌をいちばん短い道の管で結ぶ。

音で揺らした水面は、揺れの半分の周波数で模様を描く

皿の水を50 Hzで上下に揺らすと、水面には25 Hzで上下する格子の模様が立つ。