黑体辐射和普朗克曲线，这个东西在物理课本里出现频率极高。但说实话，我在本科那会儿学热力学的时候，看到普朗克公式那密密麻麻的指数项，第一反应其实是“这玩意儿背下来能干嘛”。直到后来做红外测温相关的模拟项目，需要自己画不同温度下的黑体辐射曲线、找峰值波长，我才意识到这个公式不仅不是死记硬背的纸面理论，还是很多工程判断的底层依据。这篇文章就把我用Matlab绘制普朗克曲线并定位峰值点的完整过程写出来，包括公式选型、单位陷阱、寻峰算法对比，以及一些网上教程很少提到的细节坑。适合正在学热辐射、准备做红外系统仿真、或者写课程作业需要画黑体辐射曲线的同学参考。
1. 项目核心：普朗克曲线与维恩位移定律
1.1 从黑体辐射公式到光谱辐射出射度
先明确一个概念：普朗克曲线描述的可不是简简单单“温度越高辐射越强”这种定性关系，它给出的是黑体在特定温度T下，辐射能量在不同波长上的定量分布。常用的波长形式公式是这样的：
B_λ(T) = (2hc² / λ⁵) × 1 / (exp(hc / (λk_BT)) - 1)
其中h是普朗克常数，c是真空光速，k_B是玻尔兹曼常数，T是热力学温度，λ是波长。这个式子的物理意义可以理解为：在温度T下，单位面积的黑体表面、单位时间、单位立体角内，波长落在λ附近单位波长间隔中的辐射能量。
这里有个特别容易让初学者犯迷糊的地方：公式里的变量是波长λ，但实际使用时，λ的取值跨度非常大。常规情况下我们画曲线，横轴要么取纳米（可见光到近红外），要么取微米（热红外波段）。单位不一致会直接影响后面寻峰代码的结果，这部分我在第2章里单独说。
1.2 为什么要找峰值点：峰值波长背后的工程意义
找到峰值点不是单纯为了完成一个“找极值”的编程题。峰值波长λ_max与温度T之间满足维恩位移定律：λ_max = b / T，其中b约等于2898 μm·K。这个关系在工程上的应用非常广：
- 红外测温仪选型：测常温物体（约300K）时，辐射峰值在9.7μm附近，属于长波红外，那就需要选择相应响应波段的探测器，而不是随便拿一个可见光传感器去测。
- 太阳光谱分析：太阳表面温度约5800K，峰值波长大概在0.5μm，这正好解释了为什么可见光波段是太阳辐射能量最集中的区域。
- 照明与显示设计：色温这个概念本质上就和黑体辐射峰值位置强相关。暖色调光源的峰值偏长波，冷色调光源的峰值偏短波，背后的物理逻辑就是维恩位移定律。
所以，能把普朗克曲线画出来，并且能够准确定位峰值的坐标，本身就等于掌握了一个跨学科的分析工具。
1.3 整体实现思路
基于Matlab的实现路线可以拆成四个关键环节：
- 构建波长数组和物理常数，注意单位统一。
- 根据普朗克公式计算对应温度下的辐射出射度数组。
- 用plot绘图，观察曲线形状和峰值的大致位置。
- 用不同的数值方法求峰值坐标，并与维恩位移定律的理论值互相验证。
我实际操作下来，最耗时反而不是写代码，而是处理“看起来很简单但结果对不上”的问题，比如波长单位换算错误、指数项溢出、寻峰函数受到数组离散化影响导致定位不准等。后面我逐条讲清楚。
2. 任务拆解与工具选择细节
2.1 科学计算的常量设置与初始条件
Matlab本身并没有内置“普朗克公式”这样一个现成函数，但科学计算的优势在于公式可以自己写，关键是常量别抄错。我用的是CODATA推荐值：
- h = 6.62607015e-34 J·s
- c = 299792458 m/s
- k_B = 1.380649e-23 J/K
这里我特别强调一下：所有常量使用国际标准单位，波长也建议先统一成米，等绘图时再换算显示单位，这样可以避免公式中指数项的单位错乱。很多新手一上来就用“微米”代入公式，觉得数值写起来方便，结果算出来曲线峰值位置差了好几个数量级。原因很简单：公式里波长是以米为基准的，你若代入以微米为单位的数值，hc/(λkT)这个量就会凭空多出10^6倍的误差。
2.2 波长范围的选取原则
选横轴范围也是有讲究的。300K的常温黑体，峰值在9.7μm，如果你把横轴范围设定成0.4μm到0.7μm的可见光区，画出来的曲线就是单调递减，根本看不到峰。这不是公式错了，而是窗口选错了。
我的经验是：先根据目标温度T用维恩位移定律粗算一个λ_max，然后让横轴覆盖 0.1×λ_max 到 10×λ_max，这样既能观察到完整峰值，又能看到短波端和长波端各自的变化趋势。比如T=3000K时，λ_max约0.966μm，那我通常取λ范围从0.05μm到10μm，用对数坐标或线性坐标都行。文章里的图我一般线性坐标看峰值，对数坐标看拖尾。
2.3 绘图与寻峰的关键代码实现
下面给出一段我实际调试通过的完整Matlab代码。它包含了三个功能：计算普朗克分布、绘制曲线、用不同方法寻找峰值点。
% 普朗克黑体辐射曲线绘制与峰值点寻找 % 物理常量（国际单位制） h = 6.62607015e-34; % 普朗克常数 J·s c = 299792458; % 光速 m/s kB = 1.380649e-23; % 玻尔兹曼常数 J/K % 参数设定 T = 5800; % 温度，单位 K，这里用太阳表面温度 lambda_min = 1e-9; % 波长下限，1 nm lambda_max = 3e-6; % 波长上限，3 μm N = 50000; % 采样点数，越多峰值定位越准 lambda = linspace(lambda_min, lambda_max, N); % 普朗克公式（波长形式） % 注意使用点乘、点除 exp_arg = h * c ./ (lambda * kB * T); % 防止 exp 溢出，可以先判断 exp_arg(exp_arg > 700) = 700; % 超过阈值时，exp(result) 在双精度下会溢出 I = (2 * h * c^2 ./ lambda.^5) ./ (exp(exp_arg) - 1); % 绘图：横轴换算成微米，更方便观察 lambda_um = lambda * 1e6; figure; plot(lambda_um, I, 'b-', 'LineWidth', 1.5); xlabel('波长 (μm)'); ylabel('光谱辐射出射度 B_λ (W·sr^{-1}·m^{-3})'); title(['普朗克曲线 T = ', num2str(T), ' K']); grid on; % 方法1：用 max 函数找离散数组最大值 [maxI, idx] = max(I); peak_lambda_1 = lambda(idx); % 单位：米 peak_lambda_1_um = peak_lambda_1 * 1e6; fprintf('方法1（max函数）峰值波长: %.4f μm\n', peak_lambda_1_um); % 方法2：用 findpeaks 函数（需要 Signal Processing Toolbox） % 注意先用 small window 平滑，避免局部抖动干扰 Ismooth = smoothdata(I, 'gaussian', 200); [pks, locs] = findpeaks(Ismooth, lambda_um, 'SortStr', 'descend'); if ~isempty(pks) peak_lambda_2_um = locs(1); fprintf('方法2（findpeaks）峰值波长: %.4f μm\n', peak_lambda_2_um); else fprintf('方法2：在给定区间内未找到峰值\n'); end % 方法3：利用维恩位移定律的数值优化（fminbnd） % 因为求极大值等价于求负函数极小值 negI = @(lam) -((2 * h * c^2 / lam^5) / (exp(h * c / (lam * kB * T)) - 1)); [peak_lambda_3, ~] = fminbnd(negI, 0.1e-6, 2e-6); % 搜索区间需包含真实峰值 peak_lambda_3_um = peak_lambda_3 * 1e6; fprintf('方法3（fminbnd）峰值波长: %.4f μm\n', peak_lambda_3_um); % 理论对照：维恩位移定律 b = 2.897771955e-3; % 单位 m·K peak_theory_um = b / T * 1e6; fprintf('维恩位移定律理论峰值波长: %.4f μm\n', peak_theory_um);这段代码直接复制就能跑出结果。以T=5800K为例，三种方法得到的峰值基本上都在0.4995μm附近，与理论值2898/5800≈0.4997μm符合得很好。
2.4 为什么我建议用波长形式而不是频率形式
普朗克公式除了以波长为变量B_λ(T)，还有以频率为变量的形式B_ν(T)。两者在数学上是等价的，但峰值位置并不相同。这一点非常搞人：B_λ(T)的峰值波长对应维恩位移定律的λ_max，但如果你用频率形式去求峰值频率，再换算成波长，得到的数值和λ_max是不一致的。
根本原因在于：B_λ(T)表示单位波长间隔内的能量，B_ν(T)表示单位频率间隔内的能量。λ和ν是倒数关系，导致同一个光谱在不同横轴度量下呈现出不同的“形状”。做天文学或遥感数据处理时，经常要先确认自己用的是波长谱还是频率谱，否则对比文献数据就会闹笑话。我的建议是：工程应用中，优先用波长形式，因为维恩位移定律的经典表述和大多数探测器的响应波段都是按波长来定义的。
3. 实操过程与核心环节实现
3.1 参数的物理意义与敏感性分析
代码中的T是唯一需要频繁改动的物理量。很多朋友拿到代码后会先试一个温度，比如300K、3000K、5800K，然后对比曲线。这种操作没有问题，但要注意T对曲线形状的影响不仅仅是“峰值位置移动”：
- 温度升高，峰值波长变短，也就是蓝移。
- 温度升高，峰值辐射出射度显著增加，而且是按T^4增长，这就是斯特藩-玻尔兹曼定律的含义。
- 温度升高，短波段的辐射占比明显抬升，这就是为什么高温物体会从暗红变成橙黄再到白亮。
为了直观展示，我习惯在同一坐标系画多条温度曲线，并用归一化处理——把每个温度下的I都除以该温度下的峰值I_max。归一化后，能明显看到曲线峰值左移、形状“变瘦”，这比光看数值更有教学冲击力。
figure; Tlist = [300, 1000, 3000, 5800]; colors = lines(length(Tlist)); for i = 1:length(Tlist) T = Tlist(i); I_temp = (2 * h * c^2 ./ lambda.^5) ./ (exp(h * c ./ (lambda * kB * T)) - 1); normI = I_temp / max(I_temp); semilogx(lambda_um, normI, 'Color', colors(i, :), 'LineWidth', 1.5); hold on; end hold off; xlabel('波长 (μm)'); ylabel('归一化辐射出射度'); legend(arrayfun(@(x) sprintf('%d K', x), Tlist, 'UniformOutput', false)); grid on;这里用了semilogx对数横轴，因为300K曲线的峰值在9.7μm附近，若用线性轴，它和5800K曲线几乎挤成一团，观察效果很差。
3.2 数值寻峰的边界条件与搜索区间
虽然max函数是最简单粗暴的求峰值方式，但它的精度受采样点数N限制。假设我取λ范围从1nm到3μm，采样50000个点，每个点间距约60pm，那么这个分辨率对峰值定位的误差大约在几十皮米量级，肉眼完全够用。但如果想追求更高精度，或者要求程序自动返回一个平滑的峰值坐标，fminbnd这类连续优化方法更合适。
使用fminbnd时，必须注意搜索区间要包裹住真实峰值，而且不能太大。我在代码里写的[0.1e-6, 2e-6]是针对5800K的太阳光谱的合理区间。如果你把温度改成300K，峰值跑到9.7μm，那么搜索区间必须相应调整，否则fminbnd会在边界处截断并给出错误结果。这是新手最容易踩的坑——换了一个温度，却没换搜索区间。
为了避免这种“埋参数”的问题，我实际操作时会先动态计算一个粗略峰值位置，然后以它为中心构造搜索区间：
lambda_guess = b / T; % 维恩定律估算峰值 search_low = lambda_guess * 0.5; search_high = lambda_guess * 1.5; [peak_lambda, ~] = fminbnd(negI, search_low, search_high);这样无论温度怎么变，搜索区间都会自动跟随，不需要每次手改。
3.3 如何验证寻峰结果是否可信
数值找峰之后别急着往下走，拿理论值对照是最快、也最稳妥的验证手段。维恩位移定律的常数b = 2898 μm·K，那么理论峰值波长就是b/T。我在前面代码里已经把这部分写好了，运行结果基本能在小数点后三位吻合。
如果出现明显偏差，问题往往出在这几个地方：
- 波长单位在计算过程中混用，导致指数项数量级错乱。
- 采样点太少，导致离散max值的索引位置偏移。
- 寻峰区间设置不当，fminbnd收敛到边界。
- 没有考虑普朗克公式中分母的-1项，误用瑞利-金斯近似。
其中瑞利-金斯近似导致的错误比较隐蔽，在长波范围（hc/λkT << 1）时，exp(x) - 1约等于x，公式退化成2ckT/λ⁴，这个近似在长波端没问题，但用它来找峰值就完全偏离了。所以代码里必须保留完整的括号项。
4. 常见问题与排查技巧实录
4.1 指数项溢出导致曲线出现NaN或Inf
当λ特别小或者T特别低时，hc/(λkT)这个数会非常大。Matlab中exp(700)以上就已经接近double类型的表示上限，再往上就是Inf，导致I的计算出现Inf或NaN。
解决办法就是在计算前先做一次钳位，把所有超过700的指数值强行设成700。从物理角度看，这些短波区的辐射本来就低到可以忽略，钳位不会影响可视化效果，也不影响峰值定位，因为峰值区不会落在这些极端位置上。
exp_arg = h * c ./ (lambda * kB * T); exp_arg(exp_arg > 700) = 700;4.2 峰顶“抖动”导致findpeaks找错点
如果你用的原始数据没有经过平滑，而采样点数又特别多，曲线在峰顶附近可能有极其微小的数值波动。findpeaks函数可能会把这些局部波动识别成许多个小峰值，返回的locs就有多个，排序后还不一定第一个就是全局最大。
我实际踩过的坑是：在太阳光谱附近取50000个点，findpeaks返回了七八个“峰”，主峰反而排在第二三个。最后我加了一步高斯平滑，窗口选200左右，再调findpeaks，问题就解决了。注意：平滑窗口不能太大，否则会把峰顶真正削平，导致峰值波长系统性地偏移几个纳米。
另外，如果Matlab环境没有Signal Processing Toolbox，直接用max会简单得多，精度也能接受。并不是所有场景都必须上findpeaks。
4.3 横轴单位带来的标签和坐标误解
很多朋友第一次画普朗克曲线，写完代码后发现曲线尖峰在某个奇怪的位置，或者发现y轴数值大得离谱。我排查过几次，基本都是因为横轴用了微米、y轴却还保留着基于米单位的辐射出射度，导致坐标标签和数值完全对不上。
建议养成一个习惯：计算全部用国际单位制，绘图前统一换算显示单位。横轴要显示微米，就乘1e6；要显示纳米，就乘1e9。标签上明确写出单位，避免别人看你的图时还得猜。
4.4 维恩位移定律常数有时候记不准怎么办
如果你不想硬背b=2898μm·K，有个很实用的方法：直接在代码里用h、c、kB自己算。维恩常数实际上来自对普朗克公式求导等于零的数值解，理论表达式大约为：
b ≈ hc/(4.9651kB)
这意味着你完全可以不依赖那个常数，而是在代码里定义：
b_calc = h * c / (4.965114 * kB);这样不仅能减少记忆负担，数值也能对得上。
4.5 用matlab的App Designer或者实时脚本能操作吗
如果你不喜欢写脚本，想在交互环境里拖动滑块观察温度对曲线的影响，可以用Matlab的实时脚本（Live Script）或者App Designer做个简单界面。核心思路是：一个Slider控件控制T的值，回调函数里重新计算普朗克公式并更新plot数据。这个功能我实际做过，效果很好，特别适合用来给学生演示“温度升高，峰值蓝移”的直观感受。
% 在App Designer中 Slider 的 ValueChangedFcn 回调示例 T = app.Slider.Value; I = (2*h*c^2 ./ lambda.^5) ./ (exp(h*c ./ (lambda*kB*T)) - 1); plot(app.UIAxes, lambda_um, I);难点在于确保app中的lambda数组和常量在所有回调中都能访问，一般通过组件属性或在整个app的startup回调里初始化公共变量。
5. 从曲线到自动找峰：封装成一个可复用函数
5.1 为什么建议封装成函数
如果你只是临时画一次曲线，脚本完全够用。但如果你和我一样，需要在多个仿真任务中反复计算不同温度下的峰值波长，甚至要把峰值信息喂给其他程序，那么封装成函数就是必要的。
我给出的方案是这样：输入参数为温度T和可选的波长范围，返回值为峰值波长、峰值辐射出射度和绘图的句柄。
function [peak_wave_um, peak_value, figHandle] = planck_curve_peak(T, varargin) % 输入：T - 温度(K) % 可选参数：波长范围 [lam_min, lam_max]（米） % 输出：峰值波长(μm)，峰值辐射出射度，图形句柄 h = 6.62607015e-34; c = 299792458; kB = 1.380649e-23; if nargin < 2 lambda_min = 0.1e-9; lambda_max = 5e-6; else lambda_min = varargin{1}(1); lambda_max = varargin{1}(2); end lambda = linspace(lambda_min, lambda_max, 50000); exp_arg = h*c ./ (lambda*kB*T); exp_arg(exp_arg > 700) = 700; I = (2*h*c^2 ./ lambda.^5) ./ (exp(exp_arg) - 1); [peak_value, idx] = max(I); peak_wave_um = lambda(idx) * 1e6; figHandle = figure; plot(lambda*1e6, I, 'r-', 'LineWidth', 1.5); xlabel('波长 (μm)'); ylabel('光谱辐射出射度'); grid on; end封装完成之后，调用一行就出结果：
[peak_6000, val_6000] = planck_curve_peak(6000);然后我还可以连续算几个温度，把峰值点标在一张图上，自动观察维恩位移趋势。
5.2 后续扩展思路
很多做光学设计的同行还会继续做黑体辐射积分，比如计算特定波段内的辐射功率比例。这时候就需要对普朗克公式做数值积分，Matlab的integral函数可以直接上手。甚至可以把积分比例写成温度的函数画出来，看不同温度下可见光波段占比的变化。
还有一个比较实用的小扩展：把峰值位置标在图上。用plot函数在峰值点画一个圆圈标记，再配合text标注秒出坐标值，这样图会显得更专业：
hold on; plot(peak_wave_um, peak_value, 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 8); text(peak_wave_um*1.05, peak_value*0.9, sprintf('%.2f μm', peak_wave_um));5.3 这件事做完之后还能用在哪
画普朗克曲线这件事看着很基础，但它是很多实际项目里的底层模块。比如我在做光伏电池材料吸收谱匹配时，就需要先算特定温度黑体光源的峰值波段，判断应该重点关注哪段波长的量子效率。学材料或者测控的同学，往后做红外探测器性能仿真、大气辐射传输估算，也避不开这段公式。
摸透这一个函数，等于把“物理公式转成数值代码”、“数值寻峰”、“单位换算”、“结果验证”四个能力全部练了一遍。如果能顺着这个思路再做一版用Python实现的，那以后就算换语言也不用慌。
最后再分享一个实用技巧：如果只是想快速体验一下峰值点和温度之间的变化规律，可以在Matlab命令行里直接执行下面这段极简代码，体会维恩位移从微妙到显著的过程：
for T = [300, 600, 1000, 2000, 3000, 6000] lam = linspace(0.05e-6, 20e-6, 100000); I = (2*h*c^2 ./ lam.^5) ./ (exp(h*c ./ (lam*kB*T)) - 1); [~, idx] = max(I); fprintf('T=%5d K, peak=%.3f μm\n', T, lam(idx)*1e6); end这个循环可以帮助你快速建立“温度→峰值位置”的直观手感，比单纯背诵维恩位移定律的常数牢靠得多。我当年把这些数值算过一遍之后，再看红外测温仪的参数表，心里基本就有底了。
