简介面向光学与大气科学领域研究者这份Matlab资源围绕Mie散射理论提供计算Mie散射系数an、bn、cn、dn的完整实现可服务于球形颗粒散射特性建模与分析。内容针对需要理解Mie理论、开展消光效率、散射效率模拟的中高级MATLAB用户也适合课堂教学演示与初学者对照学习。压缩包仅2个文件含1个m脚本和1个txt说明整体大小1KBm脚本为核心计算代码细致体现系数求解流程txt文件用于参数定义与调用说明。已有517人学习下载。通过学习该资源读者可以获取可直接运行的Mie系数计算脚本输入波长、粒径及折射率等参数即可获得an、bn等系数同时结合txt说明理解各系数的物理含义与推导关系比如an、bn决定前向/后向散射强度而cn、dn可视为对参数的导数用于相位结构分析。在此基础上还能进一步扩展散射光强分布、极化特性等可视化仿真或与实验数据对比验证。代码结构清晰简洁二次修改门槛低方便嵌入自定义光学仿真流程适合科研与教学场景快速复现理论结果。1. 做光学仿真绕不开的 Mie 散射这套 MATLAB 代码能直接救你的消光光谱在颗粒光学仿真里Mie 散射是绕不开的基础计算环节。不管你是表征金纳米球的局域表面等离激元峰还是估算气溶胶的消光系数第一步都得用 Mie 理论算一个均匀球在平面波下的散射。我这里拆的这份 MATLAB 代码包就是一套采用 ABCD 传输矩阵结构的 Mie 散射计算代码用相对精简的函数体量实现了复数折射率、任意介质背景以及粒径/波长批处理。适合正在做光学仿真、材料表征或者想细扣数值稳定性的读者——新手可以直接调函数出图熟手可以拿它做二次开发不需要额外工具箱这也是它比很多网上流传的旧版 mie.m 更顺手的原因。2. 拆开代码的骨架Mie 系数怎么算ABCD 步骤到底做了什么2.1 理论计算绕不开的参数组合折射率、尺寸参数和结果坐标Mie 理论并不神秘它是把平面波在球坐标系里做级数展开利用球 Bessel 函数加边界连续条件求散射场的展开系数。最后得到的结果不是一条拟合曲线而是一组可以直接用于验证实验的物理量。严格来说计算前要准备三个物理量入射介质折射率 n_med空气取 1.0水取 1.33 左右粒子折射率 n_part复折射率虚部代表吸收以及半径 r 和波长 λ。注意 r 和 λ 必须用同一长度单位代码包第一步做的就是下面这件事% 计算归一化尺度参数 x 和粒子内部尺度参数 y x 2 * pi * n_med * r / lambda; % 外场尺度参数无量纲 m n_part / n_med; % 相对折射率复数 y x * m; % 内场尺度参数已包含折射率影响这里 x 是决定后续级数截断阶数的核心参数。x 很小比如 0.01时粒子处于瑞利散射区偶极项 a1 占绝对主导x 超过 10 之后高阶多极项必须保留否则计算结果会明显偏小。代码包里取截断阶数 Nmax 的方式是经验式的Nmax x 4*x^(1/3) 2这个经验公式在 Bohren 和 Huffman 的经典教材里也是默认做法。输出部分最常见的三个效率是 Qext消光效率、Qsca散射效率、Qabs吸收效率三者满足 Qext Qsca Qabs 的守恒关系。这个关系是后面排查结果的基石也是我每次跑完数据后第一眼要看的东西。效率值都是归一化到粒子的几何投影面积 πr² 上的无量纲量和截面 Cext 的换算公式是 Cext Qext * πr²。2.2 Mie 系数 a_n 和 b_n 的物理含义以及为什么要用递推Mie 散射的核心是解出两组展开系数电系数 a_n 和磁系数 b_n。a_n 对应电多极子对散射场的贡献b_n 对应磁多极子贡献。在小粒子极限下a1 占主导这就是瑞利散射模型能近似成立的原因粒子变大的过程本质上是高阶 a_n、b_n 逐渐参与贡献的过程。其实这两个系数直接按定义写并不难公式里就是几组 Bessel 函数和 Hankel 函数的比值但直接算数值上很容易翻车。原因在于 Bessel 函数 j_n(x) 在 n 远大于 x 时会剧烈振荡并趋向零Hankel 函数 h_n^(1)(x) 则趋向无穷这两个量做除法会直接得到 NaN。所以工程实现里普遍用对数导数法。定义 D_n(z) ψ_n(z) / ψ_n(z)其中 ψ_n(z) z * j_n(z)。这个 D_n 可以按递推关系从高阶向低阶稳定地求出% 从 Nmax 向下递推 D_n(y)这是数值稳定的方向 D zeros(1, Nmax 1); for n Nmax : -1 : 1 D(n) n / y - 1 / (D(n1) n / y); end代码里的方向注释是重点。Bessel 函数的递推分向上和向下两种向上递推在 n 增大时很容易指数爆炸向下递推配合对数导数则能把误差抑制在可控范围。我见过不少人在网上找的 Mie 代码算到 x 20 以上就出现锯齿原因基本都是递推方向选反了。温度、压强这些模型参数不影响这段递推它只依赖折射率 m 和尺度参数 x。2.3 为什么叫 ABCD 传输矩阵版它和单文件老代码的差别网上流传的老版 Mie 代码通常只有一个 mie.m逻辑全部堆在一个函数里算完球体就结束想扩展成多层球或者核壳结构得大改。而这套代码包的设计思路不同它用传输矩阵的视角把球内外电磁场的耦合关系拆成矩阵连乘的形式。ABCD 四个子块分别对应球内场系数与外场系数之间的线性映射边界条件被翻译成矩阵元素最后 a_n、b_n 只是这个矩阵两个分量的比值。这样做的好处很实际如果后续要做核壳结构只需在矩阵连乘链上增加一层传输矩阵而不是重写整个 Bessel 递归如果要做非均匀折射率粒子也只需要把粒子层离散成若干薄层每层算一次矩阵再连乘。这个代码包作为基础版本已经把单层球的矩阵形式写清楚了我看它的函数注释时能明显感受到是为扩展留了接口的。2.4 文件结构先摸一遍动手之前知道每个脚本干嘛文件类型用途mie_abcd.m核心函数计算单球 Qext/Qsca/Qabs 与散射振幅 S1/S2mie_spectrum.m封装脚本波长扫描循环调用核心函数mie_size_dist.m封装脚本粒径分布加权计算run_demo.m演示脚本一键出图适合第一次验证README.txt说明文档参数单位约定与调用示例动手前我建议先打开 run_demo.m 跑一遍看到图再读核心函数。很多人一上来就啃 mie_abcd.m 的递推细节结果被 D_n 的递推公式绕晕其实先看输出格式反而更容易理解输入参数的意义。3. 跑通第一组结果从加载代码到画出效率曲线3.1 环境准备代码包放进 MATLAB 路径避免脚本找不到函数拿到代码包后第一步不是双击运行而是把整个文件夹加入 MATLAB 的搜索路径。否则 MATLAB 会在当前工作目录找不到 mie_abcd 函数直接报 Undefined function。常见做法是右键文件夹选择“添加到路径”或者用命令% 将代码包所在目录加入搜索路径familyfolder 替换为你的实际路径 addpath(genpath(D:\work\mieabcd_matlab));这里用 genpath 而不是 addpath是因为代码包里可能还有子目录genpath 会递归添加所有子孙目录。如果以后又要清理路径用 rmpath(genpath(...)) 就能整条移除。接着检查一下 MATLAB 版本。这套代码只用到了基础函数 besselj、besselh以及复数数组运算所以 R2016b 以后的版本基本都能跑。我没有在太老的版本上验证过但代码逻辑里没有用 matlab 2026b 才有的新特性这一点可以放心。确认版本之后直接运行 run_demo.m% 运行演示脚本观察输出的三条曲线是否正确 run_demo;如果顺利你应该看到一个窗口弹出横轴是波长单位纳米纵轴是 Qext/Qsca/Qabs三条曲线在可见光范围有明显的峰位。看到这个图说明代码包本身已经跑通接下来可以开始改参数。3.2 单球计算详解金纳米球在空气中的 Mie 散射demo 跑通之后我建议你从单球场景开始熟悉参数。下面这段是我平时最常用的基准测试——金纳米球半径 50 nm在空气中波长从 400 nm 扫到 800 nm% 金纳米球基准测试波长扫描 lambda linspace(400, 800, 200); % 波长单位 nm r_gold 50; % 半径单位 nm n_air 1.0; % 空气折射率 n_gold interp1(...); % 金的色散数据需外部导入 Qext zeros(size(lambda)); for i 1 : numel(lambda) n_part n_gold(i); % 当前波长的复折射率 [Qext(i), ~, ~] mie_abcd(n_air, n_part, r_gold, lambda(i)); end这段代码里最容易出错的地方是折射率数据的对应关系。金的介电常数是波长相关的不能用一个常数代替。如果你手头有 Johnson 和 Christy 的实验数据直接线性插值到你的波长网格即可。我一般会额外做一个检查n_gold 的虚部在可见光区必须是正的表示吸收如果你导入的数据虚部是负的后面所有吸收峰都会变成负值。还有一点值得提醒mie_abcd 函数中 r 和 lambda 的单位保持一致。我上面的代码都用了纳米所以 n_med 与 n_part 是无量纲折射率尺度参数 x 不受单位影响。如果你在另一处用了微米这里也要跟着改成微米否则 x 会差 1000 倍结果自然是错的。3.3 读取输出参数到底返回了哪些量怎么变成可发表的图mie_abcd 的完整返回参数是 [Qext, Qsca, Qabs, S1, S2, g]。前三个是效率后两个是散射振幅矩阵元素最后是不对称因子 g。S1 和 S2 通常不是直接用因为在某一个散射角 θ 下这两个复数值描述了散射光两个偏振分量的振幅和相位。绘制散射相位函数时用的是角强度 I(θ) (|S1|² |S2|²) / (2 * x²)这个量才是实验上看到的角分布。我建议第一次跑代码时就把 S1、S2 的结构打印出来确认它们是随角度变化的复数数组而不是标量。很多新人在这里踩坑以为返回的是强度直接取模平方结果数量级对不上。不对称因子 g 是另一个有用的输出它描述散射光朝前还是朝后偏。g 接近 1 表示强前向散射接近 -1 表示强后向散射。这个量在生物组织光学里用得很多如果你算的是血红细胞散射g 值一眼就能判断结果是否合理。4. 实战计算波长扫描和粒径谱的构建4.1 波长扫描的批处理写法避免新手常见的预分配错误单球单波长跑通后实际工作中更常见的是整条光谱的扫描。前面 demo 里的 for 循环写法虽然能跑但在 MATLAB 里性能一般。如果扫描点从 200 个增加到 2000 个最好用向量化思路。% 向量化波长扫描一次调用计算整条谱线 lambda linspace(400, 900, 500).; % 列向量 500 个波长点 n_med 1.0 * ones(size(lambda)); % 介质折射率保持不变 n_part interp1(lam_data, n_data, lambda, pchip); % 粒子复折射率插值 [Qext, Qsca, Qabs] mie_abcd(n_med, n_part, r, lambda);这段代码的关键在于把第一个输入 n_med 写成与 lambda 同尺寸的数组。mie_abcd 内部会自动按元素计算每个波长点的结果从而避免 for 循环。我拆这个包时发现它的向量化接口处理得比较好但前提是你传入的 n_part 是复数数组且尺寸与 lambda 一致。如果你只传一个标量折射率函数内部会自动扩展但这也容易隐藏单位错误所以我自己习惯全部传数组。4.2 粒径分布的计算从单粒径到分布加权的思路实际样品通常不是单一粒径而是服从某种分布比如对数正态分布。此时效率值不能只算一个 r要做分布加权。常见做法是生成一组粒径网格对每个粒径算 Mie 效率再按概率密度加权% 对数正态粒径分布的加权计算 r_grid linspace(10, 200, 200); % 粒径范围单位 nm mu 80; sigma 0.3; % 对数正态参数 pdf 1 ./ (r_grid * sigma * sqrt(2*pi)) .* ... exp(-(log(r_grid) - log(mu)).^2 / (2 * sigma^2)); Qext_dist zeros(size(r_grid)); for i 1 : numel(r_grid) [Qext_dist(i), ~, ~] mie_abcd(n_med, n_part, r_grid(i), lambda_fix); end Qext_weighted trapz(r_grid, Qext_dist .* pdf) / trapz(r_grid, pdf);这里用 trapz 做数值积分而不是直接 sum因为 r_grid 是均匀网格trapz 更接近连续积分的定义。加权之后得到的 Qext_weighted 就是分布平均后的等效消光效率。做材料表征时这个值和 DLS 测出的粒径分布结合能帮你反推样品的真实光学响应。采样粒径网格的密度也有讲究。如果你的分布很窄200 个点足够如果分布很宽且出现拖尾网格间距要小一些否则峰值会被平滑掉。4.3 大尺寸参数下的稳定递推把 Nmax 往上推的两个硬条件当粒子半径大或介质折射率高时x 可能超过 100此时 Nmax 会到 100 以上Bessel 函数的递推稳定性就成了生死线。代码包内部已经做了向下递推但在极端情况下还有两个硬条件需要你注意。第一内部计算全部要用双精度浮点在 MATLAB 里默认就是 double不要手贱转成 single。我见过有人为了省内存把整个数组转单精度结果在 x 50 时 Qabs 直接变成负值。第二递推公式中 D(n) n / y - 1 / (D(n1) n / y) 在 y 的实部接近整数时会遇到接近奇点的位置。这个情况并不常见但如果你算的折射率正好让 y 的实部落在整数附近结果里会出现孤立尖峰。排查手段是扫描粒径时观察结果有没有不连续跳变如果有就把粒径步长加密重跑确认尖峰是否仍然存在。4.4 输出图的保存和后续处理别让数据只活在弹出窗口里跑完扫描之后我习惯把原始数据直接存成 .mat 文件避免下次重新计算。保存命令很简单% 保存计算结果方便后续绘图或对比 save(mie_result_au50.mat, lambda, Qext, Qsca, Qabs, g);如果你要把图导出到论文里用 exportgraphics 比 print 更清晰而且能避免中文乱码% 导出高清图分辨率 300dpi适合直接放入论文 exportgraphics(gcf, mie_au50_spectrum.png, Resolution, 300);这里提一句matlab 图像处理方向的同学还经常会把 Mie 计算得到的相位函数作为点扩散函数模板用于散射成像的仿真。代码包里 S1、S2 的输出格式正好支持这个用法不用额外转换。5. Mie 计算避坑总结数值不稳定、单位符号等常见问题Mie 计算的坑比想象中多我自己调试这个代码包时也翻过几次车下面几条是最具代表性的按现象、原因、解决的方式列出来希望能帮你省下几个晚上的调试时间。5.1 现象算大粒径时 Qext 出现 NaN 或锯齿状跳变有一次我算直径 2 μm 的二氧化硅球在可见光下的消光光谱结果在 600 nm 附近 Qext 直接变成了 NaN。检查代码逻辑没问题又回退到 x 40 以下一切正常。原因x 超过 60 之后Bessel 函数的动态范围过大双精度浮点的 15 位有效数字不够用递归过程中误差被放大最后出现上溢或下溢。解决设置实验性的输入截断和缩放因子。常见做法是在 Bessel 函数计算中引入对数标度或者把 x 的超出部分单独处理。这个代码包里用向下递推已经缓解了一部分但如果你需要在 x 大于 100 的情况下计算建议改用 Wiscombe 的稳定化版本或增加更高精度的底层函数。5.2 现象吸收峰形状对但符号反了Qabs 是负的我做金纳米球的时候第一次算出的吸收光谱在 520 nm 附近有个峰但 Qabs 在峰位处是负值。原因折射率虚部的符号约定搞反了。很多文献里折射率写为 m n - ik但 MATLAB 中的指数约定是 exp(-iωt)对应的吸收要求虚部为正即 m n ik。如果你把虚部输成负的计算出的吸收就变负。解决在调用 mie_abcd 之前强制检查虚部符号% 检查并修正折射率虚部符号 if imag(n_part) 0 warning(折射率虚部为负代码已自动取反); n_part conj(n_part); end从那以后我写任何 Mie 脚本开头都会加一段类似的符号检查。5.3 现象Qext 大于 Qsca Qabs能量守恒不成立有次我算完一组数据发现 Qext 和 QscaQabs 差了 5%不是数值噪声能解释的。原因我看结果时直接把三个量分别调用了三次函数每次调用时参数类型不一致其中一次 r 传成了微米另一次传成了纳米导致三个量对应的 x 不同。解决一次调用只算一组参数并且把结果打包输出。不要为了省事在命令行分别调用三次再去手动对准这不是函数的问题是我自己手滑。把调用封装成一个函数传入同一个 r 和 lambda保证三个效率从同一组内部变量计算得到。5.4 现象散射角分布图非常平滑但数值小得离谱画散射相位函数图时角度分布曲线形状合理但纵轴数量级从 10 掉到了 0.001。原因: 我把角度单位搞混了。mie_abcd 默认要求角度用弧度但我在外部用角度制绘图直接把 0 到 180 的度数数组传进去了三角函数全部失真。解决统一用弧度或者在图内转换别在数据层面混用。我的习惯是创建一个角度向量然后用 deg2rad 转成弧度再传给函数。如果你需要极坐标图polarplot 也能直接用弧度数据省去转来转去。5.5 现象粒子在液体里算出来谐振峰位置整体偏移测试水中的金纳米球结果等离激元峰位置比实验文献偏了 30 nm 左右怎么调粒径都补不回来。原因相对折射率没算对。我只把粒子的折射率传进去介质折射率仍然用了空气的 1.0等于假设粒子在空气里。解决介质折射率必须填实际环境值水取 1.33油取 1.45 左右然后函数内部会自动算相对折射率。这一步代码包没有自动判别纯粹是输入参数的问题。换了介质之后峰位偏移立即消失和文献对上了。6. 最后再做一遍验证用光学定理和几何极限校验你刚跑出的数据代码跑通只是第一步数据对不对需要系统校验。我最后分享一个自己固定执行的验证套路三分钟就能做完。第一重验证是能量守恒。Qext 应该严格等于 Qsca Qabs。如果在某个波长差出超过 1%基本可以判定是数值问题而不是物理问题。代码包里三个量是在同一次计算中得到的理论上不会差但当我修改了递推阶数后偶尔会出现微小偏差。只要偏差在 1e-6 以内就可以忽略超过就回头检查 Nmax 是否取够。第二重验证是瑞利极限。当 x 小于 0.01 时Qsca 应该与 x 的四次方成正比Qabs 与 x 成正比。你可以取一组极小粒径跑一遍把 Qsca 和 x 做双对数图如果斜率接近 4说明系数实现正确。这一步能暴露递推公式中符号错误因为瑞利极限对系数非常敏感。第三重验证是几何光学极限。当 x 远大于 1比如大于 1000且粒子无吸收时Qext 应该趋近于 2这就是所谓的消光悖论——大粒子消光截面是几何截面的两倍。如果你算出来一个很大的值偏离 2那说明递推到高阶项时出了问题。% 三重验证脚本的骨架按顺序跑 % 1. 能量守恒 assert(max(abs(Qext - (Qsca Qabs))) 1e-6, 能量守恒不满足); % 2. 瑞利极限x 很小时 Qsca ~ x^4 x_small 2 * pi * 1.0 * 0.001 / 633; % 极小球 [~, Qsca_small, ~] mie_abcd(1.0, 1.5, 0.001, 633); ratio Qsca_small / x_small^4; % 应接近常数 % 3. 大粒子极限无吸收大球 Qext 趋近 2 [Qext_large, ~, ~] mie_abcd(1.0, 1.50i, 50000, 633); disp(Qext_large); % 期望值在 2 附近这套验证跑完我对计算结果就比较有信心了。从那以后我每次换一个新的 Mie 代码包或者改折射率数据源都会强制走一遍这三重验证十几分钟换几个小时的安心。数值计算这东西不怕慢就怕错得莫名其妙。希望这次的拆解能帮到你少走一点我走过的弯路。本文还有配套的精品资源点击获取