1. 普朗克公式里的单位才是画对曲线的第一道坎前几天晚上一个做光学测量的朋友给我发来一张Matlab画出来的曲线图说自己的普朗克曲线怪怪的温度设成6000K峰值波长出现在4800 nm附近和维恩位移定律给的483 nm差了整整十倍。我扫了一眼代码就发现问题不在公式而在单位。这种“用Matlab绘制普朗克曲线并寻找峰值点”的场景听上去不过是几行plot的事可实际做起来坑几乎都藏在单位换算、坐标尺度和搜索区间这些不起眼的地方。这篇文章就把我平时处理这类需求的完整流程整理出来从公式拆解、Matlab代码骨架、峰值点定位方法到绘图排版细节和一个个实测踩坑记录。适合刚接触热辐射模拟的本科生、需要画黑体辐射谱的研究生以及任何想快速用Matlab出一张靠谱谱线图的人。1.1 波长形式还是频率形式先确认你画的是哪条谱普朗克公式常见的写法有两种一种按波长给出一种按频率给出。大多数人在讨论“黑体辐射曲线”时默认用的是波长形式因为实验上更容易按波长去测量横轴直接是纳米或者微米看起来直观。波长形式的光谱辐射出射度表达式是M(λ, T) 2πhc² / [λ⁵ * (exp(hc / (λkT)) - 1)]其中 h 是普朗克常数k 是玻尔兹曼常数c 是真空光速。这个式子的物理含义是温度为 T 的黑体在波长 λ 附近单位波长间隔内、单位表面积上辐射出去的功率。注意这里的“单位波长间隔”如果 λ 用米做单位那么 M 的单位就是 W·m⁻²·m⁻¹也常写成 W·m⁻³。频率形式则是M(ν, T) 2πhν³ / [c² * (exp(hν / (kT)) - 1)]频率形式和波长形式不能直接画在同一张图上比较数值因为坐标间隔的度量方式不一样。从波长形式换到频率形式需要做变量替换而且纵轴的量纲也会改变。我在实际中见过不少同学把两个公式的横轴都标成“频率”然后把曲线叠加结果形状对不上原因是波长形式的横轴频率是 c/λ但纵轴没有做 Jacobian 变换。所以动手之前先问自己一句我要表达的是“随波长分布的辐射”还是“随频率分布的辐射”绝大多数项目用到的是前者下面所有内容也都围绕波长形式展开。1.2 从米制公式到纳米坐标一个换算因子引发的误差Matlab 里写普朗克公式本身没有难度真正容易翻车的是单位。教科书上公式里的 λ 默认是米但你画图的时候大概率希望横轴显示纳米因为可见光、近红外这些波段用纳米更舒服。我建议的做法是计算时严格用米制画图坐标再单独转成纳米两个变量分开存放。很多人图省事直接让 λ 变量以纳米为单位代入公式然后指数项里的 hc/(λkT) 因为 λ 大了 1e9 倍指数变小整个曲线被拉宽峰值位置也会错得离谱。下面是我常用的代码骨架% 物理常数 h 6.62607015e-34; % 普朗克常数单位 J·s k 1.380649e-23; % 玻尔兹曼常数单位 J/K c 2.99792458e8; % 真空光速单位 m/s % 波长范围横轴用纳米 lam_nm 200:1:3000; % 200 nm 到 3000 nm lam lam_nm * 1e-9; % 换算成米用于公式计算 % 设定温度 T 6000; % 单位 K % 按米制计算谱辐射出射度单位 W·m^-2·m^-1即 W·m^-3 E_m 2 * pi * h * c^2 ./ (lam.^5 .* (exp(h * c ./ (lam * k * T)) - 1)); % 转成“每纳米”的谱辐射出射度单位 W·m^-2·nm^-1 E E_m * 1e-9; % 画图 plot(lam_nm, E, LineWidth, 1.5); xlabel(波长 \lambda / nm); ylabel(M_\lambda / W·m^{-2}·nm^{-1}); grid on;关键点在于最后那个* 1e-9。米制算出的是每一米波长间隔对应的辐射功率而横轴坐标是每纳米同样一段曲线在 1 纳米间隔内的功率自然只有 1 米间隔的十亿分之一所以要乘以 1e-9 才是纵轴上一一对应的值。不做这一步曲线形状和峰值位置完全不变但纵轴数值会大得离谱比如 6000K 的峰值附近会出现 4e13 这种天文数字。对于只想看曲线形状的人来说乘不乘无所谓但如果你要把数据放进论文或者和实验值对比单位不统一基本没法看。1.3 一次算完所有温度向量化写法实际项目里很少只画一条曲线更常见的是画一组温度来展示峰值移动规律。这时候用 for 循环逐个温度计算并不丢人Matlab 跑几百个点也不慢但代码会显得啰嗦。更推荐的做法是建一个温度数组再把温度也当成向量一次性算出所有曲线存成矩阵后面做任何后处理都方便。Ts [300 500 800 1200 2000 3000 4000 6000]; cmap parula(numel(Ts)); figure(Position, [100 100 700 500]); hold on; for i 1:numel(Ts) T Ts(i); E_m 2 * pi * h * c^2 ./ (lam.^5 .* (exp(h * c ./ (lam * k * T)) - 1)); E E_m * 1e-9; plot(lam_nm, E, Color, cmap(i, :), LineWidth, 1.2); end hold off; xlabel(波长 \lambda / nm); ylabel(M_\lambda / W·m^{-2}·nm^{-1});注意我这里依然在循环内部重新计算每条曲线而不是把所有温度一下塞进一个向量。原因是每条曲线要单独使用plot循环写法最简单直接。如果你想把所有 E 存成一个矩阵也可以先预分配E_all zeros(numel(lam_nm), numel(Ts));再逐个填充。用parula颜色映射的好处是色带从深蓝到浅黄温度高低对应颜色深浅视觉上很自然比默认的 lines 色系更适合表达连续变化的温度序列。2. 找峰值点的三种姿势穷举、fminbnd、维恩定律对拍曲线画出来之后下一个核心需求就是找峰值点。峰值点就是某温度下光谱辐射出射度取最大值对应的波长这个波长就是维恩位移定律里的那个“位移”结果。找峰值看似简单但不同精度需求要选不同方法我下面按“最简单、最精确、最物理”三条路线分别说。2.1 穷举法采样密了就够用如果只是想要一个视觉效果正确的峰值点穷举法永远是最稳的。曲线已经算出来了直接用max函数取最大值和对应下标就行T 6000; E_m 2 * pi * h * c^2 ./ (lam.^5 .* (exp(h * c ./ (lam * k * T)) - 1)); E E_m * 1e-9; [maxE, idx] max(E); peak_nm lam_nm(idx);这段代码没有什么技术含量但有一个前提必须满足你的波长采样步长要足够小。如果 lam_nm 只取了 200:100:3000 这种 100 nm 一档的稀疏网格峰值位置很可能偏掉几十甚至上百纳米。我一般取步长 1 nm 或者 0.5 nm对于 300K6000K 这个范围的普朗克曲线来说已经足够平滑峰值定位精度完全够用。当然穷举法的精度受限于采样间隔没办法直接给出比 1 nm 更细的峰值。但只要不是做精密拟合1nm 的误差对于绝大多数物理实验来说已经可以接受。这个方法的另一个好处是同时能拿到maxE就是峰值处的纵坐标值后面标注峰值点的时候直接用省得再查一遍。2.2 fminbnd数值优化让峰值精度再进一步如果你对峰值位置的要求到了“微米级别还不够还得精确到纳米小数点”的程度那就用fminbnd。Matlab 的fminbnd是单变量函数在给定区间内的无约束极小值求解器它基于黄金分割搜索稳定性很好。由于我们要求最大值所以要对函数取负号。% 定义一个负的普朗克函数输入 x 的单位是米 M_neg (x) -(2 * pi * h * c^2 ./ (x.^5 .* (exp(h * c ./ (x * k * T)) - 1))); % 搜索区间单位米 x_low 1e-7; x_high 3e-5; options optimset(TolX, 1e-12); [x_peak, ~] fminbnd(M_neg, x_low, x_high, options); peak_nm x_peak * 1e9;这里最容易翻车的是搜索区间。fminbnd是局部搜索如果区间内没有峰值它会收敛到区间边缘给你一个完全没意义的结果。比如你搜 6000K 的峰值峰值在 483 nm区间却写成 1e-7 到 3e-5 没问题峰值在里面但如果你用同一个固定区间去搜 300K 的峰值峰值在 9660 nm 左右还在区间内也没问题。真正的问题是当你要从 100K 到 6000K 自动批量处理时100K 峰值在 29 μm固定区间上限 3e-5 就把它切掉了。我的习惯是先用穷举法粗定位再把粗定位结果附近的一个窄区间喂给fminbnd。这样既避免区间不含峰值的问题也降低了搜索区间过宽带来的浮点误差累积。2.3 维恩位移定律解析结果用来给前两种方法“判卷”数值方法算出来的峰值到底对不对最好的验证方式是跟维恩位移定律的解析结果对拍。维恩位移定律说的是峰值波长和温度的乘积是一个常数λ_peak * T b其中 b ≈ 2.898e-3 m·K这个常数不是经验凑出来的而是对普朗克公式求导令导数为零后得到的方程的解。令 x hc / (λ·k·T)则峰值对应 x ≈ 4.965114因此 b hc / (k * 4.965114) ≈ 2.897771955e-3 m·K。我在代码里一般直接用wien_nm 2.897771955e-3 / T * 1e9;然后把wien_nm和上面两种方法算出来的peak_nm做个差值对比。如果差值在 1 nm 以内说明你的采样足够密、搜索区间设置也合理如果差值明显偏大优先去检查单位换算。顺带一提峰值处的谱辐射出射度也可以用这个 x 值算出来公式是 E_max (2πh c² / λ_peak⁵) * 1/(exp(x) - 1)。我经常在脚本里同时输出峰值波长和峰值高度后者用来做相对比较比如看不同温度下峰值高度之比是否接近温度五次方之比。3. 横轴尺度与坐标习惯为什么你的曲线总是“看不出峰值”很多人在画多温度普朗克曲线时会遇到一个尴尬情况300K 的曲线几乎贴着横轴6000K 的曲线高耸入云所有细节都堆在一起峰值点根本看不清楚。这不是公式写错而是横轴尺度和纵轴归一化的选择问题。3.1 线性坐标与半对数坐标各有用武之地普朗克曲线的波长范围从紫外一路延伸到远红外跨度往往超过两个数量级。如果横轴用线性坐标300K 的峰值在 9.7 μm 附近6000K 的峰值在 483 nm 附近两者在横轴上的位置差了几十倍画在一张图里低温曲线的峰会非常“扁”几乎看不见。这时候有两种处理方式。第一种是用半对数坐标横轴取对数semilogx(lam_nm, E)。半对数坐标可以让紫外、可见、红外在一个图里都展开温度曲线之间不会挤在一起。但代价是峰值位置的视觉直观性变差因为对数坐标下 483 nm 和 9660 nm 的距离被压缩了读者不太容易直接从图上读出峰值点在哪。我的经验是如果目的是“展示峰值随温度移动”用线性坐标更合适因为峰值波长和温度反比线性坐标下曲线形态最直观如果目的是“展示从紫外到红外的全波段形态”半对数坐标更好。两者没有绝对的优劣取决于你准备把图画给谁看。3.2 归一化对比与半高全宽的数值提取另一种让多温度曲线可以放在同一张图里比较的办法是归一化。把每条曲线除以自己的最大值纵轴变成 0 到 1这样所有曲线都在同一个高度范围内读者可以一眼看出形状变化温度越高峰越窄曲线越“瘦”温度越低峰越宽曲线越“胖”。Matlab 里归一化很简单E_norm E / max(E); plot(lam_nm, E_norm);归一化之后还可以顺手算半高全宽FWHM。做法是先找到峰值对应的索引然后分别向左向右找到第一个纵坐标小于 0.5 的位置再用线性插值精确求交点half maxE / 2; left_idx find(E(1:idx) half, 1, last); right_idx idx find(E(idx:end) half, 1, first); left_nm interp1(E(left_idx - 1:left_idx), lam_nm(left_idx - 1:left_idx), half); right_nm interp1(E(right_idx:right_idx 1), lam_nm(right_idx:right_idx 1), half); fwhm_nm right_nm - left_nm;这个数值提取在热辐射分析里比峰值点本身还常用比如判断发光谱线的单色性、估算滤光片的透过带宽等。我第一次在课题里需要 FWHM 时是用 find max 手动写的后来发现 interp1 可以把离散网格的误差降得很低就不再用“找最近点”的粗糙办法。3.3 顺带一提能量谱峰值和光子数谱峰值不是一回事如果你做的方向是光电探测、太阳电池、光催化这类和光子计数相关的领域还要注意一个容易被忽略的地方上面讨论的普朗克曲线是能量谱而光子数流密度谱的峰值位置会比能量谱峰值位置偏向长波方向。原因是光子数谱的表达式里多了 1/(hc/λ) 这个因子相当于在能量谱基础上乘一个与波长成正比的比例系数。解析上能量谱峰值满足 hc/(λkT) 4.965114而光子数谱峰值满足的方程解约在 2.821对应峰值波长更长。我在某个关于红外探测的项目里就因为这个差异踩过坑用能量谱的维恩位移定律去对实验结果发现峰值对不上后来意识到探测器响应的是光子数而不是能量。如果你不确定该用哪个谱先想清楚你的“检测器”到底在测量什么物理量。4. 把曲线画成可以直接放进报告/论文的样子代码能跑、峰值能找这只能算完成了一半。实际写报告、做组会展示或者投论文的时候图的排版、标注、字体都会影响别人对结果的判断。我见过很多人曲线画得完全正确但横轴标签里没有单位、图例密密麻麻、峰值点没有特别标注导致一张很好的图信息传达效果大打折扣。下面说说我习惯的排版套路。4.1 画布、字体、线宽一次配好Matlab 默认弹出来的图窗其实不太适合直接用于报告。我每次画正图前会先设置画布大小和字体这样后续所有子图、标注都基于同一套参数不会出现图例字体忽大忽小的情况。figure(Position, [100 100 700 500]); set(gca, FontName, Times New Roman, FontSize, 12, LineWidth, 1);画布宽度 700、高度 500 是一个比较稳妥的比例适合单张曲线图。如果后面要并排两张子图我会把宽度调到 1200。字体用 Times New Roman 是因为绝大多数期刊模板对正文和图的字体要求就是 Times 系图内文字和正文统一观感协调。坐标轴标签我建议用 LaTeX 解释器这样希腊字母和上下标都能正常显示xlabel(波长 \lambda / nm, Interpreter, latex); ylabel($M_\lambda$ / W m^{-2} nm^{-1}, Interpreter, latex);注意如果用了 latex 解释器中文字符会显示不出来。所以图内标签要么全英文要么直接把中文放到图注里。我的策略是图内一律英文标签中文含义在正文或图注中解释这样投稿时不用为字体问题再折腾一遍。4.2 峰值点标注与图例的排版思路找到峰值点后我会用散点加文字标注把它突出出来。散点用填充圆点颜色和曲线本身略有区别文字标注放在峰值点的右上方避免遮挡曲线。hold on; plot(peak_nm, maxE, o, MarkerSize, 8, MarkerFaceColor, r, MarkerEdgeColor, k); text(peak_nm 30, maxE * 0.9, sprintf(\\lambda_{peak} %.1f nm, peak_nm), FontSize, 10); hold off;这里有个小细节低温度曲线的峰很宽峰值点附近曲线变化平缓文字标注即使离峰值点有一段距离也不会造成误导但在高温曲线里峰很尖锐标注稍微偏移就可能盖住旁边的曲线。我习惯把偏移量设置为横轴跨度的 1%~2%比如横轴 200~3000 nm偏移 30 nm 就是一个比较安全的值。图例方面多温度曲线时不要把所有温度都放进图例观感极差。我会在循环里只对代表性的两条曲线加DisplayName其余曲线设置HandleVisibility为off让图例只显示关键信息。比如画 300、500、800、1200、2000、3000、4000、6000K 八条曲线时图例只写 300K 和 6000K中间温度靠颜色映射去体现正文里再说明颜色从蓝到黄对应温度递增。4.3 多温度对比的三种子图布局当需要在同一次汇报里同时展示线性坐标、半对数坐标和归一化曲线时我会用一行三个子图横轴分别是线性、对数、线性纵轴分别是原始谱、原始谱、归一化谱。这样一张图能讲清楚三件事谱线绝对形状、宽波段覆盖范围、归一化后的峰宽变化。Figure figure(Position, [50 50 1200 360]); subplot(1, 3, 1); plot(lam_nm, E1, LineWidth, 1.2); xlabel(波长 / nm); ylabel(M_\lambda); title(线性坐标); subplot(1, 3, 2); semilogx(lam_nm, E1, LineWidth, 1.2); xlabel(波长 / nm); ylabel(M_\lambda); set(gca, XTick, [300 500 1000 2000 3000]); title(半对数坐标); subplot(1, 3, 3); plot(lam_nm, E1 / max(E1), LineWidth, 1.2); xlabel(波长 / nm); ylabel(归一化强度); title(峰值归一化);三个子图共享同一组数据只是表征方式不同。做组会汇报时这种图比一张大图放八条曲线更容易讲清楚因为听众可以分别看到绝对量级、宽波段覆盖和形状变化三个维度。5. 实测踩坑记录四个我反复遇到的“怪问题”这部分是我个人最想写的。普朗克曲线本身不复杂但我在不同项目里反复遇到几乎同样的几个坑每次排查都花不少时间。整理在这里希望能帮别人省点力气。5.1 峰值算出来是理论值的10倍——单位混用开篇提到的那个朋友就是这个问题。他把lam_nm直接当成米带进了公式曲线形状看起来还挺像那么回事但峰值从 483 nm 偏到了 4830 nm。为什么刚好是十倍因为 6000K 的峰值接近 483 nm他错把这个数当成了 4830 nm 附近的峰维恩位移定律一算差十倍。排查方法很简单用维恩位移定律给一个理论值再和代码输出对比。如果你发现比例正好是 1e3 或者 1e-3第一反应就应该是单位换算。这个坑这么容易踩是因为 1 nm 1e-9 m 这个换算是人人都知道的但在长公式里随手一写就容易漏掉一个因子。我的解决思路是把波长换算单独放一行变量并在注释里写清楚“lam 为米制lam_nm 为坐标用”。5.2 峰值点跑到搜索区间的边界——fminbnd 的局部性有一次我需要批量处理从 200K 到 1000K 的峰值提取偷懒给fminbnd写了一个固定搜索区间结果 800K 和 1000K 的结果全部返回区间的上边界。当时第一反应是公式写得不对检查了半天才发现是搜索区间压根没包含这些温度下的峰值。fminbnd是局部搜索它不会自己“发现”区间外的极小值。区间设置必须保证峰值在区间内部而且区间宽度不能太大否则黄金分割搜索在迭代次数有限时也可能停在不够精确的位置。我现在统一的做法是先用max粗定位拿到粗峰值后把搜索区间缩小到粗峰值附近比如[peak_nm*1e-9*0.8, peak_nm*1e-9*1.2]再喂给fminbnd。这样既快又稳还不用为每个温度手工设区间。5.3 纵轴出现 e13——谱辐出度的单位怎么换算很多人第一次画出普朗克曲线纵轴数字大到吓人6000K 的峰值附近能到 10 的 13 次方。这个数字在数学上没有错因为按米制算M_λ 的单位是 W·m⁻²·m⁻¹6000K 的峰值辐出度确实在 10 的 13 次方这个量级。但如果你在论文里写“谱辐射出射度 4.4e13”审稿人会懵因为这不符合常见的光谱数据表达习惯。解决办法就是我前面反复强调的* 1e-9。把纵轴换算成 W·m⁻²·nm⁻¹6000K 的峰值就变成了 4.4e4 左右很常见的数量级。这里特别提醒纵轴单位一定要在标签里写清楚如果写的是 W·m⁻²·nm⁻¹数据就必须对应每纳米间隔的值如果你忘了乘 1e-9标签就要写 W·m⁻²·m⁻¹。单位标签和数值必须对得上这是原则问题。5.4 温度一多图就糊——图例与颜色的取舍第一次画多温度对比时我一股脑把 20 个温度都画上去用了默认的lines色系结果相邻温度颜色区分度很低图例占了半个图幅。后来改成parula颜色映射从色带里取等间隔 20 个颜色图例只保留 4 个代表温度最低、最高、中间两个其余温度用LegendOff属性隐藏。这个取舍的意义在于多温度曲线的信息核心是“趋势”而不是“每一个温度的具体值”。读者需要看到的是峰值随温度单调移动、峰宽随温度变窄而不是逐个去读某个温度对应哪条曲线。把不重要的曲线视觉上降级把代表曲线突出图的表达能力反而更强。最后再分享一点个人经验我做了不少次普朗克曲线的绘制最深的体会是这活儿难点不在代码而在对“你画出来的东西到底代表什么物理量”保持清醒。单位换算出错峰值点偏差十倍这类问题靠调试代码往往发现不了但和维恩位移定律一对拍就漏了馅。所以我现在写任何曲线脚本都会在开头放一段注释写清楚常数、单位制、横轴变量是米还是纳米、纵轴是 W·m⁻²·m⁻¹ 还是 W·m⁻²·nm⁻¹。半年后翻出来还能直接改温度参数就用不用重新推理一遍。如果后续想继续往深了做可以试试在普朗克曲线基础上计算某段波长范围的辐射能量占比或者把实测光谱和普朗克谱做最小二乘拟合来估计光源色温。这些都是在“画曲线、找峰值”的基础上延伸出来的实际需求工具还是 Matlab思路还是同一套单位检验和数值方法。希望这篇整理能让你少走一点我走过的弯路。
