做颗粒光散射计算的人，大概率逃不开一个问题：我要算的颗粒是任意形状，不是标准球，Mie理论给不了答案，应该怎么办？如果你也卡在这，那我强烈建议你亲手写一遍离散偶极子近似（DDA）的MATLAB实现。我记得第一次跑通自己的DDA代码，把球形颗粒的消光效率和Mie理论叠在一张图上对比时，那种“这玩意儿真的能算出来”的兴奋感，比直接用现成工具强烈得多。这套方法的核心思路非常朴素：把任意形状的颗粒切成一个个小体积单元，每个单元等效成一个偶极子，然后自洽求解所有偶极子在入射光照射下的响应，最后再汇总成全空间的散射场。它对形状没有要求，球、椭球、立方体、团簇、甚至随机粗糙表面颗粒都能算，这也是它在纳米光学、大气辐射、生物医学光散射领域里被广泛使用的原因。
这篇整理会从物理原理一路讲到MATLAB工程落地，包括网格生成、极化率赋值、耦合方程构建、迭代求解、截面和远场角分布的后处理，最后再补上我调试过程中踩过的一堆坑。适合刚接触DDA、想自己从零搭一套可运行代码的人参考，也适合已经会用现成开源包、但想知道“里面到底干了什么”的朋友。
1. DDA在解一个什么样的“散射问题”
1.1 为什么不用Mie理论：形状自由是核心卖点
Mie理论是光散射分析的基石，但它只对均匀球体有解析解。一旦颗粒变成椭球、立方体、圆柱、甚至两个球粘在一起的二聚体，Mie理论就无能为力了。虽然也有针对椭球的扩展解析方法，但适用范围依然极其有限，更别说随机取向、表面粗糙、多层结构这些更复杂的场景。
DDA的思路是“以离散换形状自由”。它把颗粒占用的空间离散成规则格点，每个格点是一个小立方体，体积是d³。只要d远小于波长，颗粒内部和外部的电磁场在这个尺度上变化足够缓慢，那么每个小立方体对光场的响应就可以用一个点偶极子来等效描述。这样，再复杂的形状在程序代码里就只体现为一个“索引数组”——坐标为格点、在颗粒内部则标记为1，在外部则标记为0。只要这个掩模数组能生成出来，任意形状都能参与计算。
相比有限时域差分（FDTD）和有限元（FEM）这类通用电磁场求解器，DDA的离散方式非常“节约”：它只在颗粒体积内布点，颗粒外部不需要额外网格。颗粒在空气中时，这个优势特别明显，不浪费任何一个格点在入射波传播路径上。
1.2 偶极子离散与自洽方程
设入射场是平面波，在颗粒内部第i个格点位置r_i处的入射场为：
E_inc(r_i) = E0 exp(i k · r_i)
颗粒的折射率是m_p，周围介质折射率是m_med，通常用相对折射率m = m_p / m_med。每个偶极子的极化率为α_i，偶极矩为P_i。问题是：P_i不仅受入射场驱动，还受颗粒内部其它所有偶极子辐射场的驱动。所以每一个偶极矩都要满足一个耦合方程：
P_i = α_i [ E_inc,i + E_scat,i ]
其中E_scat,i是其他所有偶极子在第i个偶极子位置产生的散射场之和。这是典型的“多体自洽”问题，本质上是一个3N维的线性方程组，N是偶极子总数。整理成矩阵形式就是：
(1/α_1, 0, 0, ...；0, 1/α_2, 0, ...；... ) P + A P = E_inc
这里的A矩阵包含了偶极子两两之间的相互作用，对角线块是极化率的倒数，非对角线块是点偶极子辐射场的格林张量。
有个特别容易困惑的地方：为什么方程里i=j那一项用的是1/α，而不是格林函数本体？因为点偶极子在自己所在位置的辐射场是发散的，物理上偶极子自身的“退极化场”已经包含在极化率α的定义里了。换句话说，α是考虑了自场修正之后的有效极化率。所以在构建方程时，对角线不能直接套用非对角线的格林函数公式，而必须用1/α来表征每个偶极子的自响应。很多初写代码的人在这一步翻车，把对角线也填成某个格林函数表达式，算出来结果完全失真。
2. 从零搭建MATLAB求解器：五个关键模块
2.1 网格生成：从三维坐标到颗粒掩模
第一步是在三维空间中生成规则格点。我习惯用ndgrid，而不是meshgrid。因为meshgrid在三维下会把第一维和第二维的顺序搞错，导致坐标和索引对不上。ndgrid生成的X、Y、Z数组维度是[ nx, ny, nz ]，和空间坐标的直觉一致。
生成一个球形颗粒的掩模代码如下：
% 物理参数（长度单位统一用微米） lambda = 0.5; % 入射光波长 d = 0.010; % 偶极子间距，单位微米，即10 nm R = 0.050; % 颗粒半径，单位微米，即50 nm % 生成覆盖颗粒的最小矩形网格 half = ceil(R / d); xx = (-half : half) * d; [X, Y, Z] = ndgrid(xx, xx, xx); % 球体掩模 mask = (sqrt(X.^2 + Y.^2 + Z.^2) <= R); % 提取颗粒内部格点坐标 pos = [X(mask), Y(mask), Z(mask)]; N = size(pos, 1);这里有几个细节值得注意。一是边界条件用“形心是否落在球内”来判断，这简单但会让阶梯边界产生约半个格点的尺寸误差。如果颗粒很小、格点数很少，这个误差会直接影响截面结果的精度。更精细的做法是对球边界做体积加权修正，但对于刚上手的第一版代码，直接用形心判断即可。二是格点数量N和颗粒体积的关系：颗粒真实体积是4/3πR³，离散后的体积是N*d³，两者之间会有偏差。在调试收敛性时，我会输出这两个体积做对比，如果偏差超过百分之几，说明网格太粗。
2.2 极化率：把材料参数映射到每个偶极子
每个偶极子的极化率α由材料折射率和格点体积一起决定。最简单的方案是Clausius-Mossotti关系，这也是DDA里最基础的极化率模型：
m = 1.5 + 0.02i; % 相对折射率 V = d^3; % 格点体积 epsilon = m^2; alpha_CM = 3 * V * (epsilon - 1) / (epsilon + 2); alpha = alpha_CM * ones(N, 1);这个公式的物理含义是：把一个介质小球放在外场中，其内部场会因为退极化效应而修正，极化率自然和介电常数、体积强相关。CM关系形式简单，对折射率不太高、格点相对波长不太大的场景已经足够。
但如果你希望和成熟开源包的精度对齐，CM关系在高折射率或较大格点时误差会明显偏大，因为真实格点上的偶极子受到的晶格修正不只是简单的退极化场。业界常用的做法是Lattice Dispersion Relation（LDR）修正，它会根据入射波的波矢、格点间距和材料折射率对极化率做进一步校准。LDR实现的细节涉及多项展开系数，首次动手不建议自己推导，我的建议是先用CM版本跑通全流程，确认DDA整体架构没问题之后，再去查阅文献或参考成熟开源包源码，把极化率替换成LDR版本。这样分两步走，能最大程度减少把“物理问题”和“工程bug”混在一起排查的痛苦。
2.3 耦合矩阵的显式构建与内存边界
对于第一个可用版本，我建议直接显式构建整个耦合矩阵M，然后用MATLAB的“\”直接求解，或者交给bicgstab迭代。这样逻辑最清晰，也容易调试。
相互作用格林张量的构建是整段代码的核心。在MATLAB里，先把偶极子坐标pos拆成三者分量，然后双重循环遍历所有偶极子对。对于第i个偶极子和第j个偶极子，点偶极子的散射场作用矩阵有一个标准闭式公式。我当时手写的时候，喜欢把公式拆成“径向点积部分”和“垂直部分”来避免出错。
M = zeros(3*N, 3*N); E0vec = [1, 0, 0]; kdir = [0, 0, 1]; % 入射方向沿 z k = 2 * pi / lambda; for i = 1:N row0 = (i-1)*3; % 对角线块：自作用 M(row0+1:row0+3, row0+1:row0+3) = eye(3) / alpha(i); for j = 1:N if i == j continue; end rvec = pos(i,:) - pos(j,:); r = norm(rvec); rhat = rvec / r; g = 1/r^2 - 1i*k/r; % 点偶极子格林张量 Aij = exp(1i*k*r)/r * ( ... k^2 * (rhat.'*rhat - eye(3)) + ... (3*rhat.'*rhat - eye(3)) * g ); % 方程左端 M * P = E_inc，因此非对角块是 -Aij col0 = (j-1)*3; M(row0+1:row0+3, col0+1:col0+3) = -Aij; end end这段代码可以跑，但循环是O(N²)，N在几百时没有问题，N到两三千时已经很吃力，到五千以上基本没法等。再看内存：一个3N乘3N的复数稠密矩阵，每个元素占16字节，N=1000时是144MB，N=2000时接近580MB，N=5000时直接超过3.6GB。所以显式矩阵方案的上限大约在N=1000到1500左右，超过这个量级，必须换成“矩阵向量算子+迭代求解”的进阶模式。这部分我在第五章展开。
需要特别说明的是坐标和索引的分量排布约定。我习惯把未知向量P打包成[px(1:N); py(1:N); pz(1:N)]的结构，这样第i个偶极子的三个分量分别在第i、N+i、2N+i的位置。矩阵构建时二维块的行列索引要严格对应这个排布。另一种常见排布是[px1, py1, pz1, px2, py2, pz2, ...]的交错式，两种都可以，但后处理提取分量时容易把自己绕晕。建议从一而终。
2.4 求解器：直接解还是迭代解
小N场景直接用M\b，简单可靠。当N变大、矩阵无法显式存储时，必须走迭代路线。两个推荐：gmres和bicgstab。
MATLAB内置的gmres只需要提供矩阵向量乘的函数句柄，不需要转置矩阵向量乘，这对DDA特别友好，因为格林张量矩阵不是对称矩阵，手写转置算子容易出错。bicgstab收敛通常更快，但它同时需要A函数和A'函数，第一版代码不建议上来就搞这个。所以我的顺序是：N小于1500用M\b，N更大用gmres配一个Afun函数句柄。
迭代求解的一大痛点是收敛判据。DDA矩阵的条件数会随着颗粒尺寸、折射率虚部的变化剧烈波动。金属颗粒在等离子体共振波长附近，极化率很大，矩阵条件数明显恶化，迭代很容易“看起来收敛了，结果却是错的”。我自己遇到过两次这种伪收敛，后来都会额外检查残差对应的物理量，比如用求出的P向量重算一次C_ext，看它和迭代中间值是否稳定，而不是单纯相信求解器返回的flag。
2.5 后处理：截面与远场方向图
解出P向量后，最常见的物理量是消光截面、吸收截面和远场角分布。
消光截面C_ext可以用光学定理计算：
% 入射场各偶极子处的值 Einc = zeros(N, 3); Einc(:,1) = exp(1i * k * (pos * kdir.')); % 这里的kdir是行向量...细节略 for i = 1:N Einc(i,:) = E0vec * exp(1i * k * dot(pos(i,:), kdir)); end % 偶极矩P向量拆分 px = Pvec(1:N); py = Pvec(N+1:2*N); pz = Pvec(2*N+1:3*N); % C_ext = 4*pi*k/|E0|^2 * sum( imag( P_i . conj(Einc_i) ) ) Cext = 4 * pi * k * sum(imag(px .* conj(Einc(:,1)) + ... py .* conj(Einc(:,2)) + ... pz .* conj(Einc(:,3))));注意这里已经假设入射场振幅|E0|=1，并采用了空间相位因子exp(i k·r)的约定。如果你换用了其他约定，公式中的符号会变化，这是跨程序对比时最容易出问题的地方。
远场散射方向图则是把每个偶极子作为辐射源，对方向r̂_s做相参叠加：
E_s(r̂_s) ∝ k² * exp(ikr)/r * Σ_i exp(-i k r̂_s · r_i) * [ r̂_s × (r̂_s × P_i) ]
要对一组给定的方位角(θ, φ)计算角分辨散射强度，只需要再写一个循环遍历方向采样点。散射截面C_sca等于远场强度对全立体角积分再除以入射强度。吸收截面则用C_abs = C_ext - C_sca得到，这样比直接用近似公式更稳健，也避开了“吸收公式里要不要含偶极子相干项”的坑。我第一次写吸收截面时就踩过这个坑：直接对每个偶极子单独算辐射功率再求和，忽略了偶极子之间的相干项，算出来的C_sca和远场精确积分差了一大截。从那以后我后处理一律走远场积分的路子。
3. 上手跑通示例：50nm介质球
3.1 参数选择与网格密度估算
一个稳妥的起步算例是：半径R=50nm、相对折射率m=1.5+0.02i的介质球，入射波长λ=500nm。这样颗粒尺寸远小于波长，散射处于瑞利-米氏过渡区，Mie理论可以给出精确参考值，便于验证。
网格间距d的选择经验是：|m| k d < 1。这里k = 2π/λ，d = 10nm时代入得到m约1.5、λ=500nm，|m|kd约0.19，满足条件。N约为419。这个规模用显式矩阵求解在普通笔记本上几秒钟就能完成，非常适合当第一个测试用例。
如果你希望更高精度，可以把d减小到5nm。但那样N会涨到接近4200，显式矩阵内存达到约2.5GB，普通机器已经很难受，必须切到算子模式和迭代求解。一种折中的中间档是d=6nm左右，N约1200到1500，显式矩阵约250MB，还能一战。
这个趋势说明一个现实：DDA的精度是靠网格密度换来的，而网格密度直接决定内存和计算量。所以好的DDA实现必须在“显式便利”和“算子高效”之间选择合适路径。
3.2 与Mie理论对比验证
跑完代码，下一步一定是拿着结果和Mie理论做对比。这一步建议不要只对比一个波长，而是扫一个波长范围，把C_ext随λ变化的曲线叠在一起看。
我的对比流程是这样的：固定颗粒半径和折射率模型，在λ=400nm到900nm之间取20到30个点，每个点都调用DDA计算C_ext，同时调用一个Mie理论函数算出C_ext_ref。然后把两条曲线画在一起，计算相对误差。你会发现大部分区域的误差都在可接受范围内，但在某些特定波长可能会出现系统性偏高或偏低。
出现偏差先不要急着怀疑代码，优先检查三件事：网格d是否满足|m|kd<1的条件；颗粒边界阶梯效应对体积的误差是否超过2%；折射率虚部符号和相位约定是否和Mie参考程序一致。我印象中至少有一次查了一整天，最后发现是Mie参考函数使用的折射率虚部是负号约定，而DDA代码是正号约定，导致吸收项符号反了，误差大得离谱。
3.3 从球体扩展到任意形状
DDA真正的价值在非球形。把网格生成部分的mask表达式换一下，就能立刻从球变成其他形状。比如椭球：
% 椭球半轴 a, b, c mask = (X.^2/a^2 + Y.^2/b^2 + Z.^2/c^2 <= 1);立方体：
mask = (abs(X) <= Lx/2) & (abs(Y) <= Ly/2) & (abs(Z) <= Lz/2);两个球组成的二聚体可以先用两个mask判断再做逻辑或。我做过一个随机取向椭球的代码，只是把pos重新随机旋转一下再生成掩模，整个DDA核心逻辑一行没改。这就是DDA对初学者的最大友好之处：形状只决定掩模，掩模在物理计算里只表现为格点列表和数量。你写好了球体，就写好了所有形状。
4. 常见问题与避坑实录
4.1 误差偏大到底该调哪个参数
这是后台被问得最多的问题。误差来源可以粗暴分成两类：空间离散误差和边界逼近误差。
空间离散误差由d过大引起。颗粒内部场的空间变化由波长和折射率决定，折射率越高，电磁场在颗粒内部变化越剧烈，需要更小的d。经验判据|m|kd<1是底线，实际工程我一般取d <= λ / (15 |m|)，这样C_ext误差通常能控制在百分之几以内。
边界逼近误差则由球形边界上格点“阶梯化”引起。对同一个d，颗粒越小、边界附近的格点占比越高，相对误差越大。比如50nm半径的球用10nm格点，直径方向只有10层格点，边界格点比例很高，误差明显；如果换成200nm半径的球，同样的d，边界占比降低，误差就会小很多。所以小颗粒想要同样精度，必须付出更细网格的代价。
我的习惯是做一个网格收敛测试：把d从20nm逐步减到8nm，观察C_ext是否单调逼近某个稳定值。如果结果还在明显漂移，那就继续加密。做这个测试时不要忘了记录每个d对应的N和内存占用，因为你会很快撞到显式矩阵内存上限。
4.2 迭代求解器卡死或“伪收敛”
用gmres或bicgstab时，最常见的问题是迭代几百步不收敛，或者flag为0但物理结果明显错误。
先说卡死。高折射率金属颗粒、或者颗粒尺寸较大时，矩阵对角占优性下降，迭代收敛会很慢。一个有效的预处理是对角预处理：把整个矩阵按对角块做归一化。对角块恰好是各偶极子的1/α，物理上相当于把偶极子的自响应剥离，剩余的相互作用是短程+辐射项，条件数会改善很多。MATLAB中可以用gmres(Afun, b, restart, tol, maxit, L, U)传入预处理矩阵L和U。简单的Jacobi预处理就够用。
再说伪收敛。gmres返回的残差是基于矩阵方程M P = E_inc的，但物理上我们更关心截面量。有时候残差降到1e-8，C_ext还在随迭代步数缓慢漂移，这是因为C_ext是偶极矩和入射场的点积，某些高频误差分量可能在残差范数里权重很小，却在物理量里权重很大。我处理的方法是：每隔一定迭代步数，把当前的P向量取出来算一次C_ext，画一条C_ext随迭代步数的曲线，看它是否平台。平台震荡剧烈，说明迭代参数没调好；平台单向漂移，说明收敛判据太低。这个习惯后来救了我好几次。
4.3 内存不足的工程化解法
显式矩阵方案在N到1500以上就开始吃力，2万甚至10万偶极子的大规模运算必须彻底放弃A矩阵。这里给出的工程路径是：
一是在不改变物理模型的前提下，把数据类精度从double切到single。但DDA的复平面运算里精度损失可能让消光截面误差放大，所以只对极度内存受限的场景用，且必须和double版本做对比验证。
二是把“显式矩阵”换成“矩阵向量算子”。写一个函数，输入P向量，输出P在方程左端的像。算子内部不需要构建完整矩阵，而是按偶极子对关系循环累加。N在几千时，MATLAB的向量化循环还能扛住。
三是更进一步，利用DDA相互作用矩阵的平移不变性，用FFT加速矩阵向量乘。这一步是大型DDA性能的关键，也是成熟开源包的核心做法。
4.4 单位、符号约定与跨程序对比
单位问题几乎是每个写DDA的人都会踩的坑。我的建议是代码内统一用微米作为长度单位，所有涉及频率、波矢、截面的量都据此推导。不要混用纳米和微米，否则一个10^-3的系数错误会让你怀疑整个物理模型。
符号约定上，DDA有两种主流选择：空间因子exp(-i k·r)和exp(+i k·r)。它们会显著改变插值、格林函数和吸收项的符号。我使用exp(+i k·r)约定，对应的折射率虚部为正代表材料吸收。在和外部Mie程序对比前，先确认双方约定一致。一个小技巧是：算一个已知解析解的小球在真空中的消光效率，虚部设为零且折射率近1，看结果是否趋近于零。这样就可以快速锁定符号方向。
还有折射率应该是相对折射率还是绝对折射率。颗粒在空气里两者相等，但颗粒在水里就完全不同。DDA里的极化率和格林函数都假设背景介质是均匀且无源的，所以入射场、波矢、极化率都必须在背景介质折射率下重新定义。我接触过几个初学者把颗粒绝对折射率直接当成相对折射率用，导致共振波长整体偏移，排查时非常迷惑。
5. 性能优化进阶：不做FFT就慢到没朋友
5.1 矩阵向量算子的实现思路
当N超过两三千，显式矩阵开始不可行，第一步还是实现Afun算子：即不存储A矩阵，而是在函数内动态计算A*x。这样内存占用从O(N²)降到O(N)，时间则是每次矩阵向量乘都是O(N²)循环。
Afun的基本结构是输入x向量，先拆成三个分量px、py、pz，然后对每个偶极子i，循环所有j≠i，累加格林张量作用在P_j上的结果。虽然还是O(N²)，但因为没有保存3N乘3N的复数矩阵，N到5000时内存不再是瓶颈，只是时间涨得厉害。
MATLAB里写的第一个Afun很容易写成三重循环：外层i、内层j、每层里面再算3乘3矩阵乘。这样在N=3000时一次矩阵向量乘可能就要好几分钟，迭代几百步根本无法忍受。必须改成向量化：对每个i，一次性用数组运算算出所有j的位移和格林项。核心是避免在j循环里使用norm函数，改为对整个j向量作减法、求模、广播点积。
我建议第一次写Afun时，先用显式矩阵的结果做验证：随机生成P，分别用A*P和Afun(P)的结果对比，最大误差应小于1e-10。通过这一步再进入迭代求解，能省掉后面好几个小时的排查。
5.2 平移不变结构与FFT加速
DDA的相互作用矩阵M只依赖于两个格点之间的相对位移Δr = r_i - r_j。这意味着，如果所有偶极子都排在一个完整、规则的矩形网格上，那么格林张量矩阵就具有Toeplitz结构，矩阵向量乘可以用三维FFT加速到O(N log N)。
注意这里有个容易忽略的细节：颗粒形状通常是不规则的，偶极子只占据矩形网格的一部分。FFT加速时不能在稀疏的子集上直接做卷积，而必须先把偶极子响应扩展到完整矩形网格上，非颗粒区域的偶极矩设为0，再做三维卷积，最后从结果中提取颗粒区域内格点的部分。这个过程相当于把空外部区域的贡献显式算成了零，物理上完全一致，但FFT必须做一次三维尺寸的卷积，内存占用取决于包围盒尺寸而不是颗粒体积本身。
实际实现FFT版Afun时，要处理的三个格林张量分量分别对三个偶极矩分量做卷积，输出三个散射场分量，总共是9次三维卷积。我建议分阶段完成：先实现一个朴素FFT版本，在颗粒近似填满矩形网格的场景下验证，再用细长颗粒验证稀疏场景下的加速度和内存代价。
FFT加速后，N到数万甚至数十万都成为可能。颗粒出现凹槽或尖锐结构时，包围盒体积可能远大于颗粒体积，FFT的内存优势会被削弱，这也是新型DDA变体试图解决的方向之一，但对大多数应用场景，标准FFT加速已经足够用。
5.3 给新手的优化优先级
如果你想把代码从“能跑”升级到“能算大颗粒”，我建议按照以下优先级推进，不要一上来就攻坚FFT。
第一步，把显式矩阵改成Afun算子，配合gmres跑通N=3000到5000的算力边界。这一步能帮助你理解迭代求解和内存瓶颈。
第二步，优化Afun的向量化，把O(N²)循环做到尽可能快。这里通常能做到一个量级的加速，并且Afun本身可以作为后续FFT版本的验证基准。
第三步，实现FFT矩阵向量乘。这也是工作量最大的一步，建议把格林张量公式拆开写清楚，分块测试卷积结果，不要一步到位。
第四步，加LDR极化率修正。当你已经确信DDA算法本身没问题，再用LDR去抠那百分之几的精度，边际收益最大。
优化要随时回到物理验证：每改一个模块，都要用上一章的50nm球在某个波长下重算一遍，和Mie理论对比。不然你可能在加速的同时悄悄引入了一个符号错误，直到算一个大颗粒时才暴露，那才是最难受的调试场景。
最后再分享一个我个人的小习惯：每次改完网格生成逻辑，我都顺手打印一行体积信息，比较离散体积N*d³和理论体积。这个检查几乎不花时间，却能提前预警很多边界上的隐蔽问题。DDA这个方向看起来公式多、矩阵大，但真正把它跑通之后，你会明显感觉到，所谓“任意形状散射计算”并没有想象中那么高不可攀。如果你也被某个形状奇怪的颗粒困扰着，强烈建议照着上面的代码框架，从最小算例开始，一步步把它算出来。
