企业官网建设流程全解析

一提到普朗克曲线，很多人第一反应是物理课本上那条“两头低、中间高”的曲线，觉得离代码很远。其实用 Matlab 把普朗克公式画出来，再把峰值点精确找出来，是热辐射分析、色温估算、光电检测甚至 LED 光谱评估里非常常用的一套小流程。我见过不少人在这一步反复折腾：要么曲线画出来全是 0，要么峰值点根本不在预期的位置，问题大多出在单位制、数值范围和求峰值的策略上。

这篇文章我就以“用 Matlab 绘制普朗克曲线并寻找峰值点”为主线，把从物理模型到代码实现、从粗略峰值到高精度峰值提取的完整过程拆开讲。适合刚接触 Matlab 的物理、光电类专业学生，也适合需要快速做热辐射估算的工程师。你不需要懂很深的光谱学，只要知道公式长什么样、单位怎么换算，跟着步骤走就能把图做出来。

1. 画图前必须想清楚的两件事：公式版本和单位制

很多人一上来就写B = 2*h*c^2 ./ lam.^5 ./ (exp(...)-1)，逻辑没错，但画出来经常是一条水平线或者一片“NaN 警告”。问题往往不在公式本身，而在你选的公式版本和单位之间的换算。

1.1 普朗克公式的两个版本，差别不只是变量

普朗克公式按自变量不同有两个常见写法。一种是按波长表示，单位是 W·m⁻³·sr⁻¹；另一种是按频率表示，单位是 W·m⁻²·Hz⁻¹·sr⁻¹。如果你手上的数据是波长坐标，直接用波长版本最顺手：

B(λ, T) = (2·h·c² / λ⁵) · 1 / (exp(h·c / (λ·k·T)) − 1)

其中 h=6.626e-34 J·s 是普朗克常数，c=3.0e8 m/s 是光速，k=1.381e-23 J/K 是玻尔兹曼常数，λ 是波长，T 是绝对温度。

如果改用频率版本，按 ν = c/λ 换过去也能得到同样的物理量，但峰值点位置和曲线形状会落在完全不同的坐标区间，很容易把自己绕晕。我的建议很直接：做热辐射可视化和峰值分析，固定用波长版本，不要混合用。

1.2 单位换算决定生死：纳米、微米和米

Matlab 的数值计算默认是国际单位制，如果你直接写lam = 500，那等于设定了波长 500 米，算出来的结果毫无意义。我做首次验证时，习惯把所有长度统一成米，最后画图时再手动除以 1e-6 转成微米显示，刻度数字好看得多。

比如一个常见场景：太阳表面温度约 5800 K，辐射峰值应该在可见光附近。按维恩位移定律估算峰值约 0.5 μm，也就是 500 nm。如果我用linspace(0.1e-6, 3e-6, 10000)生成 0.1 μm 到 3 μm 的网格，坐标单位是米，横轴标签写成Wavelength (μm)时就要把数组除以 1e-6。不要图省事把波长直接写 500，那样后续所有计算都会出问题。

1.3 归一化是为了让曲线能“看”

普朗克公式算出来的 B 值是非常大的数还是极小的数，取决于波长区间和温度。例如近红外区间、几千开尔文温度下，数值可能在 1e13 这个量级；如果只看相对光谱形状，谁大谁小无所谓，但直接绘图时容易出现“全是一条直线”的观感，因为最大值太突出，细节全被压扁了。

我常用的处理是把曲线除以它自己的最大值，把纵轴压到 0 到 1 之间，也就是所谓“归一化普朗克曲线”。这样一条曲线画出来，你一眼就能看出峰值在哪、下降趋势是什么形状，也能把不同温度的多条曲线放在同一张图里对比，不会因为数值量级差异导致高温曲线把低温曲线压成一条横线。

2. 一版能直接用的 Matlab 绘图代码

确认了模型和单位之后，代码本身并不复杂。但为了让不同基础的人都能跑通，我把代码分成“基础绘图”和“局部细化”两层，先保证能出图，再做峰值提取。

2.1 最小可用版本：10 行出图

这是一个我经常发给同事的模板：

clear; clc; h = 6.626e-34; % 普朗克常数，J·s c = 3.0e8; % 光速，m/s k = 1.381e-23; % 玻尔兹曼常数，J/K T = 5800; % 黑体温度，K lam = linspace(0.1e-6, 3e-6, 10000); % 波长范围，单位米 B = (2*h*c^2 ./ lam.^5) ./ (exp(h*c ./ (lam*k*T)) - 1); B_norm = B / max(B); plot(lam*1e6, B_norm, 'LineWidth', 1.5); xlabel('Wavelength (\mum)', 'FontSize', 12); ylabel('Normalized Radiance', 'FontSize', 12); title(['Blackbody Radiation at T = ', num2str(T), ' K']); grid on;

这里我故意用了lam*1e6再画图，这样横轴直接显示微米，标签用\mum让 Matlab 渲染出希腊字母 μ。第一次跑通后，你会发现峰值大概在 0.5 附近，对应 500 nm 左右，和维恩位移定律的估算一致，这个“对得上”本身就是验证代码是否写对的第一步。

2.2 为什么用linspace而不是步长写法

有些教程里用0.1e-6:1e-8:3e-6这种步长写法，不是不行，但你要自己去数数组里到底有多少个点，不够直观。用linspace(起点, 终点, 点数)的好处是清晰：第 3 个参数直接控制采样密度，后续求峰值时也能知道横轴每个间隔的长度，方便做细化和误差分析。

采样点数我建议至少 5000 起步。原因很简单：往后你要做峰值提取，如果采样太稀疏，最大值索引对应的横坐标离真实峰值差得很远，误差可能大到几十纳米。10000 点在多数场景下是一个性价比很高的选择，算得很快，精度也够日常用了。

2.3 画图前先检查三个“雷区”

我踩过几个坑，提醒你绕开：

第一，exp函数溢出。普朗特公式的指数部分是 h·c/(λ·k·T)，当 λ 很小或者 T 很大的时候，指数参数可能超过 700，Matlab 的 double 会直接把它当成 Inf，结果整条曲线变成 NaN。我处理这个问题的办法，是先把指数参数算出来，看一眼最大值，超过 600 就提醒自己要考虑换波长起点或改用化简形式。

第二，max(B)如果 B 里有 NaN，那B/max(B)会得到一堆 NaN。写代码时最好在归一化前用any(isnan(B))自检一次，有 NaN 先解决区间问题，而不是硬着头皮画图。

第三，title字符串里用num2str(T)拼接，要注意空格。这个不影响计算，但会影响图面美观，属于经验问题。

3. 峰值点提取的三种方法，按精度递进

峰值点提取是这篇博文的核心。很多人画完图就直接用光标在图上点一下，觉得“差不多就行了”。如果只是看趋势当然无所谓，但如果你要拿这个波长去做测温换算、材料光谱匹配或者光学设计，就必须用代码精确算出来。下面是三种由易到难的做法。

3.1 方法一：暴力搜索，max加find

最简单的方法是直接从数组里找最大值和下标：

[B_max, idx] = max(B_norm); lam_peak = lam(idx);

% 显示结果 fprintf('峰值波长: %.4f μm\n', lam_peak*1e6);

这种做法的精度完全取决于采样间隔。比如一开始用linspace(0.1e-6, 3e-6, 10000)，采样间隔是 290 nm / 9999 ≈ 0.029 μm，也就是 29 纳米。直接用原始索引取峰值，误差最多可能有 29 nm 的一半左右，约 15 nm。对很多工程估算来说，这个误差不算小。尤其温度较低时曲线比较平坦，顶部区域变化缓慢，索引位置经常偏出真实峰值好几个点。

暴力搜索适合快速看看峰值在哪里，但不适合精度要求高的场景。如果你想提高暴力搜索的精度，简单手段是局部加密采样：在粗峰值附近再生成一个很小的精细区间，比如lam_fine = linspace(lam(idx-5), lam(idx+5), 100000)，重新算一遍 B 再找 max。这个思路很简单，但计算量会增加，适合临时使用。

3.2 方法二：解析求导，直接套维恩位移定律

如果你不想依赖采样密度，可以用数学方法直接求峰值。普朗克公式对 λ 求导后令导数为零，会得到一个超越方程，最终化简成：

x / (1 − e^(−x)) = 5

其中 x = h·c / (λ·k·T)。这个方程的解 x ≈ 4.965114。因此峰值波长满足：

λ_max = h·c / (4.965114·k·T)

把常数代入，h·c/k ≈ 1.4388e-2 m·K，因此 λ_max·T ≈ 2.8978e-3 m·K，这就是维恩位移定律。我其实不建议你记那么多位小数，记住“峰值波长乘温度约等于 2898 μm·K”就够了。

所以如果只是找峰值波长，直接用维恩定律是最高效的：

T = 5800; lam_peak = 2.8978e-3 / T; % 单位：米 fprintf('理论峰值波长: %.4f μm\n', lam_peak*1e6);

这种方法不依赖采样，还特别适合做理论估算时的交叉验证。刚做完方法一的同学，可以立刻和这个理论值对比，差异如果在几十纳米内，说明代码流程正确；如果差了几百纳米，通常是区间设置或者公式写错了。

3.3 方法三：局部多项式拟合，亚采样精度提取

如果既不想加大采样量，又想要比暴力搜索更准确，可以用局部拟合的思路。核心逻辑是：在最大值附近的几个采样点，用多项式拟合曲线形状，再对拟合多项式求导解根。因为普朗克曲线在峰值附近非常光滑，三到五次多项式就能拟合得很好。

操作步骤是：先用max找到峰值索引 idx，然后取 idx-2 到 idx+2 这 5 个点，做三次多项式拟合，得到系数 p，然后对多项式求导并求根，找到导数由正变负的根，这个根对应的横坐标就是更精细的峰值波长。

seg = idx-2:idx+2; lam_seg = lam(seg); B_seg = B(seg); p = polyfit(lam_seg, B_seg, 3); dp = [3*p(1), 2*p(2), p(3)]; % 导数多项式系数 r = roots(dp); % 筛选落在区间内的实根 r = r(abs(imag(r)) < 1e-10); r = r(real(r) >= lam_seg(1) & real(r) <= lam_seg(end)); lam_peak_fine = real(r(1)); fprintf('拟合峰值波长: %.6f μm\n', lam_peak_fine*1e6);

这里有两个容易被忽视的细节。一是polyfit在区间比较窄时，数值稳定性一般没问题，但如果区间跨越的波长范围太大，拟合反而会失真，所以只取峰值附近 5 到 7 个点比较合适。二是roots可能给出复数根，一定要增加筛选步骤，只保留落在采样区间里的实根。

我用太阳表面温度 5800 K 测试过，暴力搜索得到的峰值可能在 0.501 或 0.497 μm 附近跳变，而局部拟合法基本能稳定到 0.4999 μm 左右，和维恩定律的理论值 0.4996 μm 非常吻合。这个精度已经足够应对绝大多数科研和工程场景了。

4. 多条曲线对比与真实场景应用

单条曲线只是热身。实际场景里，你往往要对比不同温度下的黑体辐射，或者根据一个未知光谱反推温度。所以这节我把两种常见应用都展开说。

4.1 一图展示不同温度下的曲线移动

温度越低，峰值波长越长，这既是维恩位移定律的推论，也是“烧红的铁先变红、再变白”的物理本质。用 Matlab 把 3000 K、4000 K、5000 K、6000 K 四条曲线画一起，能直观看到峰值从左往右移动。

参考代码如下：

T_list = [3000, 4000, 5000, 6000]; colors = lines(numel(T_list)); legend_str = cell(1, numel(T_list)); figure; hold on; for i = 1:numel(T_list) T_i = T_list(i); B_i = (2*h*c^2 ./ lam.^5) ./ (exp(h*c ./ (lam*k*T_i)) - 1); B_i_norm = B_i / max(B_i); plot(lam*1e6, B_i_norm, 'LineWidth', 1.5, 'Color', colors(i, :)); [~, idx_i] = max(B_i_norm); plot(lam(idx_i)*1e6, B_i_norm(idx_i), 'o', 'MarkerSize', 6, ... 'MarkerFaceColor', colors(i, :), 'MarkerEdgeColor', 'k'); legend_str{i} = [num2str(T_i), ' K']; end hold off; xlabel('Wavelength (\mum)'); ylabel('Normalized Radiance'); legend(legend_str, 'Location', 'northeast'); grid on;

画好以后能清晰看到，3000 K 的曲线峰值已经接近 1 μm，而 6000 K 的曲线峰值明显移到 0.48 μm 附近。峰值点用圆圈标记出来，比纯看曲线方便得多。

这里有个容易被误解的物理事实：3000 K 的温度对白炽灯来说已经算很亮了，但它的峰值辐射在近红外，不在可见光区，所以我们看到的是暖黄色，而不是蓝色或白色。很多初学者觉得“温度越高越偏蓝”，其实就是维恩位移定律的体现：峰值波长和温度成反比。

4.2 实际数据反推温度的思路

如果你手里有一份实测的相对光谱强度数据，想知道这个辐射源等效黑体温度大概是多少，最简单的做法是：找到光谱数据的峰值波长，然后用维恩定律反推 T = 2898 / λ_max(μm)。

但要注意，实测数据往往包含仪器响应、介质吸收，峰值位置可能和真实黑体有偏差。更稳妥的方式是做一个最小二乘拟合：把普朗克公式归一化之后和实测归一化数据对比，扫描不同的 T 值，找到误差最小的那一个。这个思路不复杂，但代码里要处理的范围搜索和插值逻辑比较多，我会在下一节讲几个调试心得。

5. 常见问题速查与避坑记录

这部分是我觉得最有价值的内容。下面这些情形，都是实际执行中大概率会遇到的问题，我按“现象-原因-方案”整理成了表格。

5.1 关于曲线“看起来全为 0”的补充

我遇到过一位同学，把普朗克公式在 300 K 室温下画出来，纵轴数量级是 1e-10 比 1e-15，图上看起来确实是零。他自己一度以为是公式错了。其实公式没错，只是温度太低，辐射强度太弱，显示上完全被忽略。解决办法就是归一化，或者只看相对强度。

另外，如果温度很低、波长又很长，指数项exp(h*c/(λ*k*T))会非常接近 1，直接相减容易出现严重的灾难性抵消。这时可以用数值上更稳定的等价形式：B = (2hc^2 ./ λ.^5) · 1/(expm1(hc/(λk*T))) 吗？Matlab 里有expm1函数，专门计算 e^x - 1，在 x 很小时精度比exp(x)-1好得多。遇到低温长波场景时，直接换成expm1：

B = (2*h*c^2 ./ lam.^5) ./ expm1(h*c ./ (lam*k*T));

这个细节看起来不起眼，但对算 300 K、10 μm 这种场景，精度差距非常明显。

5.2 别忘了用维恩定律做“交叉验证”

我自己的习惯是，不管用哪种方法求峰值，都会先在代码里算一行：

lam_peak_theory = 2.8978e-3 / T;

然后跟数组计算结果做一次对比。如果误差小于 1%，代码基本没问题；如果误差大于 5%，赶紧回头查单位。这个习惯帮我省了很多排查时间，因为很多错误不是逻辑错误，而是量纲错误。量纲错误在普朗克公式里特别隐蔽，因为指数项可以因为单位不对而完全合理，但峰值偏到完全离谱的地方。

5.3 图形标注的小技巧

画完曲线后，我通常会在峰值点附近加一个文本注释，标注具体的波长值，这样图拿出去可以直接放在报告里，不用另外解释。

text(lam_peak*1e6 + 0.05, 0.95, ... ['\lambda_{max} = ', num2str(lam_peak*1e6, '%.3f'), ' \mum'], ... 'FontSize', 10);

再配合xline(lam_peak*1e6, '--', 'Peak')画一条垂直虚线，视觉上会很清晰。如果你对 Matlab 的图形对象熟，还能用 annotation 做更精细的箭头标注，但对大多数场景，text加xline已经足够了。

6. 从绘图到实际项目：我的一点个人体会

说回到项目本身。用 Matlab 绘制普朗克曲线并寻找峰值点，本质上是一个“模型-实现-验证”的经典流程。模型很简单，公式就一行；实现也不复杂，但真正容易出问题的恰恰是单位、数值稳定性和采样精度这类“细节中的魔鬼”。我之所以在文章里反复强调单位制和维恩定律验证，就是因为这些地方一旦出错，不是报错，而是给你一张“看起来完全正常但数值全错”的图。

我自己的经验是，做这类物理仿真时，永远不要相信一次跑通的结果。先画单条曲线，和理论值手算对比；再画多条曲线，看趋势是否符合维恩定律；最后再进入峰值提取和拟合。每增加一个步骤，就做一次交叉检查。这样从源头规避了大量低效的调试时间。

如果你想继续扩展这套代码，可以考虑把不同温度下的峰值波长自动提取出来，画一条“峰值波长-温度”曲线，看看它和 1/T 的反比关系是否完美贴合；也可以把自己采集的光谱数据读进来，经过插值、归一化后与不同温度的理论曲线做拟合，反推等效黑体温度。这些延伸场景都建立在今天这篇博文的基础能力上，后面的路会很顺。