均匀圆球Mie散射的MATLAB仿真实现

Mie散射理论描述了电磁波遇到球形粒子时的散射、吸收与消光行为，该理论是研究气溶胶、水滴、粉尘等微粒光学特性的重要基础。下面这套MATLAB程序，将Mie散射理论完整转化为可执行的代码，并集成了可视化模块，方便用户直观理解不同参数对散射行为的影响。

%% 均匀圆球Mie散射计算程序% 描述：计算均匀圆球的Mie散射参数（散射效率、吸收效率、消光效率等）clear; clc; close all;%% 1. 参数设置lambda = 0.55e-6; % 入射光波长 (m) - 550nm绿光r = 1e-6; % 粒子半径 (m) - 1μmm = 1.33   0.01i;% 粒子复折射率 (水: 1.33@550nm)n_medium = 1.0;% 周围介质折射率 (空气)theta = linspace(0, 180, 361); % 散射角范围 (度)%% 2. Mie散射计算核心函数[mie] = calculateMie(lambda, r, m, n_medium, theta);%% 3. 结果可视化% 效率因子随尺寸参数变化plotEfficiencyFactors();% 散射相函数plotPhaseFunction(mie);% 散射强度分布plotScatteringIntensity(mie);% 电场分布plotElectricField(mie);%% 4. 输出关键参数fprintf('===== Mie散射计算结果 =====');fprintf('波长: %.1f nm', lambda*1e9);fprintf('粒子半径: %.1f μm', r*1e6);fprintf('粒子折射率: %.2f   %.2fi', real(m), imag(m));fprintf('尺寸参数 x = %.4f', mie.x);fprintf('相对折射率 m = %.4f   %.4fi', real(mie.m), imag(mie.m));fprintf('===== 效率因子 =====');fprintf('消光效率 Q_ext = %.4f', mie.Qext);fprintf('散射效率 Q_sca = %.4f', mie.Qsca);fprintf('吸收效率 Q_abs = %.4f', mie.Qabs);fprintf('散射不对称因子 g = %.4f', mie.g);fprintf('===== 其他参数 =====');fprintf('散射截面 σ_sca = %.4e m²', mie.sigma_sca);fprintf('吸收截面 σ_abs = %.4e m²', mie.sigma_abs);fprintf('消光截面 σ_ext = %.4e m²', mie.sigma_ext);%% Mie散射计算核心函数function [mie] = calculateMie(lambda, r, m, n_medium, theta)% 计算基本参数k = 2 * pi * n_medium / lambda;% 波数x = k * r; % 尺寸参数m_complex = m / n_medium;% 相对复折射率% 确定求和项数n_max = round(max(x   4*x^(1/3)   2, 10));% 经验公式% 初始化变量an = zeros(n_max, 1);% Mie系数a_nbn = zeros(n_max, 1);% Mie系数b_ncn = zeros(n_max, 1);% 其他系数dn = zeros(n_max, 1);% 其他系数% 计算Riccati-Bessel函数及其导数[psi, psi_prime, xi, xi_prime] = riccatiBessel(x, n_max);[psi_m, psi_m_prime] = riccatiBessel(m_complex*x, n_max);% 计算Mie系数for n = 1:n_max% 系数a_nnumerator_a = m_complex * psi_m(n) * psi_prime(n) - psi(n) * psi_m_prime(n);denominator_a = m_complex * psi_m(n) * xi_prime(n) - xi(n) * psi_m_prime(n);an(n) = numerator_a / denominator_a;% 系数b_nnumerator_b = psi_m(n) * psi_prime(n) - m_complex * psi(n) * psi_m_prime(n);denominator_b = psi_m(n) * xi_prime(n) - m_complex * xi(n) * psi_m_prime(n);bn(n) = numerator_b / denominator_b;end% 计算效率因子Qext = 0;Qsca = 0;Qabs = 0;asym = 0;for n = 1:n_maxterm_ext = (2*n   1) * real(an(n)   bn(n));term_sca = (2*n   1) * (abs(an(n))^2   abs(bn(n))^2);term_asym = (n*(n 2)/(n 1)) * real(conj(an(n))*an(n 1)   conj(bn(n))*bn(n 1)) ...  ((2*n 1)/(n*(n 1))) * real(an(n)*conj(bn(n)));Qext = Qext   term_ext;Qsca = Qsca   term_sca;asym = asym   term_asym;endQext = (2/(x^2)) * Qext;Qsca = (2/(x^2)) * Qsca;Qabs = Qext - Qsca;g = (4/(x^2*Qsca)) * asym;% 不对称因子% 计算散射截面sigma_sca = Qsca * pi * r^2;sigma_abs = Qabs * pi * r^2;sigma_ext = Qext * pi * r^2;% 计算散射振幅函数S1 = zeros(size(theta));S2 = zeros(size(theta));for i = 1:length(theta)theta_rad = deg2rad(theta(i));for n = 1:n_maxpi_n = legendreP(n, cos(theta_rad));tau_n = n * cos(theta_rad) * pi_n - (n 1) * legendreP(n-1, cos(theta_rad));S1(i) = S1(i)   (2*n 1)/(n*(n 1)) * (an(n)*pi_n   bn(n)*tau_n);S2(i) = S2(i)   (2*n 1)/(n*(n 1)) * (an(n)*tau_n   bn(n)*pi_n);endS1(i) = S1(i) * sin(theta_rad);S2(i) = S2(i) * sin(theta_rad);end% 计算散射强度I_sca = (abs(S1).^2   abs(S2).^2) / (k^2 * r^2);% 存储结果mie = struct();mie.x = x;mie.m = m_complex;mie.Qext = Qext;mie.Qsca = Qsca;mie.Qabs = Qabs;mie.g = g;mie.sigma_sca = sigma_sca;mie.sigma_abs = sigma_abs;mie.sigma_ext = sigma_ext;mie.S1 = S1;mie.S2 = S2;mie.I_sca = I_sca;mie.theta = theta;end%% Riccati-Bessel函数计算function [psi, psi_prime, xi, xi_prime] = riccatiBessel(x, n_max)% 初始化数组psi = zeros(n_max, 1);psi_prime = zeros(n_max, 1);xi = zeros(n_max, 1);xi_prime = zeros(n_max, 1);% 初始值 (n=0)psi(1) = sin(x);psi_prime(1) = cos(x);xi(1) = sin(x) - 1i*cos(x);% ξ_0 = ψ_0 - iχ_0, χ_0 = -cos(x)xi_prime(1) = cos(x)   1i*sin(x);% 递推计算 (n=1,2,...)for n = 1:n_max-1% Riccati-Bessel函数 ψ_npsi(n 1) = ((2*n 1)/x) * psi(n) - psi(n-1);% 导数 ψ_n'psi_prime(n 1) = psi(n) - (n 1)/x * psi(n 1);% Riccati-Hankel函数 ξ_n = ψ_n   iχ_nxi(n 1) = ((2*n 1)/x) * xi(n) - xi(n-1);% 导数 ξ_n'xi_prime(n 1) = xi(n) - (n 1)/x * xi(n 1);endend%% Legendre多项式计算function P = legendreP(n, x)% 递归计算Legendre多项式if n == 0P = ones(size(x));elseif n == 1P = x;elseP = ( (2*n-1)*x.*legendreP(n-1, x) - (n-1)*legendreP(n-2, x) ) / n;endend%% 效率因子随尺寸参数变化绘图function plotEfficiencyFactors()lambda = 0.55e-6; % 波长 (m)r = linspace(0.01e-6, 5e-6, 100); % 粒子半径 (m)m = 1.33   0.01i;% 粒子复折射率n_medium = 1.0;% 介质折射率x = 2 * pi * n_medium * r / lambda; % 尺寸参数Qext = zeros(size(r));Qsca = zeros(size(r));Qabs = zeros(size(r));for i = 1:length(r)mie = calculateMie(lambda, r(i), m, n_medium, [0]);Qext(i) = mie.Qext;Qsca(i) = mie.Qsca;Qabs(i) = mie.Qabs;endfigure;semilogy(x, Qext, 'b-', 'LineWidth', 2); hold on;semilogy(x, Qsca, 'r-', 'LineWidth', 2);semilogy(x, Qabs, 'g-', 'LineWidth', 2);xlabel('尺寸参数 x');ylabel('效率因子');title('Mie散射效率因子随尺寸参数变化');legend('消光效率 Q_{ext}', '散射效率 Q_{sca}', '吸收效率 Q_{abs}');grid on;% 标记瑞利散射区和米氏散射区line([0.1, 0.1], ylim, 'Color', 'k', 'LineStyle', '--');text(0.12, max(ylim)*0.9, '瑞利散射区', 'FontSize', 10);line([10, 10], ylim, 'Color', 'k', 'LineStyle', '--');text(10.2, max(ylim)*0.9, '米氏散射区', 'FontSize', 10);line([100, 100], ylim, 'Color', 'k', 'LineStyle', '--');text(102, max(ylim)*0.9, '几何光学区', 'FontSize', 10);end%% 散射相函数绘图function plotPhaseFunction(mie)theta_rad = deg2rad(mie.theta);% 计算相函数 P(θ) = (|S1|^2   |S2|^2)/(σ_sca * 2πk^2)P = (abs(mie.S1).^2   abs(mie.S2).^2) / (2 * pi * mie.sigma_sca * (2*pi/mie.x)^2);figure;polarplot(theta_rad, P, 'b-', 'LineWidth', 2);title('散射相函数 P(θ)');rlim([0 max(P)*1.1]);figure;plot(mie.theta, P, 'b-', 'LineWidth', 2);xlabel('散射角 θ (度)');ylabel('相函数 P(θ)');title('散射相函数');grid on;end%% 散射强度分布绘图function plotScatteringIntensity(mie)figure;polarplot(deg2rad(mie.theta), mie.I_sca, 'r-', 'LineWidth', 2);title('散射强度分布 I(θ)');rlim([0 max(mie.I_sca)*1.1]);figure;plot(mie.theta, mie.I_sca, 'r-', 'LineWidth', 2);xlabel('散射角 θ (度)');ylabel('散射强度 I(θ)');title('散射强度分布');grid on;% 前向散射和后向散射特写figure;subplot(2,1,1);idx_forward = find(mie.theta <= 30);plot(mie.theta(idx_forward), mie.I_sca(idx_forward), 'r-', 'LineWidth', 2);xlabel('散射角 θ (度)');ylabel('前向散射强度');title('前向散射 (0°-30°)');grid on;subplot(2,1,2);idx_backward = find(mie.theta >= 150);plot(mie.theta(idx_backward), mie.I_sca(idx_backward), 'b-', 'LineWidth', 2);xlabel('散射角 θ (度)');ylabel('后向散射强度');title('后向散射 (150°-180°)');grid on;end%% 电场分布绘图function plotElectricField(mie)% 计算电场分布 (简化模型)theta_rad = deg2rad(mie.theta);E_parallel = real(mie.S2 .* exp(1i*kr));% 平行分量E_perpendicular = real(mie.S1 .* exp(1i*kr)); % 垂直分量figure;polarplot(theta_rad, abs(E_parallel), 'b-', 'LineWidth', 2); hold on;polarplot(theta_rad, abs(E_perpendicular), 'r-', 'LineWidth', 2);title('电场强度分布 |E|');legend('平行分量', '垂直分量');rlim([0 max([abs(E_parallel), abs(E_perpendicular)])*1.1]);figure;plot(mie.theta, abs(E_parallel), 'b-', 'LineWidth', 2); hold on;plot(mie.theta, abs(E_perpendicular), 'r-', 'LineWidth', 2);xlabel('散射角 θ (度)');ylabel('电场强度 |E|');title('电场强度分布');legend('平行分量', '垂直分量');grid on;end

程序功能概览

1. 核心计算模块

首先聚焦核心计算模块，该模块由三个函数协同完成：

• calculateMie函数——这是Mie散射计算的主力函数。它从最基础的参数起步：尺寸参数 \( x = \frac{2\pi r n_m}{\lambda} \) 和相对折射率 \( m = \frac{n_p}{n_m} \)，然后求解Riccati-Bessel函数及其导数，从而得到Mie系数 \( a_n \) 和 \( b_n \)，进一步导出效率因子 \( Q_{\text{ext}}, Q_{\text{sca}}, Q_{\text{abs}} \) 以及不对称因子 \( g \)，最后给出散射振幅函数 \( S_1(\theta), S_2(\theta) \) 和散射强度分布。
• riccatiBessel函数——专门用于计算Riccati-Bessel函数系列。它通过递推关系依次计算 \( \psi_n(x), \psi'_n(x), \xi_n(x), \xi'_n(x) \)，为Mie系数的分母与分子提供关键的数值支撑。
• legendreP函数——负责计算Legendre多项式 \( P_n(\cos\theta) \)。同样采用递归实现，思路清晰且效率较高。

2. 可视化模块

该程序提供了四类可视化手段，让抽象的数据变得直观易懂：

• 效率因子随尺寸参数变化——绘制 \( Q_{\text{ext}}, Q_{\text{sca}}, Q_{\text{abs}} \) 随 \( x \) 的变化曲线，并自动标注瑞利散射区、米氏散射区以及几何光学区的分界线，方便用户快速区分不同散射区域。
• 散射相函数——同时采用极坐标和直角坐标展示 \( P(\theta) \)，后者更有利于观察角度细节。
• 散射强度分布——除全角度分布图外，还专门截取前向（0°-30°）和后向（150°-180°）两个关键区段进行特写，这两个方向在实际应用中具有重要意义。
• 电场强度分布——分别绘制平行分量和垂直分量的幅值分布，有助于深入理解偏振特性。

3. 关键物理量

• 效率因子：\( Q_{\text{ext}} \)（消光效率，代表散射+吸收的总效果）、\( Q_{\text{sca}} \)（散射效率）、\( Q_{\text{abs}} \)（吸收效率）
• 不对称因子 \( g \)：描述散射的方向倾向，正数表示前向散射占优
• 散射截面：\( \sigma_{\text{sca}}, \sigma_{\text{abs}}, \sigma_{\text{ext}} \)
• 散射振幅函数：\( S_1(\theta), S_2(\theta) \)

物理原理与算法

1. Mie散射基本理论

Mie散射理论本质上是麦克斯韦方程组在球坐标系下的精确解。其核心思路是将散射场用矢量球谐函数展开，展开系数 \( a_n, b_n \) 包含了粒子尺寸、折射率和入射波长的全部信息。散射振幅函数最终可表示为：

2. 效率因子计算公式

3. 特殊区域的散射特性

• 瑞利散射区 ( \( x \ll 1 \) )：
  
  此时散射效率正比于 \( x^4 \)，因此蓝光散射强度远高于红光——这正是天空呈现蓝色的物理根源。
• 米氏散射区 ( \( x \sim 1 \) )：需要完整的Mie级数计算，不能做任何近似处理。
• 几何光学区 ( \( x \gg 1 \) )：散射主要由反射、折射和衍射主导，消光效率趋近于2，即几何光学极限。

计算结果分析

1. 效率因子随尺寸参数变化

• 瑞利散射区（\( x < 0.1 \)）：\( Q_{\text{sca}} \propto x^4 \)，蓝光散射明显强于红光
• 米氏散射区（\( 0.1 < x < 50 \)）：出现共振峰，散射效率可达2-4倍
• 几何光学区（\( x > 50 \)）：\( Q_{\text{ext}} \approx 2 \)

2. 散射相函数特征

• 前向散射（\( \theta \approx 0^\circ \)）：强度最大，随角度增加迅速衰减
• 后向散射（\( \theta \approx 180^\circ \)）：强度次之
• 90°散射：强度最小
• 不对称因子 \( g \) 为正时，表明前向散射占据主导地位

3. 典型应用场景

• 大气科学：用于计算气溶胶光学厚度、云滴散射特性
• 海洋光学：分析海水散射行为
• 生物医学：实现细胞散射测量
• 纳米光子学：研究金属纳米颗粒的局域场增强效应

扩展功能

1. 多分散体系计算

实际样品往往不是单一粒径，而是服从某种分布。下面的示例展示了如何利用粒径分布权重来计算平均散射效率：

% 计算多分散体系的散射特性radii = logspace(-9, -5, 100); % 粒径分布 (0.001-10 μm)weights = rayleigh_pdf(radii); % 瑞利分布权重Qsca_a vg = 0;for i = 1:length(radii)mie = calculateMie(lambda, radii(i), m, n_medium, [0]);Qsca_a vg = Qsca_a vg   weights(i)*mie.Qsca;end

2. 偏振特性计算

利用散射振幅函数可以方便地计算偏振度：

% 计算偏振度Polarization = (abs(S1).^2 - abs(S2).^2) ./ (abs(S1).^2   abs(S2).^2);

3. 吸收特性分析

通过改变波长并扫描吸收效率，即可获得吸收光谱：

% 计算吸收截面谱wa velengths = linspace(0.3e-6, 1.0e-6, 100); % 0.3-1.0 μmQabs_spectrum = zeros(size(wa velengths));for i = 1:length(wa velengths)mie = calculateMie(wa velengths(i), r, m, n_medium, [0]);Qabs_spectrum(i) = mie.Qabs;endfigure;plot(wa velengths*1e9, Qabs_spectrum, 'b-', 'LineWidth', 2);xlabel('波长 (nm)');ylabel('吸收效率 Q_{abs}');title('吸收效率光谱');grid on;

总结

本MATLAB程序完整实现了均匀圆球的Mie散射计算，归纳起来具有以下特点：

• 物理完备：涵盖完整的Mie级数展开，以及Riccati-Bessel、Legendre等特殊函数的计算
• 可视化丰富：效率因子、相函数、强度分布、电场分布一应俱全
• 参数灵活：波长、粒径、折射率均可自由调整，适应不同研究需求
• 计算高效：采用递推算法，避免高阶函数直接求值，提升运行速度
• 扩展性强：稍加修改即可应用于多分散体系、偏振分析、吸收光谱等多种场景