当前位置：首页 > news >正文

用Matlab复现普朗克黑体辐射曲线：从公式推导到一键出图的保姆级教程

news 2026/6/4 5:10:55

用Matlab复现普朗克黑体辐射曲线：从公式推导到一键出图的保姆级教程

在热辐射研究中，普朗克黑体辐射定律是理解物体电磁辐射特性的基石。对于物理、光学或遥感领域的研究者而言，能够准确绘制不同温度下的黑体辐射曲线，不仅是验证理论的重要手段，也是分析实验数据的基础技能。本文将手把手带你从公式推导开始，逐步实现Matlab代码编写，最终生成可直接用于学术发表的精美图表。

1. 普朗克公式的物理内涵与数学表达

黑体辐射的普朗克定律描述了理想黑体在热平衡状态下，单位表面积向半球空间发射的光谱辐射出射度（spectral radiant exitance）与波长λ、温度T的关系：

$$ M_{\lambda}(T) = \frac{c_1}{\lambda^5} \frac{1}{e^{c_2/(\lambda T)} - 1} $$

式中两个关键常数需要特别注意单位转换：

第一辐射常数$c_1 = 2\pi h c^2 = 3.7418 \times 10^4$ W·μm⁴/cm²
第二辐射常数$c_2 = hc/k = 1.4388 \times 10^4$ μm·K

注意：在Matlab实现时，必须确保所有物理量的单位统一。常见错误是忽略波长单位（微米）与常数单位的匹配。

理解公式中每个参数的物理意义至关重要：

参数	物理意义	典型单位	Matlab变量名
λ	辐射波长	微米(μm)	`lambda`
T	黑体温度	开尔文(K)	`T`
M_λ	光谱辐射出射度	W/(cm²·μm)	`M_lambda`

2. Matlab实现的核心代码解析

2.1 普朗克公式的函数封装

创建一个独立的函数文件planck_law.m，实现公式计算：

function M_lambda = planck_lambda(lambda, T) % 计算普朗克黑体辐射定律 % 输入: % lambda - 波长数组(μm) % T - 温度标量或数组(K) % 输出: % M_lambda - 光谱辐射出射度(W/cm²/μm) c1 = 3.7418e4; % 第一辐射常数 (W·μm⁴/cm²) c2 = 1.4388e4; % 第二辐射常数 (μm·K) % 处理除零问题 exp_term = exp(c2./(lambda.*T)); M_lambda = c1./( (lambda.^5) .* (exp_term - 1) ); % 处理可能出现的NaN值（当lambda=0时） M_lambda(isnan(M_lambda)) = 0; end

代码优化技巧：

使用.*和./实现数组运算，避免循环
添加isnan检查增强鲁棒性
清晰的输入输出注释方便后续调用

2.2 主程序架构设计

主程序需要完成以下任务：

定义温度范围和波长范围
调用普朗克函数计算辐射曲线
绘制多温度对比曲线
标记峰值波长位置
美化图形输出

%% 黑体辐射曲线绘制主程序 clear; close all; clc; % 1. 参数设置 temps = [300, 500, 800, 1200, 2000, 3000]; % 温度数组(K) lambda = logspace(-1, 2, 500)'; % 波长范围0.1-100μm，对数均匀分布 % 2. 预分配数组并计算 M = zeros(length(lambda), length(temps)); peak_lambda = zeros(size(temps)); peak_M = zeros(size(temps)); for i = 1:length(temps) M(:,i) = planck_lambda(lambda, temps(i)); % 寻找峰值位置 [peak_M(i), idx] = max(M(:,i)); peak_lambda(i) = lambda(idx); end % 3. 绘图设置 figure('Position', [100, 100, 800, 600]); colors = lines(length(temps)); % 使用 distinguishable colors

3. 专业级图形绘制技巧

3.1 多曲线对数坐标绘制

使用loglog命令实现双对数坐标显示：

% 绘制辐射曲线 hold on; for i = 1:length(temps) h(i) = loglog(lambda, M(:,i), 'LineWidth', 1.8, ... 'Color', colors(i,:), ... 'DisplayName', sprintf('%d K', temps(i))); % 标记峰值点 stem(peak_lambda(i), peak_M(i), '--', ... 'Color', colors(i,:), 'LineWidth', 1.2, ... 'MarkerFaceColor', colors(i,:)); % 添加温度标注 text(peak_lambda(i)*1.5, peak_M(i)*0.7, ... sprintf('(%.1f μm, %.2e)', peak_lambda(i), peak_M(i)), ... 'Color', colors(i,:), 'FontSize', 9); end hold off;

3.2 图形美化关键参数

一张可直接用于发表的图形需要精心调整以下元素：

% 坐标轴与标题设置 set(gca, 'XScale', 'log', 'YScale', 'log'); xlim([0.1 100]); ylim([1e-4 1e6]); grid on; box on; % 标签设置 xlabel('波长 \lambda (\mum)', 'FontSize', 12, 'FontWeight', 'bold'); ylabel('光谱辐射出射度 M_{\lambda} (W·cm^{-2}·\mum^{-1})', ... 'FontSize', 12, 'FontWeight', 'bold'); title('不同温度下黑体的光谱辐射曲线', ... 'FontSize', 14, 'FontWeight', 'bold'); % 图例设置 legend(h, 'Location', 'northeastoutside', 'FontSize', 10); % 导出设置 set(gcf, 'Color', 'w'); exportgraphics(gcf, 'blackbody_radiation.png', 'Resolution', 300);

常见问题解决：

曲线显示不完整？调整xlim和ylim范围
峰值标记重叠？使用*1.5和*0.7等系数偏移文本位置
颜色区分不明显？尝试parula、jet或其他colormap

4. 维恩位移定律的验证与应用

普朗克曲线峰值波长位置遵循维恩位移定律：

$$ \lambda_{max} T = b \quad (b \approx 2898 \ \mu m \cdot K) $$

我们可以在Matlab中验证这一关系：

%% 维恩位移定律验证 wien_constant = peak_lambda .* temps'; figure; plot(temps, peak_lambda, 'o-', 'LineWidth', 2); xlabel('温度 (K)'); ylabel('峰值波长 (\mum)'); title('维恩位移定律验证'); grid on; % 显示平均常数 avg_b = mean(wien_constant); text(mean(temps), mean(peak_lambda), ... sprintf('λ_{max}·T ≈ %.1f μm·K', avg_b), ... 'FontSize', 12, 'HorizontalAlignment', 'center');

实际应用技巧：

红外测温：通过测量辐射峰值波长反推物体温度
恒星分类：分析恒星光谱确定表面温度
热成像优化：根据目标温度范围选择最佳探测器波段

5. 高级功能扩展

5.1 交互式温度调节GUI

创建可交互调节的图形界面：

function interactive_blackbody() % 创建图形窗口 fig = figure('Position', [200, 200, 850, 650]); % 添加滑动条 uicontrol('Style', 'slider', ... 'Min', 300, 'Max', 6000, 'Value', 3000, ... 'Position', [100, 20, 600, 20], ... 'Callback', @updatePlot, ... 'Tag', 'tempSlider'); % 添加温度显示 uicontrol('Style', 'text', ... 'Position', [350, 45, 100, 20], ... 'Tag', 'tempText'); % 初始绘图 lambda = logspace(-1, 2, 500)'; updatePlot(); function updatePlot(~,~) temp = get(findobj('Tag', 'tempSlider'), 'Value'); set(findobj('Tag', 'tempText'), 'String', sprintf('温度: %d K', round(temp))); % 计算辐射曲线 M = planck_lambda(lambda, temp); % 更新图形 if isempty(findobj('Tag', 'blackbodyLine')) loglog(lambda, M, 'r-', 'LineWidth', 2, 'Tag', 'blackbodyLine'); xlabel('波长 (\mum)'); ylabel('光谱辐射出射度 (W·cm^{-2}·\mum^{-1})'); title(sprintf('黑体辐射曲线 (T = %d K)', round(temp))); grid on; xlim([0.1 100]); ylim([1e-4 1e6]); else set(findobj('Tag', 'blackbodyLine'), 'YData', M); title(sprintf('黑体辐射曲线 (T = %d K)', round(temp))); end end end

5.2 辐射总量计算

通过积分普朗克公式可以得到斯特藩-玻尔兹曼定律：

% 计算总辐射出射度（斯特藩-玻尔兹曼定律） sigma = 5.670374419e-12; % W/cm²/K⁴ total_radiation = sigma * temps.^4; % 数值积分验证 calculated_radiation = zeros(size(temps)); for i = 1:length(temps) calculated_radiation(i) = trapz(lambda, M(:,i)); end % 结果显示 disp('温度(K) 理论值(W/cm²) 计算值(W/cm²) 相对误差'); disp([temps', total_radiation', calculated_radiation', ... abs(total_radiation-calculated_radiation)./total_radiation*100]);