简介：面向随机分布二维柱散射问题研究者的MATLAB程序包，基于多级散射理论计算反射与透射特性。该程序将复杂散射系统分解为多级散射网络，逐步求解每个节点的散射矩阵，并通过统计平均获得随机柱排列下的反射率和透射率，适用于纳米光学、光子学、声学等领域中类似随机介质问题的快速模拟。压缩包仅含1个.m脚本，大小约1KB，核心代码简洁，便于阅读、修改与嵌入其他数值实验。当前已有173人浏览学习，对正在接触多级散射或蒙特卡洛统计方法的研究生、科研人员具有一定参考价值。通过运行该脚本，用户可复现随机二维柱散射的反射/透射计算流程，理解散射网络的构建思路与统计平均的实现方式，在此基础上调整柱体参数、入射条件或扩展为更大规模随机样本，辅助后续课题研究。
1. 多级散射理论算随机二维柱阵反射透射，先从“多次散射”而不是“几何叠加”开始
随机分布二维柱散射的反射和透射算不准，十有八九是因为把多根柱体当成了独立散射体做叠加。实际场景里，柱间距小到几个波长时，一根柱的散射场会成为另一根柱的入射场，来回迭代的结果既不能忽略也不能只看前两阶。多级散射理论（MST）的正路是把每根柱的响应压缩成一个 T 矩阵，再用加法定理把“柱间传输”拼成一组线性方程组，解一次就得到所有柱的散射系数，最后用远场积分给出反射率和透射率。下面按可复现的 Python/NumPy 路径讲：单柱 T 矩阵、多柱相互作用矩阵、随机样本统计，以及程序写完必须过的三个检验。适用人群是做计算电磁、随机介质、光子结构或薄膜涂层设计的人，知道 Mie 散射会更快。
2. 二维柱散射的数学化：T 矩阵、展开阶数与截断误差
在进入多级散射理论之前，先把单根柱的问题钉死。柱体沿 z 轴无限长，入射平面波在 xy 平面内传播，问题退化为二维标量波。电场或磁场只剩一个分量，因此可以直接按 TE/TM 两种极化分开算。这样的好处是每根柱的散射性质可以用一个对角 T 矩阵描述，散射系数和入射系数只差一个复数因子，多体方程写起来就紧凑很多。
2.1 场展开：入射波、散射波和局部柱坐标
在柱 j 的局部坐标 (r_j, φ_j) 中，入射平面波 exp(i k·r) 可以展开为柱函数 J_m(k r) e^{i m φ} 的无穷级数，散射场则用第一类汉克尔函数 H_m^{(1)}(k r) 展开。J_m 在原点有限，代表入射和驻波部分；H_m^{(1)} 在无穷远处呈向外传播的行波，代表从柱体出去的散射波。把场写成这两种基底，边界条件就能逐项匹配。
2.1.1 为什么用柱函数做基底
二维圆截面天然匹配柱坐标，贝塞尔函数是圆域上的本征解。更重要的是，任意一根柱的散射场在其他柱的局部坐标里，可以通过加法定理重新展开成 J_m，这样“柱 l 的散射波落到柱 j 上”就变成了一个耦合矩阵元素。如果改用平面波展开，圆边界条件会变成卷积积分，程序复杂度和截断误差都更难控制。
2.1.2 平面波展开系数的相位约定
对于一个从方向 φ0 入射、幅度为 1 的平面波，在柱 j 中心的局部展开系数一般写成：
c_jm = exp(i k·r_j) × i^m × exp(-i m φ0)
其中 exp(i k·r_j) 是坐标原点移到柱 j 带来的相移，i^m 来自平面波的柱函数展开，exp(-i m φ0) 是入射方向约定。不同程序可能把符号反过来，导致散射幅度绕 x 轴镜像对称。写程序前先把这组约定固定下来，否则和后续远场公式对不上。
2.2 单柱 T 矩阵：Mie 系数的程序实现
圆形介质柱的 T 矩阵是对角的，第 m 个角动量的散射系数只和它自己有关。TM 极化（电场沿柱轴）下，单柱 T 矩阵元由两对贝塞尔函数组成：
import numpy as np from scipy.special import jv, yv, jvp, yvp def mie_tmatrix_tm(ka, mr, nmax): """ 单根介质柱在 TM 极化下的 T 矩阵对角元。 ka : 背景介质中的波数 k0 乘以柱半径 a mr : 柱体折射率与背景折射率之比 nmax: 角动量截断阶数，返回数组长度为 nmax """ t = np.zeros(nmax, dtype=np.complex128) for n in range(nmax): jx = jv(n, ka) jp = jvp(n, ka, 1) jmx = jv(n, mr * ka) jmp = jvp(n, mr * ka, 1) hx = jx + 1j * yv(n, ka) hp = jp + 1j * yvp(n, ka, 1) num = jmx * jp - mr * jx * jmp den = jmx * hp - mr * hx * jmp t[n] = num / den return t逻辑上，分子是柱内场与背景场在边界上的导数匹配，分母是柱内场与向外辐射场的匹配。t[n] 的模接近 1 表示柱体对大尺寸共振散射体作用强，模接近 0 表示该阶角动量几乎不参与散射。代码里用 jvp(n, x, 1) 求一阶导数，第二个参数取 1 表示求导阶数。TE 极化把 mr 换成 1/mr 再套同样结构即可，因为磁边界条件里折射率比值是反的。
2.3 截断阶数 nmax 的选取：经验公式与误差表
二维柱散射的截断阶数主要由尺度参数 ka 决定。经验上取 nmax = ceil(ka + 4 × (ka)^(1/3)) + 5，再额外加几阶作为安全余量，基本能让截断误差落到 1e-6 以下。 | ka | 推荐 nmax | 实数方程组宽度（单个柱） | 备注 | |----|-----------|--------------------------|------| | 0.5 | 9 | 19 | 小半径柱，低阶占优，但相位精度敏感 | | 2 | 13 | 27 | 常见微米颗粒/太赫兹波段 | | 10 | 24 | 49 | 需要开始关注矩阵条件数 | | 20 | 36 | 73 | 高阶模参与，矩阵变大，内存翻倍 |
nmax 取多了不会立刻出错，但矩阵规模按 2nmax+1 线性增长，N 根柱的总维度是 N×(2nmax+1)。程序里建议把截断阶数做成函数，允许外部传入覆盖值，方便后面做截断收敛性扫描。
注意：以上只是初始截断，随机分布高填充率时，柱间多次散射会让高阶模再次被激发，最终应以 R/T 随 nmax 不再变化作为判定标准。
3. 多柱耦合：用加法定理组装矩阵，一次求解全部散射系数
多级散射理论的核心是：任取一根柱 j，它感受到的入射场等于外部平面波加上其他柱 l 散射到 j 位置的场。把每个“其他柱”的散射波用柱 j 的局部坐标重新展开，就出现加法定理里的汉克尔函数和相位因子。将所有柱、所有角动量模放在一起，得到一个维度为 N×(2*nmax+1) 的线性方程组。
3.1 耦合方程与矩阵结构
定义 f_jm 为作用在柱 j 上的第 m 个模的有效入射展开系数，a_jm = t_m f_jm 为对应的散射系数。柱 l 贡献给柱 j 的耦合项在局部坐标下是：
G_{jl}^{mn} = H_{m-n}^{(1)}(k0 × R_lj) × exp(-i (m-n) φ_lj)
其中 R_lj 是柱 l 到柱 j 的距离，φ_lj 是该距离矢量的方位角。把贡献累加起来，方程组写成：
f_jm = c_jm + Σ_{l≠j} Σ_n G_{jl}^{mn} × t_{|n|} × f_ln
对角线是 1，矩阵 A 的维度为 N×M，M=2*nmax+1。注意 c_jm 是 2.1.2 节里的平面波展开系数，只有外场那一项，不包含其他柱的贡献。G 矩阵是稠密的，因为多级散射在原理上允许任意模之间耦合。
3.2 组装程序：位置、波数、阶数到系统矩阵
下面的函数把柱位列表、背景波数、入射角和 T 矩阵组装成 A 矩阵和右端项：
from scipy.special import hankel1 def build_msa_matrix(pos, k0, phi0, nmax, t_matrix): """ pos : (N, 2) 柱心坐标数组 k0 : 背景波数，2*pi/波长 phi0 : 入射角，单位弧度 nmax : 截断阶数 t_matrix : 单柱 T 矩阵，长度为 nmax，索引取绝对值对应阶数 返回 A, rhs，求解后得到每个柱的有效入射系数 f """ n_scat = len(pos) m_size = 2 * nmax + 1 a_mat = np.eye(n_scat * m_size, dtype=np.complex128) rhs = np.zeros(n_scat * m_size, dtype=np.complex128) kx = k0 * np.cos(phi0) ky = k0 * np.sin(phi0) for j in range(n_scat): phase = np.exp(1j * (kx * pos[j, 0] + ky * pos[j, 1])) for mi, m in enumerate(range(-nmax, nmax + 1)): row = j * m_size + mi rhs[row] = phase * (1j ** m) * np.exp(-1j * m * phi0) for l in range(n_scat): if l == j: continue dx = pos[l, 0] - pos[j, 0] dy = pos[l, 1] - pos[j, 1] r_lj = np.hypot(dx, dy) phi_lj = np.arctan2(dy, dx) for ni, n in enumerate(range(-nmax, nmax + 1)): nu = m - n h_val = hankel1(abs(nu), k0 * r_lj) if nu < 0: h_val *= (-1) ** nu g_val = h_val * np.exp(-1j * nu * phi_lj) col = l * m_size + ni a_mat[row, col] -= g_val * t_matrix[abs(n)] return a_mat, rhs组装顺序是先遍历所有源柱 l，再遍历所有角动量 n，把加法定理系数累加到目标柱 j 上。代码里对负阶汉克尔函数做了符号修正，否则相位会整体差一个 (-1)^n。R_lj 为 0 时汉克尔函数发散，所以随机位置生成时必须有最小距离约束，这会在下一章处理。
3.3 求解器选型与矩阵体检
对中等规模构型，直接用 numpy.linalg.solve。矩阵是稠密复数非对称的，维度增大时内存增长很快。 | N | nmax | M=2*nmax+1 | 总自由度 | 复数矩阵内存 | 直接求解预期 | |----|------|-------------|----------|--------------|--------------| | 50 | 10 | 21 | 1050 | 约 18 MB | 秒级 | | 200 | 10 | 21 | 4200 | 约 282 MB | 秒到分钟级 | | 200 | 24 | 49 | 9800 | 约 1.5 GB | 分钟级，需关注内存 |
求解之前先检查 A 矩阵的条件数。条件数超过 1e12，后面 R/T 的小数后几位往往不可信。常见原因是柱间距离过近、nmax 不足，或者柱体材料参数导致 T 矩阵接近奇点。应对方法是剔除重叠随机位置，并把 nmax 提高两到三阶再重算。
提示：如果 N 超过几百，直接组装稠密矩阵会非常吃力。可以先把加法定理矩阵按距离截断，或用快速多极方法加速矩阵向量积，再配合 GMRES 迭代求解。
4. 随机分布与反射透射的统计计算：样本、平均和误差条
真实随机介质里没有“唯一”的反射率。同样的体分比、同样的膜厚，换一组随机位置结果会抖。程序要做的不是算一个构型交差，而是随机位置生成、逐构型求解、系综平均三步走。这也是随机分布二维柱散射程序里最容易漏掉的环节。
4.1 生成不重叠的随机柱位
随机柱位必须满足最小间距约束，这个约束不是几何洁癖，而是数学上的刚需：距离过小的柱体之间，加法定理中的汉克尔函数会让矩阵病态，6.3 节的条件数检查会直接报警。一个可用的拒绝采样函数如下：
def random_positions_no_overlap(n_scat, lx, ly, min_gap, seed): rng = np.random.default_rng(seed) pos = np.zeros((n_scat, 2)) placed = 0 tries = 0 max_tries = 5000 * n_scat while placed < n_scat and tries < max_tries: tries += 1 cand = rng.uniform([0.0, 0.0], [lx, ly]) if placed == 0 or np.min( np.hypot(pos[:placed, 0] - cand[0], pos[:placed, 1] - cand[1])) > min_gap: pos[placed] = cand placed += 1 tries = 0 if placed < n_scat: raise RuntimeError("排布失败：减小填充率或扩大区域") return posmin_gap 一般取柱直径的 1.1 倍以上，具体由柱半径和波长共同决定。拒绝采样在低填充率下效率很高，但超过 40% 填充率时容易长时间跑不出新位置，这时改用抖动网格或逐步挤压算法更合适。随机种子用可复现的 default_rng，方便调试和对照。
4.2 从散射系数算反射和透射功率
多级散射方程解出的散射系数 a 是柱 j 第 m 个角动量的复数幅度。远场某一方向 φ 的散射幅度可以按柱函数渐近式合成：
from scipy.integrate import trapezoid def far_field_phi(phi_grid, pos, a_vec, nmax, k0): n_scat = len(pos) m_size = 2 * nmax + 1 f_phi = np.zeros(len(phi_grid), dtype=np.complex128) for ip, phi in enumerate(phi_grid): s = 0j for j in range(n_scat): phase_shift = np.exp(-1j * k0 * ( pos[j, 0] * np.cos(phi) + pos[j, 1] * np.sin(phi))) for mi, m in enumerate(range(-nmax, nmax + 1)): a_val = a_vec[j * m_size + mi] s += (phase_shift * a_val * (1j ** m) * np.exp(1j * m * phi)) f_phi[ip] = np.sqrt(2.0 / (np.pi * k0)) * np.exp(1j * np.pi / 4) * s return f_phi代码里的 phase_shift 把远场参考点从柱中心平移到整体坐标原点，保证各个柱的散射波到达观察点时相位一致。得到全角度的 F(φ) 后，把反射半平面和透射半平面的 |F(φ)|² 分别做梯形积分，再归一化到两者之和，就得到该构型的反射率和透射率。这个定义把柱阵看成平面上的散射屏，避开了有限波束宽度归一化的歧义，适用于随机薄层的工程评估。
4.3 系综平均、样本数与停止条件
随机分布的程序必须做多构型平均，例如生成 50 到 200 个不同柱位样本，每个样本独立求解并记录 R 和 T，最后输出均值和标准差。判定停止的常用条件是滑动平均的相对变化小于某个阈值，比如最近 20 个样本的均值与总体均值差小于 1%。 | 参数 | 示例值 | 说明 | |--------------|--------|-------------------| | 柱半径 a | 0.5 | 单位统一为波长 | | 背景折射率 | 1.0 | 空气或真空 | | 柱折射率 | 1.5 | 无耗介质 | | 填充率 | 10% | 高填充率需提高 nmax | | 入射角 φ0 | π/2 | 垂直照射散射屏 | | 样本数 | 80 | 结合误差条决定 |
每个样本都要重新生成随机位置、组装矩阵、求解和积分。最小 demo 程序可以先固定 20 个样本跑通流程，正式研究时再把样本数提高到标准误差小于目标精度。注意保存每次构型的 R、T 和随机种子，后续排查奇异矩阵或异常反射率时可以定点重放。
5. 写完程序先做这三个检验：单柱对照、空介质极限和阶数扫描
程序跑出第一张 R/T 曲线后，先不要直接调整物理参数，按下面三个顺序做验证，缺一个都可能把错误当成物理结果。
5.1 单柱极限对照解析 Mie 公式
在随机分布程序中只放一根柱，关闭其他柱的耦合，把散射系数 a_m 与单柱 T 矩阵直接做对比。散射系数应当严格等于 t_m × 入射展开系数，这是方程组退化成单柱的必经之路。远场积分得到的散射截面也应当和 Mie 理论解析值一致，误差超过 1% 时优先检查加法定理组装里的相位符号。
# 单柱：N=1，A 应为 1，rhs 等于平面波展开系数 pos_one = np.array([[0.0, 0.0]]) a_mat, rhs = build_msa_matrix(pos_one, k0, phi0, nmax, t_matrix) # a_mat 应接近单位矩阵，f 应接近 rhs如果单柱都合不上，多柱统计没有任何意义。此时最可能出错的是负阶汉克尔函数符号、T 矩阵中 mr 的极化定义、或平面波展开系数里的 e^{-i m φ0} 符号。
5.2 空介质与弱散射极限
把柱半径缩小到波长的千分之一，或把柱折射率设为 1.0001，随机阵的反射率应趋近 0，透射率应趋近 1。实际操作中，弱散射极限的精度可以检验随机位置生成和远场积分是否有系统性偏差。柱半径极小时 T 矩阵接近零，耦合矩阵的贡献也接近零，此时多级散射程序被压回了几乎不散射的平凡解，这是最简单的整体链路冒烟测试。
5.3 阶数扫描与条件数扫描
固定一个随机构型，从 nmax 等于推荐值的一半开始，逐步增大到推荐值两倍，记录每个阶数下的 R、T 和矩阵条件数：
for nmax_test in [10, 14, 18, 22, 26]: a_mat, rhs = build_msa_matrix(pos, k0, phi0, nmax_test, t_matrix) f_vec = np.linalg.solve(a_mat, rhs) # 由 f_vec 计算 a_vec，再求 R、T cond_val = np.linalg.cond(a_mat) print(nmax_test, R, T, cond_val)R 和 T 的判读标准是最后两档阶数的变化小于 0.001。如果一直在摆动，说明柱间多次散射明显激发了更高阶模，需要在初始 nmax 基础上再加几阶，同时观察条件数是否急剧上升。条件数超过 1e13 时，即使 R/T 看起来收敛，单个样本的相位误差也可能被放大，程序应主动输出警告。
这三个检查都通过后，再回到系综平均流程，把样本数从 20 增加到 100 以上。R 和 T 的标准误差随样本数按平方根下降，但单构型计算时间复杂度更高，所以先用阶数扫描确定 nmax，再用单柱和空介质极限排除组装错误，最后才把计算资源投入到随机样本平均中。
本文还有配套的精品资源，点击获取
