简介面向随机分布二维柱散射问题研究者的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, dtypenp.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 取多了不会立刻出错但矩阵规模按 2nmax1 线性增长N 根柱的总维度是 N×(2nmax1)。程序里建议把截断阶数做成函数允许外部传入覆盖值方便后面做截断收敛性扫描。注意以上只是初始截断随机分布高填充率时柱间多次散射会让高阶模再次被激发最终应以 R/T 随 nmax 不再变化作为判定标准。3. 多柱耦合用加法定理组装矩阵一次求解全部散射系数多级散射理论的核心是任取一根柱 j它感受到的入射场等于外部平面波加上其他柱 l 散射到 j 位置的场。把每个“其他柱”的散射波用柱 j 的局部坐标重新展开就出现加法定理里的汉克尔函数和相位因子。将所有柱、所有角动量模放在一起得到一个维度为 N×(2*nmax1) 的线性方程组。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×MM2*nmax1。注意 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, dtypenp.complex128) rhs np.zeros(n_scat * m_size, dtypenp.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 | M2*nmax1 | 总自由度 | 复数矩阵内存 | 直接求解预期 | |----|------|-------------|----------|--------------|--------------| | 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), dtypenp.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% 时优先检查加法定理组装里的相位符号。# 单柱N1A 应为 1rhs 等于平面波展开系数 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再用单柱和空介质极限排除组装错误最后才把计算资源投入到随机样本平均中。本文还有配套的精品资源点击获取 多级散射理论计算随机二维柱阵反射与透射的Python实现
简介面向随机分布二维柱散射问题研究者的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, dtypenp.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 取多了不会立刻出错但矩阵规模按 2nmax1 线性增长N 根柱的总维度是 N×(2nmax1)。程序里建议把截断阶数做成函数允许外部传入覆盖值方便后面做截断收敛性扫描。注意以上只是初始截断随机分布高填充率时柱间多次散射会让高阶模再次被激发最终应以 R/T 随 nmax 不再变化作为判定标准。3. 多柱耦合用加法定理组装矩阵一次求解全部散射系数多级散射理论的核心是任取一根柱 j它感受到的入射场等于外部平面波加上其他柱 l 散射到 j 位置的场。把每个“其他柱”的散射波用柱 j 的局部坐标重新展开就出现加法定理里的汉克尔函数和相位因子。将所有柱、所有角动量模放在一起得到一个维度为 N×(2*nmax1) 的线性方程组。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×MM2*nmax1。注意 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, dtypenp.complex128) rhs np.zeros(n_scat * m_size, dtypenp.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 | M2*nmax1 | 总自由度 | 复数矩阵内存 | 直接求解预期 | |----|------|-------------|----------|--------------|--------------| | 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), dtypenp.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% 时优先检查加法定理组装里的相位符号。# 单柱N1A 应为 1rhs 等于平面波展开系数 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再用单柱和空介质极限排除组装错误最后才把计算资源投入到随机样本平均中。本文还有配套的精品资源点击获取