如何简单渲染一个彩虹
本文参考了我一位学长的同名项目，和 Philip Laven 教授 (?.?.1933-2.7.2024) 的 The optics of a water drop-Mie scattering and the Debye series。Laven 教授生前和合作者发了不少关于大气光学现象模拟的文章，对虹、霓、日晕等现象做了详尽解释。关于彩虹的主要成果总结在 Physically Based Simulation of Rainbows 一文。
1. 彩虹的物理形成机制-几何光学视角
1.1 仅考虑入射平面：单圆柱散射
彩虹是大气中的悬浮水珠散射、色散入射的阳光，在人眼中形成的虚像。彩虹形成于太阳向水汽团的延长线上，一般以彩虹圆弧为中心，视角在 40°~42° 的彩虹光线最为强烈。对于这种角度，阳光在水滴中进行一次全反射后偏转射出。
彩虹的形成原因非常复杂，水滴的尺寸和光的波长接近，导致了瑞利散射，使彩虹的构成异常复杂。本章节仅从几何光学视角考虑本问题。
图 1，内反射 k=1 时的示意图。k=0，即光线首次接触水滴内表面时的出射光线未画出。
假设水滴为半径为 R 的均匀球体，空气折射率 $n_{air}$，水折射率 $n_{water}$ ，二者都会随射入的光的波长而变化，即 $n_{air}(\lambda)$，$n_{water}(\lambda)$。太阳光可以近似为平行射入的圆偏振光。
以球心为原点建立直角坐标系，并假设光沿着 $-x$ 方向入射，这样只需要考虑 $x > 0$ 部分的光照。又由于对称性，只须考虑 y>=0 部分的即可，y<0 部分做镜像即可得到。
假设光线在 $(R\cos \theta, R\sin\theta)$ 位置入射，记光线入射球的位置为 $A_0$，之后球内光线每内反射一次的位置为 $A_k(k>=1)$，如下图。
我是图
假设 $A_0$ 处的入射角为 $i=\theta$，折射角为 $r$，则容易得到
$$
A_k=(R\cos\phi, R\sin\phi), \phi=\theta+k(\pi-2r)\tag{1}
$$
对于 $k>=1$，入射角大小恒定为 $r$，折射角为 $\theta$，则由折射定律，
$$
\frac{n_{water}}{n_{air}}=\frac{\sin r}{\sin \theta}\tag{2}
$$
考虑菲涅尔方程（Fresnell Equation）
对于发生在 $n_1$，$n_2$ 两个界面上的折反射，假设光线从介质 $n_1$ 射入 $n_2$，入射角 $\theta_i$，折射角 $\theta_r$，反射角 $\theta_t$，入射光功率为一个单位，容易得到
反射光 s 偏振部分功率 $$R_s=\left[\frac{\sin(\theta_t-\theta_i)}{\sin(\theta_t+\theta_i)}\right]^2\tag{3}$$
反射光 p 偏振部分功率 $$R_p=\left[\frac{\tan(\theta_r-\theta_i)}{\tan(\theta_r+\theta_i)}\right]^2\tag{4}$$
对于圆偏振光（太阳光, 即本文情况），总功率 $$R=(R_s+R_p)/2\tag{5}$$
不管什么时候，都有折射光功率 $$T=1-R\tag{6}$$
方便起见, 我们记空气到水的反射率、折射率为 $R_1,T_1$，反之为 $R_2,T_2$。光线从空气入射水滴，在水滴内部经历了 $k\geq 1$ 次反射后射出，容易知道，第 k 级射出水滴的光线功率 $dP_{out,k}$ 和最初入射光功率 $dP_{i,0}$ 的比值
$$
\frac{dP_{out,k}}{dP_{i,0}}=\frac{dP_{t,k}}{dP_{i,0}}=T_2R_2^{k-1}T_1 \tag{7}
$$
水的折射率是温度和光的波长的函数。参考refradtiveindex.info的水和空气在特定温度下的折射率插值表，简单记为
$$
n_{water}=n_{water}(\lambda)\tag{8}
$$
$$
n_{air}=n_{air}(\lambda)\tag{9}
$$
联立 $(1)$~$(9)$ 式，我们可以写一个程序看看圆柱体对光的色散效果。详细结果见仓库的cylinder.py就行。数值绘图和分析见1.3.
1.2 球的散射
让我们把问题拓展到一个球，也即真实的水滴的散射情况。我们假设光是从
$$
P(R,\theta, \varphi)=(R\sin\varphi\cos\theta,R\sin\varphi\sin\theta,R\cos\varphi)\tag{10}
$$
射入的，入射方向沿 $-x$，也即 $\theta = \pi /2$。此时P和x轴能确定唯一平面，也即本入射光线的所在平面$MPN$，记为$\Omega$，其法线方向
$$
\vec{\omega}=\hat{x}\times\overrightarrow{OP}\tag{11}
$$
在$\Omega$平面上，$\angle MOP$就是零级入射角$\theta$。由余弦定理，
$$
\angle MOP=\arccos \frac{\overline{OM}^2+\overline{OP}^2-\overline{MP}^2}{2 \overline{OM}\cdot\overline{OP}}\tag{12}
$$
容易得到，
$$
\angle MOP=\arccos(\sin\varphi\cos\theta)\tag{13}
$$
此后对于这条光线，将1.1部分中的入射角$\theta$替换成这个角就可以了。
在做进一步详细计算之前，让我们计算对入射角$\theta$的光线，第k级的出射光线和它的夹角（下面统称为偏转角$\Delta k$）是多少。
在光线传播的平面上，我们认为入射光的方向角为$\pi$,首先容易得到，
0级的反射光（即水滴上的反射光）的偏转角
$$
\Delta_0=-(\pi-2\theta)\tag{14}
$$
对于k=1的光线，它们各自在水滴内部反射了$k-1$次，方位角改变了$\pi-2r$ ；又2次穿过水-空气界面，方位角改变 $\theta-r$，于是
$$
\Delta_k=2(\theta-r)+(k-1)(\pi-2r)\
=2\theta-2kr+(k-1)\pi\tag{15}
$$
我是图
三维空间中，入射光的方向均为$-\hat{x}=(0,\pi,\pi/2)$。垂直平面的法向量为$-\hat{y}$，不在竖直平面入射的光，可以把其所在平面$\Omega$的法向量$\vec{\omega}$旋转到$-\hat{y}$方向。
由上述一系列式子得到，在$P(R,\theta, \varphi)$入射的光，第k级的射出的光的方向角是
$$
\text{我是公式}
$$
据此，终于可以完整计算水滴折射平型光后射出的光的强度分布。
1.3 数值模拟及结果分析
1.3.1 圆柱散射
我是图
能看到，1，2级的出射光的强度远大于其他的（包括0级）。
分析一下，为什么一级虹（彩虹）的红光在外紫光在内，用2级最大去分析。参考这个彩虹中文维基
1.3.2 球散射
2. 考虑散射对光的影响
实际的水滴的直径很小
（真的吗，这个需要查证）
直径接近光的波长时无法忽略散射现象。
（那么是rui’li