什么错误？
在所有过去的冰晕模拟中，当设置 Plate 取向的 c-axis 倾斜程度为正态分布的时候，采样的结果是错误的。我们过去都假设冰晶的 c-axis 倾斜程度为正态分布，听起来没什么问题，但是问题是 Plate 取向是一种“退化”的取向，c-axis 倾斜程度在两个自由度上分别是正态分布，它们的共同作用应该导致 c-axis 的分布是某种接近瑞利分布的东西。
在所有过去的冰晕模拟中，Plate + Gaussian 的生成方式全部为：
- 先放置好冰晶，c 轴朝上对应 z 轴方向，x、y 轴与冰晶底面的水平面平行；
- 随机让冰晶绕着 z 轴旋转一个随机的均匀分布角度，模拟冰晶的横滚（ψ，roll）;
- 让冰晶绕着 x 轴随机旋转一个从高斯分布采样的角度，模拟冰晶的晃动（θ，zenith）；
- 让冰晶绕着 z 轴再旋转随机的均匀分布角度，模拟方位角的随机性（φ，azimuth）；
这种做法相当于在一个自由度上生成了高斯分布，然后通过一个旋转操作试图让这个分布变成二维高斯分布。这很明显是不正确的，对概率密度函数的操作不能变成对随机变量的操作。所有过去的冰晕模拟中都存在采样的错误，包括但不限于：
HaloSim <= 3.6.1HaloRay <= 4.1.0HaloPoint <= 2.0
但是我们的数学都很差没有发现这个问题，直到最近 爪爷在 Lumice 4.1.13 版本中才修复了这个问题。
应该怎么样？
如果我们还是认为需要在球面上采样高斯分布，在球面上采样角度 θ ，那么就需要对球面的面积元进行修正：p_sphere(θ) ∝ p_distribution(θ) · sin(θ)，其中p_distribution(θ) ∝ exp(-(θ-μ)² / 2σ²)。即使不是高斯分布，也需要修正，否则采样点会过度集中在天顶附近，导致冰晶的 c-axis 倾斜程度分布存在一个奇点。
简单验证
我们假设在低雷诺数下，气动力矩将晶体恢复至平衡取向（板状晶体的 c-axis 趋向竖直，柱状晶体的 c 轴趋向水平）。布朗运动和大气湍流引入小角度涨落，冰晶的运动是一种稳定模式。为了直观地验证上述理论（即 Plate 冰晶的倾角大小不应该是一个简单的、在 0 处概率有奇点的高斯分布），我们设计了一个基于物理过程的数值模拟实验（代码参考。
1. 物理模型构建 我们没有使用纯数学的随机采样，而是模拟了一个浮在水面上的圆形平板，以此类推冰晶在空气中下落时受到的随机扰动。在这个模型中：
- 随机扰动：将空气分子的碰撞或气流的不规则扰动等效为“随机掉落的雨滴”。雨滴的下落遵循泊松过程，撞击在平板上的位置和雨滴大小（质量）均是随机的。
- 双轴独立受力：雨滴撞击会使得平板在 $x$ 轴和 $y$ 轴方向上同时产生随机的转矩（Torque）。这意味着平板在两个相互垂直的自由度上的晃动是独立发生的。
- 回复与阻尼：平板同时受到类似弹簧的回复力矩（模拟冰晶在气流中保持水平的气动力学趋势）以及粘滞阻尼（模拟空气阻力）。
2. 倾角演化与稳态分析 通过对该物理过程进行长时间步的积分，我们记录了平板绕 $x$ 轴和 $y$ 轴的倾角时间序列（$\theta_x$ 和 $\theta_y$）。 在去除了初始的瞬态过程后，我们提取了平板处于稳态时的总倾角大小 $|\theta| = \sqrt{\theta_x^2 + \theta_y^2}$。如果旧版冰晕模拟中的“Plate + Gaussian”采样方法是正确的，那么 $|\theta|$ 应该在 0 度附近出现概率最大值。实际上我们发现，在角度为 0 的区域附近，冰晶的取向是很难取到的。
3. 拟合结果与结论 我们将模拟得到的 $|\theta|$ 稳态分布直方图与三种理论分布进行了曲线拟合：
- 瑞利分布 (Rayleigh distribution)：基于 $\theta_x$ 和 $\theta_y$ 为独立同分布的高斯随机变量的数学推导。
- sin(x) * 正态分布：基于球面立体角面积元修正的分布。
- 高斯混合模型 (GMM)：作为对照组。
实验结果清晰地表明，稳态下平板的二维倾角散点 $(\theta_x, \theta_y)$ 在中心点（0, 0）附近虽然密集，但当我们计算倾角大小 $|\theta|$ 的概率密度时，其分布在 0 度处趋近于 0，并在一个大于 0 的角度上达到峰值。这与sin(x) * 正态分布高度吻合。
这一物理现象简单而直接地证明了：现实中受随机扰动的 Plate 冰晶，其总倾角大小的分布绝不是一个以 0 为中心的一维高斯分布。旧版冰晕模拟中单纯在一维生成高斯分布然后再绕 Z 轴旋转的做法，导致冰晶过度集中在天顶（0度）附近，在物理和统计学上都是完全不成立的。
为什么我们一直都没有发现？
这里出现了一个很微妙的事情，这样的错误还是有很明显的影响的。但是在为数不少的情况下，这个错误反而能带来看起来正确的结果。首先我们要来对比一下错误的高斯分布采样和正确的差异。我们以22度幻日为例。
错误的幻日，使用 HaloRay 4.1.0 生成
正确的幻日，使用 Lumice 4.1.13 生成
对于我们这些天天看着中高质量 Display 的来说，我们会认为前者的幻日是正确的，因为这个幻日呈现出典型的三角形，拖出来一个很长的尾巴以及幻日环，非常典型的高质量 Plate Display 的现象。
但这个 bug 也产生了很多错误现象。比如很多日常生活中常见的低质量幻日，可能更类似后者，整个幻日长得很宽，而且没有幻日环。在此之前有很多人怀疑这里面可能存在问题了，因为幻日环在原来的模拟中有点太容易出现了，很大的冰晶晃动角度都能有幻日环，而且这个模拟出现的幻日环很细，被包裹在上下方向模糊的光晕里面，俗称纸包鸡，现实中幻日环似乎不长这样。锥晶幻日在过去的模拟中总是出现莫名其妙的V字形状，也是同样的原因。此外这个 bug 不仅仅影响幻日，还会影响 Column 的模拟，当 Column 取向的 c 轴晃动幅度标准差很大的时候，会出现幻日现象，这也是很不准确的，而且问题的原因是同一个。但大家似乎都以各种理由忽视这些问题，比如“卷云实际冰晶质量差”，或者“锥晶幻日冰晶取向很差衍射太明显”，或者 HaloRay 本身有 bug 等等。
那么现在问题又来了，既然很多中高质量幻日 Display 里面，幻日长得更像错误的模拟呈现的结果，那么说明在这些 case 下的 Plate 冰晶取向还是很偏离高斯分布的。那么应该是什么情况呢？
峰回路转
纯正且正确的高斯分布模拟出来的幻日太“圆滑”了，缺少典型的三角形结构和拖出来的一个很长的尾巴。这意味着实际上确实存在不少的冰晶 c 轴集中在天顶附近小幅度晃动的情况，而且在此之外还有少量大角度晃动的冰晶存在使得幻日在垂直方向展宽。也就是说我们需要修改高斯分布，让分布的核心更加集中，但尾部又要拉长。
具体为什么会这样呢？可能是因为实际的卷云里面冰晶的大小、形状都不均匀，导致了冰晶的取向分布不均匀。最典型原因的就是同一片卷云里面冰晶的尺度差异比较大。当冰晶尺度特别小（比如只有几微米直径）的时候，冰晶下落速度很慢，气动力回复力矩不足，冰晶会大幅度摆动。当冰晶直径很大的时候（比如大于100-150微米直径），冰晶后方由于高雷诺数而产生涡流，会让冰晶周期性正弦摆动；另一部分冰晶可能会螺旋下落，形成一种“固定倾角”的特殊取向（这种取向早期被用于解释 Bottlinger 光环的出现）。而介于两者之间，直径几十微米的冰晶，由于下落过程比较稳定，会呈现较好的取向，c 轴天顶角特别稳定，只在几度之内摆动，呈现为一些方差未知的高斯分布。
我们这里假设对于特别小的冰晶，由于回复力矩不足导致摆动角度很大，这个范围内冰晶摆动角的标准差与冰晶的尺度呈现指数衰减的关系。特别大的冰晶则简化考虑认为取向也是高斯分布的，摆动角度的标准差与冰晶的尺度成线性增加的关系。两者之间的冰晶，由于下落过程比较稳定呈现为一些方差未知（可以人为规定）的高斯分布。我们假设冰晶尺度是 gamma 分布的，那么这个数学上怎么得到闭合的冰晶取向的 c 轴天顶角分布就较难计算了，我们就工程上近似一下。
下面的 html 应用就展示了这一个过程。其中 shape_k 和 scale_theta 代表了冰晶粒径 gamma 分布的参数，这两个参数变大都会让大冰晶变多，在冰晶尺度小于正弦摇摆的程度下都会让冰晶变得更稳定。这个应用输出的 json 文件可以直接导入 Lumice 4.1.13 版本进行冰晕模拟，你会看到用很多个高斯分布的冰晶叠加出来的三角形的幻日。
Halo Config Generator
配置冰晶参数并生成 halo_config.json 渲染配置文件
很多高斯分布的冰晶叠加出来的幻形的幻日
这个多个高斯叠加的问题还有一个更加优雅的处理方法，我们可以直接用拉普拉斯分布来模拟这种情况。当高斯分布的方差来源于一个指数分布的时候，这些高斯分布平均起来就是一个拉普拉斯分布。所以现在新的 Lumice 版本添加了拉普拉斯分布的支持。实际的卫星观测证据也支持冰晶 c 轴天顶角的指数分布。
拉普拉斯分布的冰晶的幻日
这种幻日的特点就是没有“纸包鸡”，但又保留了三角形，非常的完美，完全可以应对大多数质量中等的幻日 display 了。对于高质量的幻日，再添加一个很小幅度摆动的高斯分布，基本就能解释所有的幻日、幻日环的特征了。
总结
Lumice 是目前世界上唯一一个数学上正确的冰晕模拟器，而且计算效率很高，功能也逐渐完善了，都给我狠狠地用。
