旋转矩阵

背景

蒙特卡洛仿真中经常会遇到不同方向的界面，当运动的粒子撞到这个界面的时候，我们需要根据立体角的权重去抽样粒子被反射后的运动方向。这个时候，当我们想要统一处理不同方向的界面时，就需要用到旋转矩阵。

当这个平面处于\((x, y)\)平面时，它的法向量就是\((0, 0, 1)\)。这个平面进行黑体辐射随机发射一个粒子时，天顶角\(\theta\)和方位角\(\phi\)满足如下的抽样公式： \[ R_1, R_2\in [0, 1], \sin\theta = \sqrt{R_1}, \phi = 2\pi R_2 \] 在笛卡尔坐标里，抽样出来的粒子运动方向的就是： \[ (x, y, z) = (\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta) \] 假如现在这个平面是个斜的，它的法向量是\((n_x, n_y, n_z)\)。这个平面进行黑体辐射，应该怎么去抽样粒子的运动方向呢？我们知道，发射过程相对于各自平面的法向量来说总是等价的。于是我们可以这样构建抽样的过程，首先还是对\((x, y)\)平面进行抽样，满足上述的抽样公式，即相对法向量\((0, 0, 1)\)抽样出来了\((\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta)\)。这个时候我们把法向量和这个方向向量同时进行旋转，当法向量旋转到\((n_x, n_y, n_z)\)的时候，方向向量也就到了它应该在的位置。

旋转变换

于是现在问题就变成了，怎么把描述把\((0, 0, 1)\)旋转到\((n_x, n_y, n_z)\)的过程呢？我们知道，旋转是一个线性变换，而线性变换和矩阵乘法是一一对应的。也就是说，只要找到这个旋转变换对应的矩阵，乘上抽样出来的向量，就完成了这个旋转变换。那么怎么找到这个旋转矩阵呢？我们还知道，线性变换是具有可加性的，也就是说我们可以分成几步去把\((0,0,1)\)挪到\((n_x, n_y, n_z)\)上，计算了每一步的旋转矩阵以后，把这些旋转矩阵都乘起来，就可以得到最后的总旋转矩阵了。这个旋转矩阵乘上\((0, 0, 1)\)等于\((n_x, n_y, n_z)\)，乘上抽样出来的方向向量，也会把方向向量旋转到到它对应的位置上。

初始向量是\((0,0,1)\)，目标向量是\((r, \theta,\varphi)\)，如下图所示，我们可以把这个旋转过程分成两步。第一步，首先绕着\(y\)轴逆时针转\(\theta\)，把\((0,0,1)\)转到蓝色的虚线位置处。接下来，绕\(z\)轴逆时针旋转\(\varphi\)，把蓝色的虚线再转到棕色的目标向量处。我们只要分别得到这两步旋转的旋转矩阵就可以了。相当于 \[ R(\theta, \varphi) = R(\varphi)R(\theta) \]

第一步旋转

对于绕\(y\)轴旋转的过程，我们可以用一个一般向量来推导旋转矩阵。假设向量的初始的坐标是\((x_0, y_0, z_0)\)，旋转后的坐标是\((x^\prime, y^\prime, z^\prime)\)。显然，旋转前后\(y_0 = y^\prime\)。同时，几何关系上存在 \[ \begin{aligned} x^\prime / r &= \sin(\theta_0 + \theta) = \sin\theta_0 \cos\theta + \cos\theta_0 \sin\theta\\ z^\prime / r &= \cos(\theta_0 + \theta) = \cos\theta_0\cos\theta - \sin\theta_0\sin\theta \end{aligned} \] 整理可以得到， \[ \begin{aligned} x^\prime &= \cos\theta x_0 + \sin\theta z_0 \\ y^\prime &= y_0 \\ z^\prime &= - \sin\theta x_0+\cos\theta z_0 \end{aligned} \] 整理成矩阵的形式 \[ \left[\begin{array}{c}x^\prime\\y^\prime\\z^\prime\end{array}\right]=\begin{bmatrix}\cos\theta&0&\sin\theta\\0&1&0\\-\sin\theta&0&\cos\theta\end{bmatrix}\left[\begin{array}{c}x_0\\y_0\\z_0\end{array}\right] \] 于是就得到了第一步旋转的旋转矩阵。

第二步旋转

同理可以推导绕\(z\)轴旋转的矩阵，此时显然\(z^\prime = z_0\)，同时满足 \[ \begin{aligned} x^\prime / r &= \cos(\varphi_0 + \varphi) = \cos\varphi_0 \cos\varphi - \sin\varphi_0 \sin\varphi\\ y^\prime / r &= \sin(\varphi_0 + \varphi) = \sin\varphi_0\cos\varphi + \cos\varphi_0\sin\varphi \end{aligned} \] 整理可以得到， \[ \begin{aligned} x^\prime &= \cos\varphi x_0 - \sin\varphi y_0 \\ y^\prime &= \sin\varphi x_0+\cos\varphi y_0 \\ z^\prime &= z_0 \end{aligned} \] 整理成矩阵的形式 \[ \left[\begin{array}{c}x^\prime\\y^\prime\\z^\prime\end{array}\right] =\begin{bmatrix}\cos\varphi&-\sin\varphi&0\\\sin\varphi&\cos\varphi&0\\0&0&1\end{bmatrix}\left[\begin{array}{c}x_0\\y_0\\z_0\end{array}\right] \]

两步旋转

于是总的旋转矩阵就可以写成 \[ R(\theta, \varphi) = R(\varphi)R(\theta) = \begin{bmatrix}\cos\varphi&-\sin\varphi&0\\\sin\varphi&\cos\varphi&0\\0&0&1\end{bmatrix}\begin{bmatrix}\cos\theta&0&\sin\theta\\0&1&0\\-\sin\theta&0&\cos\theta\end{bmatrix} = \begin{bmatrix}\cos\varphi\cos\theta & -\sin\varphi & \cos\varphi \sin\theta \\ \sin\varphi\cos\theta & \cos\varphi & \sin\varphi \sin\theta \\ -\sin\theta & 0 &\cos\theta \end{bmatrix} \] 如果目标法线的方向向量是\((n_x, n_y, n_z)\)，那么存在 \[ \sin\theta = \sqrt{n_x^2 + n_y^2}, \cos\varphi = n_x / \sqrt{n_x^2 + n_y^2} \] 代入上面就可以算出来旋转矩阵\(R\)了。此时如果这个平面进行了一次黑体发射，抽样粒子运动方向的规则就应该是 \[ (x, y, z) = R \cdot \left[\begin{array}{c}\sin\theta\cos\phi\\\sin\theta\sin\phi\\\cos\theta\end{array}\right] \] 其中 \[ R_1, R_2\in [0, 1], \sin\theta = R_1, \phi = 2\pi R_2 \]

注意

旋转是有方向和顺序依赖的，不同的旋转路径会导致不同的结果。比如对于下面这个例子，把\((0,0,1)\)法向量旋转到\((-1,0,0)\)上，如果沿\(y\)轴顺时针旋转，那么原先指向\((1,2,3)\)的向量会旋转到\((-3,2,1)\)的位置。但如果先沿\(y\)轴逆时针旋转\(\pi/2\)，再沿\(z\)轴逆时针旋转\(\pi\)，此时原先指向\((1,2,3)\)的向量会旋转到\((-3,2,1)\)的位置。当然，不论怎样的旋转路径，法向量和方向向量之间的相对位置关系是不变的，对于这个例子来说，不论怎么旋转这法向量和方向向量之间的点积都是3。它不会影响粒子发射方向的相对几何关系，只要在模拟中始终一致地使用同一种旋转方式，模拟结果就应该是可靠的。