企业官网建设流程全解析

别光看公式了！手把手教你用Matlab复现米氏散射曲线（附完整代码）

当你在论文中看到那些优美的米氏散射曲线时，是否曾想过自己动手复现？作为光学测量领域的经典理论，米氏散射的计算确实复杂，但有了Matlab这个强大的工具，一切变得可能。本文将带你从零开始，完整走通米氏散射曲线的计算流程。

1. 环境准备与代码获取

工欲善其事，必先利其器。在开始计算前，我们需要确保Matlab环境配置正确。推荐使用R2018b或更新版本，这些版本对特殊函数的支持更完善。

必备工具箱：

• 核心Matlab（必须）
• Symbolic Math Toolbox（可选，用于公式推导）
• Curve Fitting Toolbox（可选，用于结果分析）

获取代码的几种途径：

1. 直接联系论文作者（最推荐）
2. 学术代码共享平台（如GitHub、ResearchGate）
3. 根据论文描述自行实现

注意：使用他人代码时务必遵守相关许可协议，尊重知识产权。

2. 参数解读与输入准备

米氏散射计算的核心参数包括：

在Matlab中，这些参数通常组织为结构体：

params.d = 100e-9; % 100nm粒子 params.lambda = 532e-9; % 532nm激光 params.m = 1.59 + 0.01i; % PS微球折射率 params.n_medium = 1.33; % 水介质

3. 核心算法实现

米氏散射计算的核心是级数求和，关键在于确定截断项数。经验公式为：

n_max = round(2 + 4*x^(1/3) + x)

其中x为尺寸参数，x = πd/λ。Matlab实现如下：

function [Qext, Qsca, Qabs] = mie_calc(params) % 计算尺寸参数 x = pi * params.d * params.n_medium / params.lambda; m_relative = params.m / params.n_medium; % 确定截断项数 n_max = round(2 + 4*x^(1/3) + x); % 初始化系数数组 a_n = zeros(1, n_max); b_n = zeros(1, n_max); % 计算散射系数 for n = 1:n_max [a_n(n), b_n(n)] = mie_coeff(n, x, m_relative); end % 计算效率因子 Qsca = (2/x^2) * sum((2*(1:n_max)+1).*(abs(a_n).^2 + abs(b_n).^2)); Qext = (2/x^2) * sum((2*(1:n_max)+1).*real(a_n + b_n)); Qabs = Qext - Qsca; end

辅助函数mie_coeff负责计算单个散射系数：

function [a_n, b_n] = mie_coeff(n, x, m) % 使用Ricatti-Bessel函数计算系数 psi_n = sqrt(pi*x/2) * besselj(n+0.5, x); psi_n1 = sqrt(pi*x/2) * besselj(n-0.5, x); chi_n = sqrt(pi*x/2) * bessely(n+0.5, x); D_n = (m*psi_n1 - (n/x)*psi_n) / (m*psi_n - (n/x)*psi_n1); G_n = (psi_n1 - (n/x)*psi_n) / (psi_n - (n/x)*psi_n1); a_n = (D_n/m + n/x)*psi_n - psi_n1; a_n = a_n / ((D_n/m + n/x)*psi_n - psi_n1 + 1i*((D_n/m + n/x)*chi_n - chi_n1)); b_n = (m*G_n + n/x)*psi_n - psi_n1; b_n = b_n / ((m*G_n + n/x)*psi_n - psi_n1 + 1i*((m*G_n + n/x)*chi_n - chi_n1)); end

4. 结果可视化与分析

计算完成后，我们需要将结果直观展示。典型的米氏散射曲线包含三个区域：

1. 瑞利散射区（x << 1）：Qext ∝ x⁴
2. 米氏散射区（x ≈ 1）：剧烈震荡
3. 夫琅禾费区（x >> 1）：趋近于2

绘制消光效率因子随尺寸参数变化的曲线：

% 生成尺寸参数范围 x_values = logspace(-1, 2, 500); % 0.1到100对数均匀分布 % 预分配结果数组 Qext_values = zeros(size(x_values)); % 循环计算 for i = 1:length(x_values) params.d = x_values(i) * params.lambda / (pi * params.n_medium); [Qext, ~, ~] = mie_calc(params); Qext_values(i) = Qext; end % 绘制曲线 figure; loglog(x_values, Qext_values, 'LineWidth', 2); xlabel('尺寸参数 x = \pid/\lambda'); ylabel('消光效率因子 Q_{ext}'); title('米氏散射消光效率曲线'); grid on;

典型输出曲线应显示：

• x<0.3时：陡峭上升（瑞利区）
• 0.3<x<10时：震荡特征（米氏区）
• x>10时：趋近于2（夫琅禾费区）

5. 常见问题排查

在实际复现过程中，你可能会遇到以下问题：

问题1：计算结果出现NaN

• 检查复折射率虚部不为零
• 验证尺寸参数x在合理范围
• 确认介质折射率不为零

问题2：曲线震荡不明显

• 增加尺寸参数采样密度
• 检查级数截断项数是否足够
• 确认复折射率实部与虚部比例适当

问题3：计算速度过慢

• 预分配所有数组
• 考虑使用parfor并行计算
• 对固定参数预先计算并存储

性能优化技巧：

% 使用数组运算替代循环 n = 1:n_max; coeff = (2*n + 1)./(n.*(n + 1)); % 预计算特殊函数值 [psi, chi] = riccati_bessel(x, n_max); % 使用GPU加速（如有） if gpuDeviceCount > 0 x = gpuArray(x); m = gpuArray(m); end

6. 进阶应用示例

掌握了基础计算后，我们可以扩展更多实用场景：

多分散体系计算：

% 定义粒径分布（对数正态分布） d_mean = 100e-9; % 100nm sigma = 0.2; d_values = logspace(log10(d_mean)-3*sigma, log10(d_mean)+3*sigma, 100); pdf_values = lognpdf(d_values, log(d_mean), sigma); % 加权平均计算 Qext_total = 0; for i = 1:length(d_values) params.d = d_values(i); [Qext, ~, ~] = mie_calc(params); Qext_total = Qext_total + Qext * pdf_values(i); end Qext_total = Qext_total / sum(pdf_values);

角度分布计算：

function S = mie_scattering_angle(theta, x, m, a_n, b_n) % theta: 散射角度数组(弧度) % 计算角分布函数 mu = cos(theta); S1 = zeros(size(theta)); S2 = zeros(size(theta)); for n = 1:length(a_n) [pi_n, tau_n] = legendre_deriv(n, mu); S1 = S1 + (2*n+1)/(n*(n+1)) * (a_n(n)*pi_n + b_n(n)*tau_n); S2 = S2 + (2*n+1)/(n*(n+1)) * (a_n(n)*tau_n + b_n(n)*pi_n); end S = 0.5*(abs(S1).^2 + abs(S2).^2); end

与实验数据对比：

% 加载实验数据 exp_data = readmatrix('scattering_data.csv'); exp_x = exp_data(:,1); exp_Qext = exp_data(:,2); % 计算理论值 theo_Qext = arrayfun(@(x) mie_calc(set_params(x)), exp_x); % 绘制对比图 figure; semilogx(exp_x, exp_Qext, 'o', 'DisplayName', '实验数据'); hold on; semilogx(exp_x, theo_Qext, '-', 'DisplayName', '理论计算'); xlabel('尺寸参数 x'); ylabel('Q_{ext}'); legend; grid on;

7. 完整代码整合

将所有功能整合为可重用的Matlab类：

classdef MieCalculator properties wavelength % 入射波长[m] diameter % 粒子直径[m] m_particle % 粒子复折射率 n_medium % 介质折射率 n_max % 级数截断项数 end methods function obj = MieCalculator(wl, d, m, n_med) % 构造函数 obj.wavelength = wl; obj.diameter = d; obj.m_particle = m; obj.n_medium = n_med; obj.n_max = obj.calculate_nmax(); end function nmax = calculate_nmax(obj) % 计算建议的级数截断项数 x = obj.size_parameter(); nmax = round(2 + 4*x^(1/3) + x); end function x = size_parameter(obj) % 计算尺寸参数 x = pi * obj.diameter * obj.n_medium / obj.wavelength; end function [Qext, Qsca, Qabs] = calculate_efficiencies(obj) % 计算效率因子 [a_n, b_n] = obj.calculate_coefficients(); n = 1:obj.n_max; Qsca = (2/obj.size_parameter^2) * sum((2*n+1).*(abs(a_n).^2 + abs(b_n).^2)); Qext = (2/obj.size_parameter^2) * sum((2*n+1).*real(a_n + b_n)); Qabs = Qext - Qsca; end function [a_n, b_n] = calculate_coefficients(obj) % 计算散射系数 x = obj.size_parameter(); m_rel = obj.m_particle / obj.n_medium; a_n = zeros(1, obj.n_max); b_n = zeros(1, obj.n_max); for n = 1:obj.n_max [a_n(n), b_n(n)] = obj.mie_coefficient(n, x, m_rel); end end function S = scattering_matrix(obj, theta) % 计算散射矩阵元素 [a_n, b_n] = obj.calculate_coefficients(); S1 = zeros(size(theta)); S2 = zeros(size(theta)); for n = 1:obj.n_max [pi_n, tau_n] = obj.legendre_deriv(n, cos(theta)); S1 = S1 + (2*n+1)/(n*(n+1)) * (a_n(n)*pi_n + b_n(n)*tau_n); S2 = S2 + (2*n+1)/(n*(n+1)) * (a_n(n)*tau_n + b_n(n)*pi_n); end S = [S1; S2]; % 返回S1和S2矩阵 end end methods (Access = private) function [a_n, b_n] = mie_coefficient(~, n, x, m) % 私有方法：计算单个散射系数 psi_n = sqrt(pi*x/2) * besselj(n+0.5, x); psi_n1 = sqrt(pi*x/2) * besselj(n-0.5, x); chi_n = sqrt(pi*x/2) * bessely(n+0.5, x); chi_n1 = sqrt(pi*x/2) * bessely(n-0.5, x); D_n = (m*psi_n1 - (n/x)*psi_n) / (m*psi_n - (n/x)*psi_n1); G_n = (psi_n1 - (n/x)*psi_n) / (psi_n - (n/x)*psi_n1); a_n = (D_n/m + n/x)*psi_n - psi_n1; a_n = a_n / ((D_n/m + n/x)*psi_n - psi_n1 + 1i*((D_n/m + n/x)*chi_n - chi_n1)); b_n = (m*G_n + n/x)*psi_n - psi_n1; b_n = b_n / ((m*G_n + n/x)*psi_n - psi_n1 + 1i*((m*G_n + n/x)*chi_n - chi_n1)); end function [pi_n, tau_n] = legendre_deriv(~, n, mu) % 私有方法：计算连带Legendre函数及其导数 P_n = legendreP(n, mu); P_n1 = legendreP(n-1, mu); pi_n = P_n; tau_n = n * mu .* P_n - (n + 1) * P_n1; end end end

使用这个类可以简化后续计算：

% 初始化计算器 mie = MieCalculator(532e-9, 100e-9, 1.59+0.01i, 1.33); % 计算效率因子 [Qext, Qsca, Qabs] = mie.calculate_efficiencies(); % 计算角分布 theta = linspace(0, pi, 180); S = mie.scattering_matrix(theta); % 绘制角分布 figure; polarplot(theta, 10*log10(abs(S(1,:)).^2), 'LineWidth', 2); title('米氏散射角分布 (S11)');