今天继续来粗略看看Sig23这篇计算表面反射的全波参考模拟器paper。

关键词：图形学入门、波动光学渲染、BEM、AIM、BRDF
Wave Optics Render, Full Wave Reference Simulator

原文：https://zhuanlan.zhihu.com/p/1471147574

前言

稍微总结一下目前的理解。

基于波动光学的计算机图形学对图形学研究者而言无疑是巨大的挑战。它不仅要在计算上处理电磁波，还涉及一些量子力学的理论，想要打通物理学和图形学之间的隔阂属实不易。在基于波动光学的图形学中，光在介质中的传播不再只是直线，而是会因为波长不同，表现出自旋、偏折和衍射等多种特性。从肥皂泡、偏振镜、油膜光泽或糖分介质导致的旋光再到无线电辐射如何在城市建筑之间传播，都离不开波动光学理论。

在知乎上也有很多大佬分享了非常详细的波动光学渲染理论基础教学。

• 波动光学渲染方程 – Heskey0的文章 – 知乎
• 高中生都能看懂的波动光学渲染（PLT） – Aegle的文章 – 知乎
• 波动光学渲染（入门1）：算法过程概述 – 索罗特的文章 – 知乎

但是奈何阅读门槛都较高，且距离实际应用太远。波动光学渲染涉及的数学物理知识太多，就算是一行一行细看闫老师和他博士生的paper，也是寸步难行。

图形学的圣杯是光线追踪，而全局波动光学渲染则是图形学圣杯中的圣杯。

在Sig21上，Shlomi提出了将路径追踪与物理光学相结合实现电磁辐射传播的真实模拟。目前已经有不少方法计算电磁波的传输，从高精度但计算密集的波动求解器到快速但不够精确的几何光学方法。如下图所示，展示了当前电磁波传输的计算方法。左边最准确，右边最快。

波动求解器（Wave Solvers）专注于求解麦克斯韦方程的精确解，但对于大型场景来说并不实用，一般用FDTD、BEM或FEM来做。

PO基于高频近似的电磁波计算，但是对于可见光这个频率来说其实勉强够用了。PS：[Xia 2023]的黑狗毛就属于Physical Optics方法。而本文的全波参考应该还是属于Wave Solvers的方式，但是在BEM的基础上使用了PO的一些思想（等效电流等）。

除了PO，还有一种方法介于Wave Solvers和Geometrical Optics之间，称为Hybrid GO-PO，我个人觉得应该叫做几何光学-物理光学混合方法。统一衍射理论（Uniform Theory of Diffraction, UTD）将衍射效应纳入几何光学来计算高频条件下的电磁波传输。个人理解，UTD通过计算绕射系数来补偿几何光学射线边界条件的不足，也就是说几何光学的射线也可以转弯了。这种操作在雷达探测天线设计领域非常实用。除了UTD，Hybrid GO-PO还涉及一种叫射线发射和反弹方法（Shooting and Bouncing Rays, SBR）的技术。这种技术模拟射线在物体表面多次反射，同样基于几何光学。

我个人认为，把光理解为一种平面正弦波其实有些局限。例如3b1b的光学系列视频就把光的电场看成是一个平面正弦波。

虽然可以解释绝大多数现象，并且非常适合入门学习，但是对于波动光学渲染的研究而言，把光的电场简单描述为平面正弦波没办法进一步解释例如本篇paper提及的高斯光束。

但是我依旧强烈推荐不了解波动光学的读者先看这个系列的视频。

视频通过右旋手性介质会使光“旋转”这个特殊波动光学现象，科普了一系列波动光学的有趣现象以及背后的原理。

非常引人深思。

最后甚至还讲到了如何利用物质波构建全息影像。

• This tests your understanding of light | Optics puzzles 1
• How wiggling charges give rise to light | Optics puzzles 2
• But why would light “slow down”? | Optics puzzles 3
• How well do you understand refraction? | Optics puzzles 4
• How are holograms possible? | Optics puzzles 5

在经典电磁场理论中，我们通常使用平面波展开电磁场，每个模式的光子在空间上是非定域的，具有无限的空间延展。这种展开方式下，一个光子处于一个频率和波矢明确的模式上，因此它在空间中“充满”了整个平面波的区域。对于自由空间中的电磁场模式非常常见，但这种平面波并不适合描述局域性。

另一种展开方式是局域的波包（wave packet）。这种理解方式，允许我们在空间中对电磁场进行位置上的排布，即形成一系列具有一定宽度的波包。

把光单纯理解为一个正弦波其实违反了不确定性原理。如果光被视为正弦波，也就说明这个光是完全的单色光，但是现实中不可能有单色光。光信号作用基本是频域作用，空域到频域要用傅立叶变换。带宽-时间不确定性原理指出，带宽越窄，信号在时间上会越长。相反，带宽越宽，信号的时间长度越短。

在经典电磁理论中，一个波包可以对应于一个能量集中的区域。

但是注意光子也不能单纯理解为一种波包，而是理解为一种概率波。量子力学中的概率波函数描述了光子的存在几率。波包越集中，粒子性就越明显。这就是量子电动力学解释的波粒二象性。一个光子并不一定只能在一个波包中，它也可以描述为多个波包的叠加。因为光子状态本质上是量子场的激发，允许在不同位置的波包上进行叠加。

也就是说，一个光子可以“跨越”多个波包，即它的波函数可以在空间上以多个波包的形式存在，而不局限于某个特定位置。当一个波包模式上有一个光子时，这个波包可以看作是最低激发态；而如果波包上有多个光子（高阶激发态），则会体现出更高的能量。

本文中提到的一种电磁波束高斯光束是一种特定的电磁波解，可以视为一种波包。高斯光束是一个具有稳定振幅和相位结构的单一波包，它在横截面上表现为高斯分布。不同于我们平常在经典电磁理论中常用的平面波解，因为高斯光束的波前是曲率变化的，不是无限延伸的平面。

波的系综（ensemble of waves）指的是多个波的集合，这些波可能具有不同的频率、相位或传播方向。如果考虑一个系统中多个独立波包（例如多个脉冲激光束）相互叠加，可以将这些波包理解为一个波的系综。换句话说，如果有多个不相关的波包，尤其是随机相位的波包叠加在一起，它们在统计上可构成波的系综。

当我们把电磁场展开为一系列波包，那么就可以将每个波包视为一个随机事件，他们的到达时间、相位都是随机变量。对于多个波包的集合，我们可以在不同时间、不同位置观察到这些波包的特征。在统计意义上，使用系综平均来分析光的能量分布和波动行为。

$$
\langle U(\vec{r}, t) \rangle = \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} U(\vec{r}, t) \, dt
$$

由于不是物理壬，微元写在积分式前我是不认可的（

上面说到，光可以视为多个波包的集合。而不同波包之间可能会存在相位差和时间差，因此需要引入互相关函数和互相干函数来描述两个不同位置的波包之间的相干性和相对相位，提供了对这些波包在时空上相似性的量化描述。如果说系综平均用于描述一个信号在统计意义下的平均行为，那么互相关函数是在时间平均的角度上，描述了不同位置和不同时间之间的相干关系和波动关联性。

具体而言，互相关函数描述了位置 $\vec{r}_1$ 和 $\vec{r}_2$ 之间的波动相干性。如果两个波相干，那么这个Gamma值就会比较大：

$$
\Gamma(\vec{r}_1, \vec{r}_2, \tau) = \langle U(\vec{r}_1, t) U^*(\vec{r}_2, t + \tau) \rangle
$$

交叉谱密度（CSD）矩阵 $W(\vec{r}_1, \vec{r}_2, \omega)$ 是互相关函数的傅里叶变换，表示在频率域中两个位置之间的相干性。

还记得我们之前在Xia2023那篇黑狗毛中讨论的自相关函数吗？互相关函数需要两个信号来描述，而自相关函数只有一个信号。另一个角度来理解，自相关函数（ACF）是互相关函数的一种特殊情况。自相关函数描述的是信号与自身在不同时刻或不同位置的相关性，而互相关函数则是两个不同信号或同一信号在不同位置的相关性。互相关函数与交叉谱密度为一对傅里叶变换对，自相关函数与功率谱密度为一对傅里叶变换对。

交叉谱密度（CSD）和互相干函数描述在不同位置间波动的相干性。而辐射互谱密度（Radiance Cross-Spectral Density, RCSD）是CSD的推广，描述的是光辐射（即能量密度）在不同位置和不同方向之间的相关性。可以理解为，RCSD在辐射度量的基础上提供了类似于CSD的相干性描述。公式中， $L(\vec{r}_1, \vec{r}_2, \omega)$ 表示在频率 $\omega$ 下，位置 $\vec{r}_1$ 和 $\vec{r}_2$ 之间的辐射强度相干性。

辐射互谱密度传输方程（SDTE）与经典的光传输方程（Light Transport Equation, LTE）相似，但是更加适合波动光学。LTE描述的是从点到点的光辐射度传输，而SDTE则使用RCSD函数来描述光辐射的传播，相当于将传输视为区域间的相干传输。

SDTE中的RCSD通过区域间的积分形式表达辐射的传播，可以理解为用RCSD矩阵和衍射算子代替了传统的反射和散射。并且注意，SDTE基于RCSD函数，而不是具体的光强值，这一点和LTE是有较大区别的。

在进一步讨论利用边界元方法（Boundary Element Method, BEM）和自适应积分法（Adaptive Integral Method, AIM）加速的全波参考模拟器之前，首先简单回顾下前文的内容。前文介绍了PO方法和SDTE/RCSD理论。这些方法用于不同的散射计算需求，但它们的基本理论和适用范围不同。本文将讨论一种通过BEM和AIM结合提供高精度的面散射模拟的方法。

用原作者的话来总结，Wave optics在基于物理渲染里面属于很新的分支。虽然波动光学现象在生活中随处可见，但是对画面的影响并不是很大。这个方向其实还有很多可以做的地方。

其他相关参考：

• [Tannor 2007] 3.2 Gauss波包的一般性质 – 笠道梓的文章 – 知乎
• 高斯光束，平面波，球面波三者间有什么关系? – 定制滤光片的回答 – 知乎
• 波动光学中的一个波包和量子光学中的一个光子有什么关系吗？ – Godfly的回答 – 知乎
• 不确定性关系即傅里叶带宽定理之证明 – Michael Lieman的文章 – 知乎

A Full-Wave Reference Simulator for Computing Surface Reflectance

Paper主页：https://blaire9989.github.io/assets/1_BEMsim3D/project.html

ACM主页：https://dl.acm.org/doi/10.1145/3592414

ACM Citations：

Yunchen Yu, Mengqi Xia, Bruce Walter, Eric Michielssen, and Steve Marschner. 2023. A Full-Wave Reference Simulator for Computing Surface Reflectance. ACM Trans. Graph. 42, 4, Article 109 (August 2023), 17 pages. https://doi.org/10.1145/3592414

 

演讲报告

与更常用的光线追踪技术相比，波动光模拟更费时间但也更精确。

比方说金属表面的微观划痕、毛发纤维结构，传统的光学模型渲染的，无法渲染出我们现实生活中观察到的五彩斑斓的色彩效果。

基于波动光学的渲染是一个难题，因为解决麦克斯韦方程组需要大量复杂的计算。现有的基于波的外观模型通常采用一些近似方法。

一种近似方法是标量场近似（scalar field approximation）。

在波散射问题（wave scattering problem）中，电场和磁场是不同的矢量场量，具有不同极化方向（polarizations）的光由指向不同方向的场量组成。

有些近似模型将这两个矢量场量替换为单一标量函数（single scalar function），因此可以计算光能的强度，但放弃了对极化的建模（modeling polarization）。

另一种近似是一阶近似（first-order approximation）。假设光线仅与模型结构的每个部分发生一次反射，忽略了多次反射。然而，许多情况下这些近似都不适用。

例如 Yu 等人与 Dr. Lawrence 团队的合作，Penn State University 制作了带有圆柱形横截面的表面，这些表面会引起多次光反射并产生结构色彩，使用近似模型无法很好地理解或预测这些现象。

作者希望通过计算双向反射分布函数（BRDF）来尽可能精确地表征表面散射。

现有模型都采用各种近似方法，比如基于光线、标量或一次近似的模型，在没有参考质量（reference quality）的BRDF的情况下，很难看出每种反射模型缺少了什么或适用于什么场景。

作者的解决方案是构建一个三维四波模拟器（3D 4-way simulation），用于计算具有明确微观几何结构的表面的BRDF。

声称不使用任何近似方法，为表面样本计算出参考质量的BRDF。

速度方面dddd。

作者接下来介绍他们的模拟如何工作以及它与BRDF的关系。

首先，在左边这幅图，Input一个表面样本（定义为高度场）以及一个入射方向（在投影半球上表示）。定义了一个向表面传播的入射场，并从表面计算出一个散射场。

中间这幅图，光束是入射场，背景中的散射场也显示在此图中。

输出是对应给定入射方向的BRDF模式，半球每个点代表一个出射方向，颜色代表相应方向上反射光的颜色。BRDF以RGB颜色表示，这些颜色是从光谱数据转换而来。

对于很多粗糙度不高的表面，反射模式围绕镜面方向对称并随入射光方向的移动而变化。

接下来讲的是如何使用边界元法来解决仅在表面上的信号散射问题，从而降低了问题的维度。

表面信号是从麦克斯韦方程解出的表面电流，在离散化后，将问题构建为一个线性系统，并求解出表面电流和散射场。

为了使计算可行，作者将线性系统对称化，并使用一个适合对称矩阵的最小残差求解器。

此外，使用自适应积分方法加速矩阵向量乘法，这是一种基于快速傅里叶变换的加速方法，最初用于雷达计算。

代码大部分使用了Cuda C++包进行加速。

接下来，展示了一些结果，说明其计算的BRDF与之前方法得出的BRDF的比较。

[Yan 2018] 使用标量场近似的BRDF模型，这些模型只考虑一次折射。

[Xia 2023] 这篇使用矢量场量但也只考虑一次折射。

最精确的方法还得是咱们的。不仅使用矢量场量（vector field quantities），而且考虑了所有次序的反射（reflections of all orders）。

上图每个入射方向对应五种BRDF，分别代表不同的计算方法。

表面相对平滑的材料，表面上覆盖了一些各向同性的凸起（a bunch of isotropic bombs）。

第一行显示的是法线入射（normal instance），反射模式基本上居中（reflection pattern is pretty much centered）。

第二行由于入射光方向从某个倾斜方向来，模式向左偏移。

由于表面不太粗糙，五种方法的结果非常相似。

另一种材质呢，有一些棱棱角角（corner cubes）。每个corner cubes的三个面会让光多次反射，使反射光沿入射方向返回。叫做逆向反射（retroreflection）。

咱们的模拟器也可以模拟这种现象。左边四种方法都败下阵来。

原因在于如果只考虑一次反射，当光线击中corner cubes的一个面后，会被预测为向下进入下半球。

最后的例子是一个表面覆盖了一些球形坑。

由于坑边缘的高坡度（high slopes of the surface at the edges of the pits），导致多次反射的出现。

不同的方法看到明显差异。

可以看到预测的额外反射峰（extra lobe predicting）。（中间偏右那个部分）

另外整体更亮了。

这些差异源于不同次序反射间的干涉效果（interference between reflections of different orders）。

最后，简单介绍高效计算非常多的密集采样方向的BRDF的技术。

如果需要模拟的表面很大，速度就会很慢，并需要大量的GPU内存。

但是计算是线性的，可以将大面积表面分解为多个较小的子区域。

每个子区域上投射入射场，首先执行较小的子区域模拟，然后将散射场整合，得出整个大面积表面的BRDF。

对不同子区域的散射场应用不同的复数值缩放因子（different complex value scale factors），就可以综合出对应不同入射方向的大面积表面的BRDF。

这是因为对每个子区域的局部入射场，应用适当的相移会在表面上产生具有不同净方向的总入射场。在这张图中，入射场垂直传播。如果将五个焦点在不同空间位置的相同场叠加，得到一个空间上更宽的场，仍沿垂直方向传播。如果线性组合这些场，并对每个子区域的场应用适当的复数值缩放因子，整体场可以以略微倾斜的方向传播。

这里解释一下为啥复数缩放因子（different complex value scale factors）可以产生不同的入射方向。

• 这个因子可以调整波的幅度和相位。就比如往水面丢石头，两颗石头同时丢，水面的波就会更强。如果一颗稍微晚一些丢，波纹就可能会相互抵消（相消干涉）。这个因子就是控制石头投掷的时间。通过调整相位控制波的叠加方式，进而“引导”波向不同的方向传播。详细的去搜索「Beamforming」，雷达、无线通信和声纳等领域应用很广泛。

这三张图代表光波的波前。即波峰。可以理解为光在传播中的波形截面。

上方的图，垂直地射向表面，入射场分布集中在中心。光场集中在中心，沿垂直方向传播（即中间的黄线方向）。

左下方，多个相同入射场的叠加效果，但叠加时相位保持一致（即没有相位差）。多个入射场相加，使得整个场在空间上分布更宽了，但传播方向还是保持垂直。

右下方，多个入射场的叠加效果，但是在叠加时加入了相位差。相当于“偏移”了入射场的方向。呈现出一个倾斜方向。

演示中，两个视频展示了BRDF模式的移动。（我懒得截GIF了）

最后是大佬比心合影。

论文

我的第一篇文章已经粗略介绍了这项工作的内容和成果，接下来直接进入理论推导(Section3-5)。

3 FULL-WAVE SIMULATION

整体方法从一个表面模型开始，表面由高度场及其材料属性（如复折射率）描述，并指定一个目标点。

为了计算给定入射方向的BRDF，定义了一个入射场，该场从特定方向传播至目标点。

这个入射场作为输入，通过表面散射模拟进行处理，从而求解出相应的散射电磁场。

在本节（FULL-WAVE SIMULATION）中，将重点介绍BEM在应用场景中的原理。下一节会讲解如何高效实现BEM算法，最后一节会介绍如何结合多个模拟结果，以合成在入射和散射方向上密集采样的BRDF。

开局先来一张符号表吓一吓你。

3.1 Boundary Element Method: The Basics

边界元法（BEM）主要解决单一频率的散射问题，即特定频率的电磁波（包括电场和磁场）如何在不同介质的边界上反射和散射。这里的边界将空间分成了两个均匀区域，两个区域的材料特性（入射场所处的介质参数）用( $\epsilon_1, \mu_1$ )和( $\epsilon_2, \mu_2$ )表示。其中， $\epsilon$ 代表介电常数（介电率），$\mu$ 代表磁导率（磁导系数）。

在这种方法中，我们处理的是复数值场量，这些场量既包含振幅信息，也包含相位信息（即波的传播状态）。为了简化公式，我们假设所有波都是“时间谐和”的——也就是波随时间按照特定的周期变化。在整个文中， $e^{j\omega t}$ 项被省略，以简化表述。

3.1.1 Maxwell’s Equations and Surface Currents

首先，麦克斯韦方程描述了电场（E）和磁场（H）是如何相互影响的，决定了光波如何在不同材料之间传播和散射。为了化简，这里用“时间谐和”的形式：

$$ \begin{align} \nabla \times \mathbf{E} &= -\mathbf{M} – j \omega \mu \mathbf{H} \\ \nabla \times \mathbf{H} &= \mathbf{J} + j \omega \epsilon \mathbf{E} \tag{1} \end{align} $$

等号左边描述电场和磁场在空间中的“旋转”程度。 $M$ 和 $J$ 是表面电流密度（假想的电流），分别表示磁流和电流的密度（electric and magnetic current densities）。这个公式可以理解为，电场在边界附近“旋转”时，产生磁流和磁场的变化；磁场的旋转也会产生电场和电流的变化。

边界元法的核心思想是：在边界上引入表面电流，用这些电流来间接描述场的分布，而不需要计算每个区域中的所有点。三维问题化简为边界上的二维问题。

3.1.2 Source-Field Relationships

表面上的假想电流（电磁波的“源”）是如何产生散射的电磁场（“场”）的呢？

如Fig2.所示，在区域 $R_1$ ，总场（入射场和散射场）分别表示为 $E_1$ 和 $H_1$ ；

$$ \begin{align} &\mathbf{E}_1 = \mathbf{E}_i + \mathbf{E}_s \\ &\mathbf{H}_1 = \mathbf{H}_i + \mathbf{H}_s \tag{2} \end{align} $$

在区域 $R_2$ ，总场表示为 $E_2$ 和 $H_2$ 。上方的散射场由上方的电磁流产生；下面的散射场由下方的电磁流产生。

在均匀介质中，麦克斯韦方程组可以写成积分的形式，描述电场和磁场的生成方式。

\begin{align} \mathbf{E}(\mathbf{r}) &= -j \omega \mu (\mathcal{L} \mathbf{J})(\mathbf{r}) – (\mathcal{K} \mathbf{M})(\mathbf{r}) \\ \mathbf{H}(\mathbf{r}) &= -j \omega \epsilon (\mathcal{L} \mathbf{M})(\mathbf{r}) + (\mathcal{K} \mathbf{J})(\mathbf{r}) \tag{3} \end{align}

等号左边分别表示电磁场在 $r$ 处的电磁场强度，即描述的是场的“作用效果”。 $\mathcal{L}$ 和 $\mathcal{K}$ 是积分算子，表示场如何从表面电流和磁化强度产生。这两个算子定义为：

$$ \begin{aligned} & (\mathcal{L} \mathbf{X})(\mathbf{r})=\left[1+\frac{1}{k^2} \nabla \nabla \cdot\right] \int_V G\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \mathbf{X}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} \\ & (\mathcal{K} \mathbf{X})(\mathbf{r})=\nabla \times \int_V G\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \mathbf{X}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} \end{aligned} \tag{4} $$

$G(r, r{\prime})$ 是用于标量亥姆霍兹方程的三维格林函数，定义为：

$$
G(\mathbf{r}, \mathbf{r}^{\prime}) = \frac{e^{-jkr}}{4 \pi r} \quad \text{where } r = |\mathbf{r} – \mathbf{r}^{\prime}|
\tag{5}
$$

这个函数把散射表面的源场转换为散射区域内电磁场的分布。

本论文这里和[Xia 2023]中的公式(11)其实是一样的，但本论文将格林函数隐含在算子中，且积分域更为广泛。本质上都是描述如何从电流密度 $\mathbf{J}$ 和磁流密度 $\mathbf{M}$ 生成散射电场 $E(r)$。

求解麦克斯韦方程时，通过格林函数来整合各处的源（如电流、电荷）对空间中电磁场的影响。假设电磁场以 $e^{j\omega t}$ 的形式随时间变化（单频），就可以得到类似于亥姆霍兹方程的形式：$(\nabla^2 + k^2) \mathbf{E} = -j \omega \mu \mathbf{J}$ ，实际上这是一种“频域”形式的麦克斯韦方程。引入格林函数建立一种源-场关系，即将电流 $\mathbf{J}$ 和电荷 $\rho$ 作为“源”来计算电磁场 $\mathbf{E}$ 和 $\mathbf{H}$ 的分布。那么格林函数是满足如下方程的：$(\nabla^2 + k^2) G(\mathbf{r}, \mathbf{r}{\prime}) = -\delta(\mathbf{r} – \mathbf{r}{\prime})$ ，$\delta$ 是狄拉克δ函数，表示“点源”在空间中产生的标准波形，描述的是在空间中由一个“点源”激发的波场。通过格林函数，可以表达任意电流分布 $\mathbf{J}$ 在空间中对场点 $\mathbf{r}$ 产生的影响了！接着，每个源点上的“电流”或“电荷”都通过格林函数扩散到整个空间，对每一个场点产生累积影响。$\mathbf{E}(\mathbf{r}) = \int_V G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{J}(\mathbf{r}{\prime}) d \mathbf{r}{\prime}$ 这公式就是电场表示为源电流的积分叠加。最后，结合麦克斯韦方程组的积分形式和格林函数的思想。比如电场的积分形式可以表示为：$\mathbf{E}(\mathbf{r}) = -j \omega \mu \int_V G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{J}(\mathbf{r}{\prime}) d \mathbf{r}{\prime} – \nabla \times \int_V G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{M}(\mathbf{r}{\prime}) d \mathbf{r}{\prime}$ ，通过卷积可以获得每个源点上的微小电流分布对场点的累积影响，这是通过格林函数的积分在空间中传播和叠加的结果。格林函数作为卷积核，将空间中各点的电流分布通过积分传播到目标场点，实现了空间中各源点电流对整个场的累积影响。

介质区域 $R_1$ 和 $R_2$ 内的电场就分别表示为公式(6)(7)，包含了两个算子。

在区域 $R_2$ 内：

$$ \begin{align} & \mathbf{E}_1(\mathbf{r}) = -j \omega \mu_1 (\mathcal{L}_1 \mathbf{J}_1)(\mathbf{r}) – (\mathcal{K}_1 \mathbf{M}_1)(\mathbf{r}) \\ & \mathbf{H}_1(\mathbf{r}) = -j \omega \epsilon_1 (\mathcal{L}_1 \mathbf{M}_1)(\mathbf{r}) + (\mathcal{K}_1 \mathbf{J}_1)(\mathbf{r}) \end{align} \tag{6} $$

在区域 $R_2$ 内：

\begin{align} \mathbf{E}_2(\mathbf{r}) &= -j \omega \mu_2 (\mathcal{L}_2 \mathbf{J}_2)(\mathbf{r}) – (\mathcal{K}_2 \mathbf{M}_2)(\mathbf{r}) \\ \mathbf{H}_2(\mathbf{r}) &= -j \omega \epsilon_2 (\mathcal{L}_2 \mathbf{M}_2)(\mathbf{r}) + (\mathcal{K}_2 \mathbf{J}_2)(\mathbf{r}) \tag{7} \end{align}

这样就得到了在不同介质区域中电场和磁场的表现形式。

总结一下，这一小节通过假想的表面电流产生散射电磁场，将麦克斯韦方程组转化为积分表达，格林函数将电流和电荷分布转化为电磁场的积分叠加，展示了源-场关系的具体实现方式。最后给出了区域 $R_1$ 和 $R_2$ 内电场和磁场的表达式，展示了不同介质参数对场的影响。

3.1.3 Boundary Conditions

当电磁波在两种不同介质的边界传播时，会发生反射和折射。此时波的能量不可能凭空消失，而是在界面上平滑过渡的。如果电场或磁场在边界上不连续，就会出现不符合实际的能量跃变（即能量突然消失或增加），这违反了能量守恒。

具体可以搜索：「Interface conditions for electromagnetic fields」。

因此，必须在介质边界上满足一定的边界条件，以确保电磁场的连续性和能量守恒。具体而言，当电磁波在两种不同介质的界面上传播时，电场和磁场的切向分量（tangential component）需要在边界上保持连续性：

$$
\begin{aligned}
& \mathbf{n} \times (\mathbf{E}_1 – \mathbf{E}_2) = 0 \\
& \mathbf{n} \times (\mathbf{H}_1 – \mathbf{H}_2) = 0
\end{aligned}
\tag{8}
$$

并且边界上的净电磁流密度为零，即边界两侧电磁流密度方向相反大小相等。

$$
\begin{aligned}
& \mathbf{J} = \mathbf{J}_1 = -\mathbf{J}_2 \\
& \mathbf{M} = \mathbf{M}_1 = -\mathbf{M}_2
\end{aligned}
\tag{9}
$$

这两个条件同时满足，才不会破坏物理守恒。

3.1.4 Integral Equations

结合上面公式 (6)、(7)、(8) 和 (9)，得到关于电场和磁场的积分方程，称之为 PMCHWT（Poggio-Miller-Chang-Harrington-Wu-Tsai）方程：

$$
\begin{aligned}
& {\left[j \omega \mu_1\left(\mathcal{L}1 \mathbf{J}\right)(\mathbf{r})+j \omega \mu_2\left(\mathcal{L}_2 \mathbf{J}\right)(\mathbf{r})+\left(\mathcal{K}_1 \mathbf{M}\right)(\mathbf{r})+\right.}\left.\left(\mathcal{K}_2 \mathbf{M}\right)(\mathbf{r})\right]{\tan } \\
&=\left[\mathbf{E}i(\mathbf{r})\right]{\tan } \\
& {\left[\left(\mathcal{K}1 \mathbf{J}\right)(\mathbf{r})+\left(\mathcal{K}_2 \mathbf{J}\right)(\mathbf{r})-j \omega \varepsilon_1\left(\mathcal{L}_1 \mathbf{M}\right)(\mathbf{r})-j \omega \varepsilon_2\left(\mathcal{L}_2 \mathbf{M}\right)(\mathbf{r})\right]{\tan } } \\
&=-\left[\mathbf{H}i(\mathbf{r})\right]{\tan }
\end{aligned}
\tag{10}
$$

这两个方程分别还有个名字：

• EFIE（Electric Field Integral Equation），电场积分方程。 $j \omega \mu_1 (\mathcal{L}_1 \mathbf{J})(\mathbf{r})$ 和 $j \omega \mu_2 (\mathcal{L}_2 \mathbf{J})(\mathbf{r})$ 表示在介质 1 和 介质 2 中，电流密度 $\mathbf{J}$ 对电场的贡献。 $\mathcal{K}_1 \mathbf{M}$ 和 $\mathcal{K}_2 \mathbf{M}$ 表示在介质 1 和介质 2 中，磁流密度 $\mathbf{M}$ 对电场的贡献。
• MFIE（Magnetic Field Integral Equation），磁场积分方程。

总的来说，这是一种边界积分方程，专门用来求介电物体引起的电磁散射问题。有了这个PMCHWT方程，就可以准确计算电磁场的分布了。

3.1.5 Solving for Current Densities

这一节，需要通过上面的PMCHWT方程计算物体表面上的“电流”和“磁流”分布，这些分布决定了电磁波碰到物体后会怎么“反射”或者“折射”。

求解表面电流密度 $\mathbf{J}$ 和磁流密度 $\mathbf{M}$ 是通过将边界元离散化。对于离散单元定义一个基函数 $f_m(\mathbf{r})$ ，用基函数展开法表示电流密度和磁流密度分布。

$$
\mathbf{J}(\mathbf{r}) = \sum_{m=1}^{N} I_{J_m} f_m(\mathbf{r}); \quad \mathbf{M}(\mathbf{r}) = \sum_{n=1}^{N} I_{M_n} f_n(\mathbf{r})
\tag{11}
$$

N 是基函数的总数；$ I_{J_m}$ 和 $I_{M_n}$ 是对应基函数的未知系数，代表了每个单元上的电流和磁流强度。

通过这种基函数展开，连续的表面电流密度和磁流密度被分解成一系列基函数的线性组合。

为了求解对应基函数的未知系数，将电场积分方程 (EFIE) 和磁场积分方程 (MFIE) 转化为一个线性方程组。具体使用伽辽金法（Galerkin Method）完成。这个方法基本思想是，将积分方程作用在每个基函数上，加权平均使得积分方程在每个基函数的投影方向上都成立。该方法简单的说就是离散化、找基底、算系数。一个高维的线性方程组可以用线性代数方法简化。

这样可以将原本连续形式的PMCHWT方程的EFIE部分转换为有限个线性方程，把问题转化为解如下矩阵方程。

$$
\begin{bmatrix} A_{EJ} & A_{EM} \ A_{HJ} & A_{HM} \end{bmatrix} \begin{bmatrix} I_J \ I_M \end{bmatrix} = \begin{bmatrix} V_E \ V_H \end{bmatrix}
\tag{12}
$$

其中

$$
A_{\mathrm{EJ}}^{m n} =\int_S \mathbf{f}_m(\mathbf{r}) \cdot\left[j \omega \mu_1\left(\mathcal{L}_1 \mathbf{f}_n\right)(\mathbf{r})+j \omega \mu_2\left(\mathcal{L}_2 \mathbf{f}_n\right)(\mathbf{r})\right] d \mathbf{r} \tag{13}
$$

$$
A_{\mathrm{EM}}^{m n} =\int_S \mathbf{f}_m(\mathbf{r}) \cdot\left[\left(\mathcal{K}_1 \mathbf{f}_n\right)(\mathbf{r})+\left(\mathcal{K}_2 \mathbf{f}_n\right)(\mathbf{r})\right] d \mathbf{r} \tag{14}
$$

$$
A_{\mathrm{HJ}}^{m n} =\int_S \mathbf{f}_m(\mathbf{r}) \cdot\left[\left(\mathcal{K}_1 \mathbf{f}_n\right)(\mathbf{r})+\left(\mathcal{K}_2 \mathbf{f}_n\right)(\mathbf{r})\right] d \mathbf{r} \tag{15}
$$

$$
A_{\mathrm{HM}}^{m n} =-\int_S \mathbf{f}_m(\mathbf{r}) \cdot\left[j \omega \varepsilon_1\left(\mathcal{L}_1 \mathbf{f}_n\right)(\mathbf{r})+j \omega \varepsilon_2\left(\mathcal{L}_2 \mathbf{f}_n\right)(\mathbf{r})\right] d \mathbf{r} \tag{16}
$$

和

$$
V_{\mathrm{E}}^m =\int_S \mathbf{f}_m(\mathbf{r}) \cdot \mathbf{E}_i(\mathbf{r}) d \mathbf{r}\tag{17}
$$

$$
V_{\mathrm{H}}^m =-\int_S \mathbf{f}_m(\mathbf{r}) \cdot \mathbf{H}_i(\mathbf{r}) d \mathbf{r}\tag{18}
$$

公式(12)中，需要求出 $I_J$ 和 $I_M$ 。

不严谨地讲，公式(13)-(16)分别表示每一小块的电流密度对产生电场的贡献、每个小块的磁流密度对产生电场的贡献、每个小块的电流密度对产生磁场的贡献和每个小块的磁流密度对产生磁场的贡献。公式(17)(18)分别表示外界入射电场对这个小块电流的“推动力”和外界入射磁场对这个小块磁流的“推动力”。强调一下矩阵里面的元素，比如 $A_{EJ}^{mn}$ ，这其实是一个二重积分。由于对源点 $\mathbf{r}{\prime}$ 的积分已经在 $\mathcal{L}_1$ 和 $\mathcal{L}_2$ 中完成了，因此导致原paper看起来是个一重积分。

虽然原论文没有写，但是建议聪明的读者自己推导一下。尝试根据上文提到的公式(4)把公式(13)展开。我这里尝试推导了一下，有错误请指正。首先把两个算子代入，注意这里梯度算子的位置：

$$
\begin{aligned}
A_{\mathrm{EJ}}^{mn} &= j \omega \mu_1 \int_S \mathbf{f}_m(\mathbf{r}) \cdot \left\{ \int_V G_1(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime} + \frac{1}{k_1^2} \nabla \left[ \nabla \cdot \int_V G_1(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime} \right] \right\} d\mathbf{r} \\
&\quad + j \omega \mu_2 \int_S \mathbf{f}_m(\mathbf{r}) \cdot \left\{ \int_V G_2(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime} + \frac{1}{k_2^2} \nabla \left[ \nabla \cdot \int_V G_2(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime} \right] \right\} d\mathbf{r}
\end{aligned}
$$

先考虑其中一个梯度项：

$$
\int_S \mathbf{f}_m(\mathbf{r}) \cdot \nabla \left[ \nabla \cdot \int_V G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime} \right] d\mathbf{r}
$$

用向量分布积分展开：

$$
\int_S \mathbf{f}m \cdot \nabla B \, d r = -\int_S B (\nabla \cdot \mathbf{f}_m) \, dr + \int{\partial S} B (\mathbf{A} \cdot \mathbf{n}) \, dr
$$

其中，散度项 $B$ ：

$$
B = \nabla \cdot \int_V G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime}
$$

但是很抱歉，这里是物理。在边界条件下，边界项直接简化为零，得到：

$$
\int_S \mathbf{f}_m \cdot \nabla B \, dS = -\int_S B (\nabla \cdot \mathbf{f}_m) \, dS
$$

对于散度项 $B$ ，可以直接展开散度：

$$
B = \nabla \cdot \int_V G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime} = \int_V (\nabla \cdot G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{f}_n(\mathbf{r}{\prime})) d\mathbf{r}{\prime}
$$

一个是标量函数，一个是矩阵函数，因此根据散度的乘积法则：

$$
\nabla \cdot (G \mathbf{f}_n) = (\nabla G) \cdot \mathbf{f}_n + G (\nabla \cdot \mathbf{f}_n)
$$

但是很抱歉，这里是物理。由于基函数满足无散度条件，因此这里直接简化：

$$
\nabla \cdot (G \mathbf{f}_n) = (\nabla G) \cdot \mathbf{f}_n
$$

同时我们注意到格林函数的对称性 $\nabla G(\mathbf{r}, \mathbf{r}{\prime}) = -\nabla{\prime} G(\mathbf{r}, \mathbf{r}{\prime})$ ，于是有：

$$
B = \nabla \cdot \int_V G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime} = -\int_V \nabla{\prime} G(\mathbf{r}, \mathbf{r}{\prime}) \cdot \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime}
$$

对这一项也进行分布积分：

$$
\int_V \nabla{\prime} G(\mathbf{r}, \mathbf{r}{\prime}) \cdot \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime} = \int_V \nabla{\prime} \cdot (G \mathbf{f}_n) d\mathbf{r}{\prime} – \int_V G (\nabla{\prime} \cdot \mathbf{f}_n) d\mathbf{r}{\prime}
$$

但是很抱歉，这里是物理。边界项再次化简：

$$
\int_V \nabla{\prime} G(\mathbf{r}, \mathbf{r}{\prime}) \cdot \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime} = -\int_V G(\mathbf{r}, \mathbf{r}{\prime}) (\nabla{\prime} \cdot \mathbf{f}_n(\mathbf{r}{\prime})) d\mathbf{r}{\prime}
$$

最终得到：

$$
B = \nabla \cdot \int_V G(\mathbf{r}, \mathbf{r}{\prime}) \mathbf{f}_n(\mathbf{r}{\prime}) d\mathbf{r}{\prime} = \int_V G(\mathbf{r}, \mathbf{r}{\prime}) (\nabla{\prime} \cdot \mathbf{f}_n(\mathbf{r}{\prime})) d\mathbf{r}{\prime}
$$

波数 $k_i$ 于介质参数的关系：

$$
k_i^2 = \omega^2 \mu_i \varepsilon_i \quad \Rightarrow \quad \frac{1}{k_i^2} = \frac{1}{\omega^2 \mu_i \varepsilon_i}
$$

另一个梯度项也是同理，最后得到最终的表达式 $A_{\mathrm{EJ}}^{mn}$ ：

$$
\begin{aligned} A_{\mathrm{EJ}}^{mn} &= j \omega \mu_1 \int_S \int_{V_1} \mathbf{f}m(\mathbf{r}) \cdot \mathbf{f}n(\mathbf{r}{\prime}) G_1(\mathbf{r}, \mathbf{r}{\prime}) d\mathbf{r}{\prime} d\mathbf{r} \\ &\quad – \frac{j}{\omega \varepsilon_1} \int_S \int{V_1} (\nabla \cdot \mathbf{f}m(\mathbf{r})) G_1(\mathbf{r}, \mathbf{r}{\prime}) (\nabla{\prime} \cdot \mathbf{f}n(\mathbf{r}{\prime})) d\mathbf{r}{\prime} d\mathbf{r} \\ &\quad + j \omega \mu_2 \int_S \int{V_2} \mathbf{f}m(\mathbf{r}) \cdot \mathbf{f}n(\mathbf{r}{\prime}) G_2(\mathbf{r}, \mathbf{r}{\prime}) d\mathbf{r}{\prime} d\mathbf{r} \\ &\quad – \frac{j}{\omega \varepsilon_2} \int_S \int{V_2} (\nabla \cdot \mathbf{f}_m(\mathbf{r})) G_2(\mathbf{r}, \mathbf{r}{\prime}) (\nabla{\prime} \cdot \mathbf{f}_n(\mathbf{r}{\prime})) d\mathbf{r}{\prime} d\mathbf{r} \end{aligned}
$$

同样的操作，得到剩下的矩阵元素，这里直接抄论文附加材料的内容了：

$$
\begin{aligned} A_{\mathrm{EM}}^{m n}= & A_{\mathrm{HJ}}^{m n}=\int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \mathbf{f}_m(\mathbf{r}) \cdot\left[\nabla G_1\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \times \mathbf{f}_n\left(\mathbf{r}^{\prime}\right)\right] d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \mathbf{f}_m(\mathbf{r}) \cdot\left[\nabla G_2\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \times \mathbf{f}_n\left(\mathbf{r}^{\prime}\right)\right] d \mathbf{r}^{\prime} d \mathbf{r} \\ A_{\mathrm{HM}}^{m n} & =-j \omega \varepsilon_1 \int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \mathbf{f}_m(\mathbf{r}) \cdot \mathbf{f}_n\left(\mathbf{r}^{\prime}\right) G_1\left(\mathbf{r}, \mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\frac{j}{\omega \mu_1} \int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \nabla \cdot \mathbf{f}_m(\mathbf{r}) G_1\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \nabla^{\prime} \cdot \mathbf{f}_n\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -j \omega \varepsilon_2 \int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \mathbf{f}_m(\mathbf{r}) \cdot \mathbf{f}_n\left(\mathbf{r}^{\prime}\right) G_2\left(\mathbf{r}, \mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\frac{j}{\omega \mu_2} \int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \nabla \cdot \mathbf{f}_m(\mathbf{r}) G_2\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \nabla^{\prime} \cdot \mathbf{f}_n\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r}\end{aligned}
$$

得到矩阵的每一个元素之后，引入平移不变函数（Shift-invariant Functions）。帮助得到在不同坐标系下表示格林函数及其梯度。

$$
\begin{aligned}
& g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right)=G_i\left(\mathbf{r}, \mathbf{r}^{\prime}\right)=\frac{e^{-j k_i r}}{4 \pi r} \\
& g_{2, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right)=\hat{\mathbf{x}} \cdot \nabla G_i\left(\mathbf{r}, \mathbf{r}^{\prime}\right)=-\left(x-x^{\prime}\right)\left(\frac{1+j k_i r}{4 \pi r^3}\right) e^{-j k_i r} \\
& g_{3, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right)=\hat{\mathbf{y}} \cdot \nabla G_i\left(\mathbf{r}, \mathbf{r}^{\prime}\right)=-\left(y-y^{\prime}\right)\left(\frac{1+j k_i r}{4 \pi r^3}\right) e^{-j k_i r} \\
& g_{4, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right)=\hat{\mathbf{z}} \cdot \nabla G_i\left(\mathbf{r}, \mathbf{r}^{\prime}\right)=-\left(z-z^{\prime}\right)\left(\frac{1+j k_i r}{4 \pi r^3}\right) e^{-j k_i r} \quad \text { where } r=\left|\mathbf{r}-\mathbf{r}^{\prime}\right|
\end{aligned}
$$

然后将矩阵元素展开为基函数的不同分量的组合，分量作用在平移不变函数上。最终矩阵的所有元素，都有如下形式。这样的形式可以加速边界元矩阵的构造和求解。

$$
\int_{\mathbf{f}m} \int{\mathbf{f}_n} \psi_m(\mathbf{r}) g\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \xi_n\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r}
$$

最终得到：

Here, $x, y, z, x^{\prime}, y^{\prime}, z^{\prime}$ are the Cartesian components of $\mathbf{r}, \mathbf{r}^{\prime}$. Now we have for $i=1,2$ :

$$ \begin{aligned} A_{\mathrm{EJ}, i}^{m n} &= j \omega \mu_i \int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \mathbf{f}_{m x}(\mathbf{r}) \, g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \, \mathbf{f}_{n x}\left(\mathbf{r}^{\prime}\right) \, d \mathbf{r}^{\prime} \, d \mathbf{r} \\ & \quad + j \omega \mu_i \int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \mathbf{f}_{m y}(\mathbf{r}) \, g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \, \mathbf{f}_{n y}\left(\mathbf{r}^{\prime}\right) \, d \mathbf{r}^{\prime} \, d \mathbf{r} \\ & \quad + j \omega \mu_i \int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \mathbf{f}_{m z}(\mathbf{r}) \, g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \, \mathbf{f}_{n z}\left(\mathbf{r}^{\prime}\right) \, d \mathbf{r}^{\prime} \, d \mathbf{r} \\ & \quad – \frac{j}{\omega \varepsilon_i} \int_{\mathbf{f}_m} \int_{\mathbf{f}_n} \nabla \cdot \mathbf{f}_m(\mathbf{r}) \, g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \, \nabla^{\prime} \cdot \mathbf{f}_n\left(\mathbf{r}^{\prime}\right) \, d \mathbf{r}^{\prime} \, d \mathbf{r} \end{aligned} \tag{S.18} $$

where $\mathbf{f}{m x}, \mathbf{f}{m y}, \mathbf{f}_{m z}$ are the $x, y, z$ components of the vector basis function $\mathbf{f}_m$ . Similarly, we have:

$$ \begin{aligned} & A_{\mathrm{EM}, i}^{m n}=A_{\mathrm{HJ}, i}^{m n}=\int_{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m z}(\mathbf{r}) g_{2, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n y}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m y}(\mathbf{r}) g_{2, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n z}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m x}(\mathbf{r}) g_{3, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n z}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m z}(\mathbf{r}) g_{3, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n x}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m y}(\mathbf{r}) g_{4, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n x}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -\int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m x}(\mathbf{r}) g_{4, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}_{n y}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \end{aligned} \tag{S.19} $$

Lastly, we also have:

$$ \begin{aligned} A_{\mathrm{HM}, i}^{m n} & =-j \omega \varepsilon_i \int_{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m x}(\mathbf{r}) g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n x}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -j \omega \varepsilon_i \int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m y}(\mathbf{r}) g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n y}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & -j \omega \varepsilon_i \int{\mathbf{f}m} \int{\mathbf{f}n} \mathbf{f}{m z}(\mathbf{r}) g_{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \mathbf{f}{n z}\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \\ & +\frac{j}{\omega \mu_i} \int{\mathbf{f}m} \int{\mathbf{f}n} \nabla \cdot \mathbf{f}_m(\mathbf{r}) g{1, i}\left(\mathbf{r}-\mathbf{r}^{\prime}\right) \nabla^{\prime} \cdot \mathbf{f}_n\left(\mathbf{r}^{\prime}\right) d \mathbf{r}^{\prime} d \mathbf{r} \end{aligned} \tag{S.20} $$

这部分代码在MVProd类中，看不懂也没问题，因为上面都是我瞎写瞎抄的，已经超出图形学的研究范畴了。

另外作者还讨论了EFIE和MFIE的对称性。正是这种对称性使得计算效率和空间利用率更低。注意到：

$$
A_{EJ} = A_{EJ}^T, \quad A_{HM} = A_{HM}^T \tag{19}
$$

$$
A_{EM} = A_{EM}^T, \quad A_{HJ} = A_{HJ}^T, \quad A_{EM} = A_{HJ}\tag{20}
$$

由于矩阵是对称的，我们不需要计算所有的矩阵元素，也不用存储所有矩阵元素。

在解出表面电流密度后，可以再使用公式(6)来计算从散射表面向外传播的散射场。

3.2 Rough Surface Scattering: The Specifics

在模拟电磁波与粗糙表面之间的相互作用时，表面的不规则几何结构会对散射特性造成显著影响。作者对粗糙表面离散化，将连续的表面划分为多个单元，每个单元都可以应用PMCHWT方程做数值计算的方法。

3.2.1 Rough Surface Samples

将粗糙表面表示为一个二维的高度场（height field），然后将该高度场进行离散化处理，划分为多个矩形单元。

每次模拟只考虑尺寸为 $$L_x \times L_y$$ 的表面样本。并且选择一个步进 $d$ 来定义离散化的网格。

$$
x_s = s \cdot d, \quad s = 0, 1, \ldots, N_x
\
y_t = t \cdot d, \quad t = 0, 1, \ldots, N_y
\tag{21}
$$

其中，$ N_x = L_x / d$ 和 $N_y = L_y / d$ ，分别表示在 $x$ 和 $y$ 方向上被划分成的单元数。在每个离散点 $(x_s, y_t)$ 都有一个高度场 $h(x_s, y_t)$ 函数。

作者总结了，粗糙表面的高度变化尺度非常小，通常只有几微米。也就是和可见光电磁波波长相当。

3.2.2 Basis Elements and Functions

每个基元有四个角，每个角有不同的高度，并且每个角对电流磁流贡献都不同。因此在每个小方块上，都定义了四个基函数，近似表示每个小方块上的情况。

在大多数模拟中，步长 $d$ 选取为波长的 $\lambda / 16$ 左右，以确保精度。

每个基元由两个参数 u 和 v 参数化，范围均在 [-1, 1]。

基元的形状通过一个双线性函数 $\mathbf{r}(u, v)$ 来表示，其中 $(s, t)$ 表示当前基元的索引，且基元由四个顶点的坐标决定：

$$ \begin{aligned} \mathbf{r}(u, v) &= \frac{(1 – u)(1 – v)}{4} \mathbf{p}_{s-1, t-1} +\\ & \frac{(1 – u)(1 + v)}{4} \mathbf{p}_{s-1, t} + \\ \frac{(1 + u)(1 – v)}{4} \mathbf{p}_{s, t-1} + \frac{(1 + u)(1 + v)}{4} \mathbf{p}_{s, t} \end{aligned} \tag{22} $$

其中， $\mathbf{p}{s, t} = (x_s, y_t, z{s, t})$ 是基元的四个顶点坐标。

在每个矩形基元上定义四个基函数 $f_1(u, v), f_2(u, v), f_3(u, v), f_4(u, v)$ ，它们的形式为：

$$
\begin{aligned}
& f_1(u, v) = \frac{(1 – u)}{J(u, v)} \frac{\partial \mathbf{r}(u, v)}{\partial u}, \quad f_2(u, v) = \frac{(1 + u)}{J(u, v)} \frac{\partial \mathbf{r}(u, v)}{\partial u} \\
& f_3(u, v) = \frac{(1 – v)}{J(u, v)} \frac{\partial \mathbf{r}(u, v)}{\partial v}, \quad f_4(u, v) = \frac{(1 + v)}{J(u, v)} \frac{\partial \mathbf{r}(u, v)}{\partial v}
\end{aligned}
\tag{23}
$$

这里的 Jacobian $J(u, v)$ 表示如下：

$$
J(u, v) = \left| \frac{\partial \mathbf{r}(u, v)}{\partial u} \times \frac{\partial \mathbf{r}(u, v)}{\partial v} \right|
\tag{24}
$$

Jacobian 的引入用于转换坐标系并确保基函数在不同的 $u, v$ 方向上具有合适的比例关系。

3.2.3 Gaussian Beam Incidence

由于不太了解光学，以下为个人理解。高斯光束是激光发射出去的时候，在行波场中间部分出现往内部凹陷的一种现象，换一句话说高斯光束是描述光在横截面上的能量分布。而平面波和球面波的重点是描述能量传播的方向。在传播过程中，高斯光束的波前形状近似为球面波。

Gouy 相位（Gouy phase）是高斯光束传播中的一种相位延迟效应，光束通过焦点后相位会额外增加。另一方面，高斯光束满足麦克斯韦方程组在傍轴条件下的一个解，可以近似为非均匀的球面波。感觉目前也不用太深入研究。

有关高斯光束的资料：

• https://en.wikipedia.org/wiki/Gaussian_beam
• 激光原理(Principle of Laser)笔记 – Senner的文章 – 知乎
• 高斯光束：基本公式和应用 – 中国光学
• 高斯光束：基本公式和应用 – Andy的文章 – 知乎
• 光子的轨道角动量 – Godfly的文章 – 知乎

回到原论文，高斯光束好处在于可以控制入射场的大小，进而能够控制一个稍微比照射区域大的区域的表面感应电流密度是非零数值。

高斯光束是一种电磁波，其振幅在垂直于传播方向的平面上呈二维高斯分布[Paschotta 2008]，它的能量主要集中在光束中心附近。看上方图(a)，高斯光束可以用焦平面 $P$ 、中心点 $o$ 、和光束腰径（beam waist） $w$ 来描述。场强随位置的衰减关系为 $e^{-r^2 / w^2}$ ，当距离中心超过 $2.5w$ 时，场强衰减至极小，可认为几乎为零。

虽然但是，高斯光束也具有一定的发散性。发散角 $\theta$ 近似与波长 $\lambda$ 成正比，与光束腰径 $w$ 成反比。公式：

$$
\theta = \frac{\lambda}{\pi \eta w} \tag{25}
$$

当光束斜着射入表面时，高斯光束在焦平面上有一个椭圆形的横截面。如下图所示。在两个垂直方向上有不同的腰斑尺寸（beam waist）：一个平行于入射方向和表面法线的平面，另一个与之垂直。

为了保证不同方向上的高斯光束在表面上的照射区域相同，这里引入了两个横向方向的束腰宽度：

$$
w_1 = w, \quad w_2 = w \cos \theta_i\tag{26}
$$

即使在不同入射角度下，表面上的照射面积保持一致。这对推导BRDF而言非常重要。

4. IMPLEMENTATION AND ACCELERATION

但是如果按上面的方法硬算，是不可能会有结果的。因此需要一些加速手段。

4.1 The Adaptive Integral Method

想要直接计算上面公式(12)的方程组，计算量是不可接受的。

$$ \begin{bmatrix} A_{EJ} & A_{EM} \ A_{HJ} & A_{HM} \end{bmatrix} \begin{bmatrix} I_J \ I_M \end{bmatrix} = \begin{bmatrix} V_E \ V_H \end{bmatrix} \tag{12} $$

按照原本的思路，用 $N$ 个基函数来表示电流磁流密度，矩阵的规模就是 $2N \times 2N$ 。如果直接求解矩阵（LU分解，Cholesky分解等），总复杂度可能在 $\mathcal{O}(N^3)$ 。就算用共轭梯度法总复杂度也在 $\mathcal{O}(N^2)$ 。在一些小规模的仿真中，基函数的规模大概在 $960*960$ ，存储需求大概就在 29.4GB 。用了Adaptive Integral Method, AIM，按8字节来算总存储需求约在 76.8 MB。究竟是什么方法这么神奇呢？小编接下来就带大家一起看看吧！

4.1.1 Approximating Matrix Elements

AIM最初是由Bleszynski等人[1996]提出的。AIM的核心思想是将每个基函数的作用近似为一组点源的作用，避免直接计算每一对基函数之间的精确相互作用，同时通过FFT将各基函数的影响传播开来，提升计算效率。

AIM中矩阵元素的计算方式是对矩阵元素的某些项的线性组合来近似计算。这就是为什么在上面公式(13)-(16)后面让读者们进一步推导，推导最终的结果使其符合AIM的形式。

$$
\int_{f_m} \int_{f_n} \psi_m(r) g(r – r{\prime}) \xi_n(r{\prime}) \, dr{\prime} \, dr
\tag{27}
$$

AIM首先会在包含电磁场和场源的空间内创建一个全局3D笛卡尔网格，如下图(6)所示。

为了进一步化简公式(27)，AIM算法会让原本的基函数近似为在这个三维笛卡尔坐标上一组网格点上的点源。说白了就是连续变离散，并且方便后续FFT。

$$ \psi_m(r) \approx \tilde{\psi}m(r) := \sum{p \in S_m} \Lambda_{mp} \delta^3(r – p) \\ \xi_n(r{\prime}) \approx \tilde{\xi}n(r{\prime}) := \sum{q \in S_n} \Lambda{\prime}_{nq} \delta^3(r{\prime} – q) \tag{28} $$

将公式(27)代入公式(28)，换句话说，就是将双重积分的形式转化为了一个双重求和的形式。

$$
\sum{p \in S_m} \sum_{q \in S_n} \Lambda_{mp} g(p – q) \Lambda{\prime}_{nq}
\tag{29}
$$

方法详细参考：

Kai Yang and Ali E Yilmaz. 2011. Comparison of precorrected FFT/adaptive integral method matching schemes. Microwave and Optical Technology Letters 53, 6 (2011), 1368–1372.

4.1.2 Base and Correction Matrices

基于公式(29)定义一组基础近似（base approximation）矩阵 $B_{EJ}, B_{EM}, B_{HJ}, B_{HM}$ 作为近似，专门处理距离较远的基函数对。这些矩阵通过引入 $\Lambda$ 矩阵和卷积等操作化简计算。同时，对于距离较近（$d_{near}$）的基函数对，再引入修正矩阵（Correction Matrices）来减少误差。 $C_{EJ}, C_{EM}, C_{HJ}, C_{HM}$ 是一种稀疏矩阵，定义如下：

$$ C_{\mathrm{X}}^{m n}=\left\{\begin{array}{ll} A_{\mathrm{X}}^{m n}-B_{\mathrm{X}}^{m n} & d_{m n} \leq d_{\text {near }} \\ 0 & \text { otherwise } \end{array} \quad \mathrm{X} \in\{\mathrm{EJ}, \mathrm{EM}, \mathrm{HJ}, \mathrm{HM}\}\right. \tag{30} $$

$A_X^{mn}$ 是原始矩阵的精确值，而 $B_X^{mn}$ 是基础矩阵的近似值。通过减去基础矩阵的近似值，得到一个较准确的修正项，用于补偿近距离基函数对的误差。

综上所述，AIM方法中每个矩阵的最终近似形式可以写成如下关系：

$$ \begin{aligned} A_{\mathrm{EJ}} \approx B_{\mathrm{EJ}}+C_{\mathrm{EJ}} ; & A_{\mathrm{EM}} \approx B_{\mathrm{EM}}+C_{\mathrm{EM}} ; \\ A_{\mathrm{HJ}} \approx B_{\mathrm{HJ}}+C_{\mathrm{HJ}} ; & A_{\mathrm{HM}} \approx B_{\mathrm{HM}}+C_{\mathrm{HM}} \end{aligned} $$

换句话说，原始矩阵可以通过基础矩阵和修正矩阵的组合来近似。

4.1.3 Fast Matrix-Vector Multiplication

快速矩阵-向量乘法（Fast Matrix-Vector Multiplication）是AIM的核心。

由于上面得到的修正矩阵 $C$ 是稀疏矩阵， $C$ 只在近距离的基函数上有非零值，因此矩阵 $C$ 的乘法操作是很快的。

利用基础近似矩阵 $B$ 的卷积特性，计算了矩阵 $B$ 与向量的乘积。计算过程分为三步：

$$
y_1 = \Lambda_2^T x, \quad y_2 = G y_1, \quad y_3 = \Lambda_1 y_2
\tag{32}
$$

第一步把向量投影到一个稀疏矩阵网格上。

第二步也是核心步骤，把网格点的数据传播到整个网络，也就是计算每个点对其他点的影响。两个点比较接近，矩阵 $G$ 中的传播函数就大。通过FFT来做加速。

$$
y_2 = \mathcal{F}^{-1} { \mathcal{F}(g) \mathcal{F}(y_1) }
\tag{33}
$$

第三步骤将结果映射回原来的基函数空间。

4.2 GPU-Accelerated Iterative Solving

在GPU上将AIM方法中的计算重点转移到快速傅里叶变换（FFT）和稀疏矩阵操作上。总结一下，目前将大型矩阵分为基础矩阵 $B$ 和修正矩阵 $C$ ，分别处理远距离和近距离的基函数对。

• cuFFT：将基础矩阵的乘法操作转换为频域中的卷积计算
• cuSPARSE：算稀疏矩阵加速修正矩阵 C

并且优化计算策略：

• 对于小规模仿真任务，只需要使用1个GPU
• 对于大规模仿真任务，将任务分配到4个GPU

4.2.1 Small-Scale Simulations

对于小规模仿真任务（例如 $12 \mu m \times 12 \mu m$）， 必须事先计算并存储传播函数的傅里叶变换值（即矩阵 $G$ 的傅里叶变换）。在小规模任务中稀疏修正矩阵 C 占用的显存不到5GB，一张GPU就可以搞掂。

4.2.2 Large-Scale Simulations

对于大规模仿真任务（例如 $24 \mu m \times 24 \mu m$），由于单个GPU的显存不足以存储所有数据，作者将计算任务分配到4个GPU上。在这个尺度的仿真上，基函数的个数会达到960 × 960个，存储所有修正矩阵的非零元素（包括行列索引和复数浮点数值）大约需要20GB显存。策略还是和小规模一样，每个GPU分配大约5GB内存来存储修正矩阵 $C$ 。

MINRES求解器在主机CPU上执行，而矩阵-向量乘积 $y = Ax$ 的计算在GPU上完成。但是不需要担心传输时间，这个向量大概只有30MB。

4.3 FFT-Accelerated Scattered Field Evaluation

用FFT加速计算散射场。在远场区域上评估从表面散射的场最终求出表面BRDF。

在求解BEM后，得到了表面的电流密度 $\mathbf{J}$ 和磁流密度 $\mathbf{M}$ 。这些密度分布定义了表面上的电磁源，可以用来计算在远场区域上的散射场。公式很简单，随着距离衰弱同时还会具有一定的相位变化：

$$
\mathbf{E_s}(r) \approx \mathbf{E}(\hat{r}) \frac{e^{-jkr}}{r}; \quad \mathbf{H_s}(r) \approx \mathbf{H}(\hat{r}) \frac{e^{-jkr}}{r}
\tag{36}
$$

公式右边的 $\mathbf{E}(\hat{r})$ 和 $\mathbf{H}(\hat{r})$ 是在远场上特定方向 $\hat{r}$ 的振幅。在不同的方向上，散射场的强度可能不同。

$$
\begin{aligned}
F_1(\hat{\mathbf{r}})=\int_V J_x\left(\mathbf{r}^{\prime}\right) e^{j k \mathbf{r}^{\prime} \cdot \hat{\mathbf{r}}} d \mathbf{r}^{\prime} ; & F_2(\hat{\mathbf{r}})=\int_V J_y\left(\mathbf{r}^{\prime}\right) e^{j k \mathbf{r}^{\prime} \cdot \hat{\mathbf{r}}} d \mathbf{r}^{\prime} \\
F_3(\hat{\mathbf{r}})=\int_V J_z\left(\mathbf{r}^{\prime}\right) e^{j k \mathbf{r}^{\prime} \cdot \hat{\mathbf{r}}} d \mathbf{r}^{\prime} ; & F_4(\hat{\mathbf{r}})=\int_V M_x\left(\mathbf{r}^{\prime}\right) e^{j k \mathbf{r}^{\prime} \cdot \hat{\mathbf{r}}} d \mathbf{r}^{\prime} \\
F_5(\hat{\mathbf{r}})=\int_V M_y\left(\mathbf{r}^{\prime}\right) e^{j k \mathbf{r}^{\prime} \cdot \hat{\mathbf{r}}} d \mathbf{r}^{\prime} ; & F_6(\hat{\mathbf{r}})=\int_V M_z\left(\mathbf{r}^{\prime}\right) e^{j k \mathbf{r}^{\prime} \cdot \hat{\mathbf{r}}} d \mathbf{r}^{\prime}
\end{aligned}
\tag{37}
$$

为了避免直接求解这些积分，作者利用先前（4.1.2）对 $\mathbf{J}$ 和 $\mathbf{M}$ 的点源近似（$\Lambda$ 矩阵），将每个积分项 $F_i(\hat{r})$ 离散化并重写为傅里叶变换的形式，如公式 (38) 所示：

$$
F_i(\hat{r}) = \sum_{p \in S} h_i(p) e^{jp \cdot k \hat{r}}
\tag{38}
$$

将连续的场强计算转化为离散的求和，这样就可以用FFT快速计算了。并且从上式可以观察到， $F_i(\hat{r})$ 实际上是 $h_i(p)$ 在空间频率 $-k\hat{r}$ 上的傅里叶分量。

The required spatial frequencies are not on the FFT grid but can be interpolated; we add zero padding prior to the FFT step, to ensure enough resolution in the frequency domain for the trilinear interpolation to be sufficiently accurate.

5 HIGH RESOLUTION BRDF GENERATION

搞了一大堆，终于回到熟悉的BRDF计算了。这里关键在于使用小尺度模拟的线性叠加来重建大尺度入射场的远场散射，而不是一口吃成大胖子。$N^2$个小尺度的高斯光束构成的网格来线性组合成近似大尺度的入射场。

这里用到“波束引导”（beam steering）技术。这种方法不用对每个方向都进行一次模拟，从而大幅降低计算成本。

5.1 Basic and Derived Incident Directions

首先，沿某方向 $\mathbf{u}$ 传播的 $N^2$ 个高斯光束组成在接收平面的一个 $N \times N$ 点的网格。这些光束组合后能够生成一个大的总场。

然后给每个高斯光束引入复数缩放因子，调整每个光束的相位，进而调整组合场的传播方向，这些方向称为desired direction。

$$
a_{st} = e^{j k \mathbf{p}{st} \cdot \omega_i} \tag{39}
$$

当目标入射方向 $-\omega_i$ 与基本方向 $\mathbf{u}$ 之间的夹角接近各高斯光束的发散角（divergence of the small beams）时，aliasing artifacts begin to appear. An example is shown in Fig. 8 (d).

In our framework, we decide on a primary waist w and choose a collection of basic incident directions. In general, 较小的腰宽意味着更大的发散角，这样可以从每个基本方向派生出更多的入射方向，从而减少所需的基本方向数量。较大的腰宽会降低每个高斯光束的发散角，使得组合后的总场发散更小，从而产生更精确的入射方向。

每个六边形的中心对应一个基本入射方向。整个半球的所有入射方向划分为若干territories，每个territories归属于一个基本入射方向。在半球投影中，这种比例关系与余弦因子相互抵消，因此可以使用大小相等的六边形来表示。

对于每个基本方向，可以通过“波束引导”产生一些偏离该基本方向的派生方向。

当光的入射角很小时（例如接近表面法线方向），基本入射方向附近的派生方向范围很小，因为小角度的光更集中，不会有很大的扩散。反之，派生方向范围会更大。总结就是，入射角越大（角度越接近水平），派生方向的覆盖范围就越大。公式作者也提到了：$1/\cos \theta_i$。

5.2 Individual Simulations and Synthesized Results

为了计算BRDF，我们需要知道这个大面积入射光的散射情况。然而，直接模拟这样一个大面积的光会需要很高的计算成本。因此，我们：

使用小尺度模拟的叠加来模拟大尺度入射场。

可以理解为，用很多个小的手电筒（高斯光束）来覆盖一个区域，而不是用一个巨大的探照灯。

首先决定手电筒的尺寸（光束的大小，即腰宽 $w$），并在表面上均匀地分布这些手电筒。写成公式，这个手电筒的排列就是网格点${x_s}, {y_t}$，代表每个高斯光束的中心位置。网格间距一般和腰宽一致，确保光的均匀覆盖，并保持较低的发散角。

让每个高斯光束在其中心区域产生相同的电磁场，只是在不同的位置上重复这一效果。

接下来，想要得到大光束的总散射场，这里需要相位因子来进行“调整”和“叠加”，详细请回看公式(39)。

最终，Combining with Eq. 39, the scattered fields in the far field region corresponding to the pair of directions $(w_i,w_o)$ are given by:

$$ \begin{aligned} \mathbf{E}\left(\omega_i, \omega_o\right) & =\sum_{s=1}^n \sum_{t=1}^n e^{j k \mathbf{p}{s t} \cdot\left(\omega_i+\omega_o\right)} \mathbf{E}{s t}\left(\omega_o\right) \\ \mathbf{H}\left(\omega_i, \omega_o\right) & =\sum_{s=1}^n \sum_{t=1}^n e^{j k \mathbf{p}{s t} \cdot\left(\omega_i+\omega_o\right)} \mathbf{H}{s t}\left(\omega_o\right) \end{aligned} \tag{41} $$

where $\mathbf{E}, \mathbf{H}$ refer to the far field quantities only associated with directions (without the $e^{-j k r} / r$ term).

Lastly, we can compute the surface BRDF value as

$$
f_r\left(\omega_i, \omega_o\right)=\frac{\frac{1}{2}\left|\mathbf{E}\left(\omega_i, \omega_o\right) \times \mathbf{H}\left(\omega_i, \omega_o\right)^*\right|}{\Phi_i \cos \theta_r}
\tag{42}
$$

where the incident power $\Phi_i$ is computed by integrating the incident irradiance over the surface:

$$
\Phi_i=\frac{1}{2} \int_S\left|\left[\mathbf{E}_i\left(\mathbf{r}^{\prime}\right) \times \mathbf{H}_i\left(\mathbf{r}^{\prime}\right)^*\right] \cdot \mathbf{n}\right| d \mathbf{r}^{\prime}
\tag{43}
$$

where $\mathbf{n}$ is the surface normal at the macro scale ( $+\mathbf{z}$ ). Note that Eq. 42 and Eq. 43 can also be applied in single simulations, where $\Phi_i$ is computed from a single Gaussian beam.