简介：面向热辐射与量子物理方向的学习者，一份Matlab绘图资源以普朗克定律为核心，系统展示不同温度下黑体光谱辐射强度随波长变化的分布规律，解决利用数值计算可视化抽象物理公式的问题。通过对比多条温度曲线，可直观呈现黑体单色辐射能力随温度的移动特征，适合本科物理实验、光学课程以及研究生数值仿真入门时参考。压缩包共3个文件，包括两个.m脚本（分别承担公式实现与图形绘制）和一张结果图，整体仅227KB，轻量便于快速复用。绘图过程涵盖波长范围设置、多曲线绘制和不均匀刻度坐标处理等典型操作，可当作Matlab科学可视化的练习案例。已有8371人学习下载，借助该资源可快速复现普朗克公式，并迁移到其他热辐射光谱分析场景，便于日常教学演示与课后自主实践。
1. 普朗克定律与 Matlab 多曲线图：同一个坐标系里的 300K 与 6000K，差距在哪里
如果把 300K 和 6000K 的两条黑体辐射曲线放进同一个线性坐标系，300K 那一条会紧贴横轴，看上去像一条没有起伏的直线。这不是算错了，而是光谱辐射强度从低温到高温跨越了约 6 个数量级，线性刻度天然会淹没低温曲线。用 Matlab 画普朗克定律的多曲线图，真正的功夫不在调用 plot，而在单位换算、常数选取和刻度处理。这篇文章按「物理公式 → 函数实现 → 多曲线绘制 → 不均匀刻度」的顺序，把从地表温度到太阳表面温度的黑体曲线变成一张能直接放进汇报里的图，适合做热红外仿真、地表温度反演或探测器波段预研的开发者直接对标复现。
2. 黑体辐射的物理底料：普朗克公式的两种形式与常数单位
2.1 波长域还是频率域：两个公式为什么不能混用
普朗克公式实际上有两个版本，区别在自变量取波长还是频率。波长形式的写法是：
B_λ(T) = 2hc² / λ⁵ × 1 / (exp(hc / λkT) − 1)
单位是 W·m⁻²·sr⁻¹·m⁻¹，含义是单位波长间隔内、单位面积、单位立体角上的辐射功率。频率形式则是：
B_ν(T) = 2hν³ / c² × 1 / (exp(hν / kT) − 1)
单位是 W·m⁻²·sr⁻¹·Hz⁻¹。问题在于，B_λ 和 B_ν 不是同一个数的两个写法，它们之间隔着坐标变换 dν = −c/λ² dλ。同一个黑体，按波长画和按频率画，峰值位置、曲线形状都不一样。如果查到文献数据后拿错版本去对齐，会出现「峰值怎么偏了」的困惑。
工程场景里，红外探测器是按波长响应的，遥感载荷的波段定义也写在中心波长和带宽上，所以画「光谱辐射强度按波长分布」时，必须用波长形式的 B_λ(T)。这也是标题里那句「按波长分布」的落点。频率形式一般只出现在理论推导或微波段分析里，做黑体曲线可视化时不要混进来。
2.2 c1/c2 常数与单位：把微米和开尔文塞进公式的正确姿势
波长形式常被简写成工程惯用结构：
B_λ(T) = c1 / [λ⁵ × (exp(c2 / λT) − 1)]
这里的两个常数来自波尔兹曼常数、普朗克常数和光速的组合。用国际单位制（λ 单位取米）时：
c1 = 2hc² ≈ 1.1910e-16 W·m⁻²·sr⁻¹·m⁻¹
c2 = hc / k ≈ 1.4388e-2 m·K
直接拿这两个数去算，λ 得先用米表示。但画图时大家习惯用微米，0.4 μm 的可见光、9 μm 的热红外读起来才直观。于是要做一次单位换算：λ 从米换成微米，公式中的 λ⁵ 产生 1e30 的因子，c1 乘上 1e24 变成 1.1910e8，c2 乘上 1e6 变成 1.4388e4。这时公式变为：
B = 1.1910e8 ./ (λ_um.^5) ./ (exp(1.4388e4 ./ λ_um ./ T) − 1)
单位对应变成 W·m⁻²·sr⁻¹·μm⁻¹。这套常数组是工程里最常用的，很多红外仿真代码里都能看到 1.1910e8 和 1.4388e4 这两个影子。
还有一个容易踩的变体：另一类教材把辐射出射度 M_λ（对半球空间积分后的量）写成 c1 = 3.7418e8，也就是 2πhc² 的微米制版本。M_λ = π × B_λ，两者差一个 π 因子。如果拿 3.7418e8 去算辐射亮度，结果整体偏大约 3.14 倍，跟实验数据对不上的时候先检查这里。
3. 从公式到曲线的 Matlab 实现：先写对再画图
3.1 函数封装：一个 planck 函数处理单位换算
不要在主脚本里堆公式，先把普朗克计算封装成独立函数。这样温度扫描、波段对比、后续验证都能复用同一个底层计算。函数里的常数用 2019 年后的 CODATA 推荐值，Matlab 双精度下完全够用。
function B = planck(lambda_um, T) % planck 计算黑体光谱辐射亮度 % 输入: % lambda_um - 波长, 单位 μm, 可为标量或列向量 % T - 温度, 单位 K, 可为标量或行向量 % 输出: % B - 光谱辐射亮度, 单位 W/(m^2·sr·μm) % 当 lambda_um 为列向量、T 为行向量时, 返回 numel(lambda)×numel(T) 矩阵 h = 6.62607015e-34; % 普朗克常数, J·s c = 299792458; % 真空中光速, m/s k = 1.380649e-23; % 玻尔兹曼常数, J/K lambda_m = lambda_um * 1e-6; % 波长单位由 μm 转为 m exp_arg = h * c ./ (k .* lambda_m .* T); % 指数内部量, 无量纲 B_per_m = (2 * h * c^2 ./ lambda_m.^5) ... ./ (exp(exp_arg) - 1); % 按每米波长计算的辐射亮度 B = B_per_m * 1e-6; % 换算到每微米波长 end这段代码里最关键的是单位路径：输入进来的微米先用 1e-6 转成米，所有物理常数保持国际单位制，最后再把结果乘 1e-6 换回微米制。这样中间不出现 1.1910e8 这种先折算好的常数，少一层出错机会。lambda_um 设计成列向量、T 设计成行向量，是为了在下一步直接利用 Matlab R2016b 之后的隐式扩展，一次生成整张二维矩阵，省掉两层 for 循环。
3.2 主脚本：用矩阵与循环生成六温度曲线
温度序列的选择要有物理意义，而不是随手写六个数。300K 是典型地表温度，500K 到 800K 覆盖工业热源和火灾场景，1500K 接近高温熔体，3000K 到 6000K 跨越恒星表面温度区间。这样一组温度画出来，正好覆盖地表温度反演和太阳辐射对比两个典型用途。
% 普朗克多曲线绘制主脚本 clear; clc; close all; % 波长范围: 0.1 μm 到 50 μm, 覆盖紫外、可见光、热红外 lambda = linspace(0.1, 50, 2000)'; % 列向量, 2000 个采样点 % 温度序列: 地表温度到太阳表面温度 T = [300 500 800 1500 3000 6000]; % 调用 planck 函数生成二维矩阵, 每列对应一个温度 B = planck(lambda, T); % 自检: 打印 300K 在 10 μm 处的辐射亮度 idx10 = find(lambda >= 10, 1, 'first'); fprintf('300K, 10μm: %.4f W/(m^2·sr·μm)\n', B(idx10, 1));linspace 采样 2000 个点是因为在 0.1 到 50 μm 范围内，最低温 300K 的曲线峰宽约 2 μm，2000 个点平均间隔约 0.025 μm，足够分辨峰值且不会让曲线出现可见折角。更多点数会增加计算量但视觉上没差别。温度序列写成行向量，与 lambda 的列向量配合，planck 函数内部自动展开成 2000×6 的矩阵，每一列对应一个温度下的光谱分布。
3.3 运行与自检：峰值和量级是否符合物理预期
跑完上面的脚本，终端会打印 300K 在 10 μm 处的辐射亮度，正常值应该在 9 到 10 W/(m²·sr·μm) 这个区间。如果偏差过大，先回头查单位换算和常数。这个自检很有必要，因为画图阶段如果曲线形态不对，人眼往往看不出是常数错还是单位错。
另一个值得手动验证的量是峰值波长。维恩位移定律给出 λ_max = 2898 / T，单位是 μm·K。300K 对应峰值约 9.66 μm，6000K 对应峰值约 0.483 μm。等会儿画完图后，直接观察曲线峰的位置，或者用 max 函数提取索引，看是否落在这两个数附近。这一步是后面第六章做定量校验的基础，先在这里建立心理预期。
4. 多曲线图与不均匀刻度图：让不同温度曲线同时现形
4.1 多曲线图的绘制：句柄、图例与配色
有了 B 矩阵，画多曲线图最直接的方式是循环 plot，并把每条线的手柄存进数组。存句柄的意义在于，后续要单独调整某条线的线宽、颜色，或者生成图例时，手里有对象而不是重新找线。
% 多曲线图: 线性坐标版本 figure('Color', 'w', 'Position', [100 100 900 600]); hold on; % 使用 parula 色图取 6 个颜色, 保证不同温度可区分 cmap = parula(numel(T)); h = gobjects(1, numel(T)); % 预分配句柄数组 for i = 1:numel(T) h(i) = plot(lambda, B(:, i), 'LineWidth', 1.8, ... 'Color', cmap(i, :), ... 'DisplayName', sprintf('T = %d K', T(i))); end grid on; xlabel('波长 λ / μm'); ylabel('光谱辐射亮度 B_λ / W·m^{-2}·sr^{-1}·μm^{-1}'); legend(h, 'Location', 'northeast'); title('普朗克黑体辐射曲线（线性坐标）');这个版本跑完后会立刻暴露问题：300K 和 500K 的曲线贴着零线，800K 也只能看到很小的峰，1500K 开始明显，3000K 和 6000K 则把坐标轴纵轴范围拉到千万级别。linear 坐标下，低温和高温曲线无法同时呈现细节。这不是代码 bug，而是数据本身跨越了 6 个数量级，线性刻度的低量程区域被压缩到视觉不可辨。所以下一步必须换成不均匀刻度。
4.2 不均匀刻度：对数纵轴与坐标轴设置
所谓不均匀刻度，指的是放弃等距刻度，让每个数量级在视觉上占据相近长度。对数坐标是最典型的实现方案。在普朗克曲线场景里，纵轴取对数后，300K 的峰值和 6000K 的峰值都能以明显可见的高度差出现，且曲线在短波段的衰减趋势变成近似直线，更容易观察维恩位移和曲线形状差异。
% 不均匀刻度版本: 纵轴对数坐标 figure('Color', 'w', 'Position', [100 100 900 600]); for i = 1:numel(T) semilogy(lambda, B(:, i), 'LineWidth', 1.8); hold on; end grid on; xlabel('波长 λ / μm'); ylabel('光谱辐射亮度 B_λ / W·m^{-2}·sr^{-1}·μm^{-1}'); legend(arrayfun(@(t) sprintf('T = %d K', t), T, ... 'UniformOutput', false), 'Location', 'northeast'); title('普朗克黑体辐射曲线（纵轴对数坐标）'); % 等价写法: 也可以先 plot 再改坐标轴属性 % set(gca, 'YScale', 'log');semilogy 绘制后，每条曲线的峰值都被抬高到可读区域，300K 曲线不再贴底。此时曲线特征非常清楚：低温曲线在长波段的斜率平缓，高温曲线峰值向短波移动，且所有曲线在长波端逐渐收拢。这里有个细节值得记住：普朗克公式在长波极限下退化为瑞利-金斯近似，B_λ 和 λ⁴ 成反比，所以对数坐标下所有曲线的长波段变成斜率相近的直线段。如果看到某条曲线在长波端被截断或者突然掉零，说明计算里有溢出或采样不足。
除了整体对数坐标，也可以手工设置不均匀的横轴刻度，只标出关键谱段。比如显示 8 到 14 μm 热红外窗口和 0.4 到 0.7 μm 可见光窗口，中间的大段红外只给少数刻度，让横轴的密度跟随业务关注度变化。
% 手工不均匀横轴刻度: 关注热红外窗口与可见光窗口 set(gca, 'XTick', [0.1 0.4 0.7 3 5 8 10 12 15 20 30 50]); set(gca, 'XTickLabel', {'0.1','0.4','0.7','3','5','8','10','12','15','20','30','50'});这种做法的价值在演示场景：如果观众关心的是 300K 地表辐射和 6000K 太阳辐射的重叠关系，横轴把 0.4、0.7、8、10、12 这几个点标出来，视线就会聚焦在这两个区间，而不是在整个 0.1 到 50 μm 上平均分布注意力。
4.3 局部放大与子图布局：同时看可见光细节与热红外全貌
对数纵轴解决的是数量级跨度，但横轴跨度同样大。0.1 μm 到 50 μm 横跨近三个数量级，线性横轴会把可见光部分挤在最左侧一小段。如果业务里要同时讨论太阳辐射的可见光峰和地表的红外峰，一种实用方案是在整图角落嵌入一个局部放大坐标轴。
% 在主图基础上叠加一个局部放大坐标轴 figure('Color', 'w', 'Position', [100 100 900 600]); for i = 1:numel(T) semilogy(lambda, B(:, i), 'LineWidth', 1.8); hold on; end grid on; xlabel('波长 λ / μm'); ylabel('光谱辐射亮度 B_λ / W·m^{-2}·sr^{-1}·μm^{-1}'); % 在右上角加入局部放大坐标轴, 显示可见光区段 ax = axes('Position', [0.55 0.55 0.32 0.3]); for i = 1:numel(T) plot(lambda, B(:, i), 'LineWidth', 1.5); hold on; end xlim([0.35 0.75]); % 可见光区段 ylim([0 6e7]); % 纵轴范围以 6000K 曲线为参考 xlabel('λ / μm'); ylabel('B_λ'); title('可见光区放大');局部放大图里能看到 300K 曲线在可见光段的数值几乎为零，尽管纵轴已经缩到 6000K 峰高附近，低温曲线仍然不可见。这个视觉结论有实际工程意义：地表自身辐射在可见光波段完全可以忽略，白天看到的可见光信号几乎全部来自太阳辐射反射，所以在做多光谱遥感时，可见光通道和热红外通道的物理机制必须分开建模。这一点对地表温度反演的波段选择很关键。
5. 避坑与排查：黑体曲线绘制最常见的五类翻车现场
5.1 短波端曲线突然消失或出现 NaN
现象：波长小于 0.01 μm 时，某些温度曲线在某一点后突然归零，或者出现断点，log 坐标下表现为一段空白。
原因：指数项 exp(hc / λkT) 在 λ 极小时会超过双精度浮点上限。当指数内部数值超过约 709，exp 返回 Inf，随后 Inf/Inf 产生 NaN，plot 遇到 NaN 直接断开。
解决：在 planck 函数里对指数参数做截断，大指数部分直接令辐射亮度为零。修改后的写法：
exp_arg = h * c ./ (k .* lambda_m .* T); exp_arg(exp_arg > 700) = Inf; B_per_m = (2 * h * c^2 ./ lambda_m.^5) ./ expm1(exp_arg); % 注意: exp_arg 为 Inf 时 expm1(Inf) = Inf, 相除后结果为 0这里把 exp(x) − 1 替换成 expm1，是因为 x 较小时直接减一会损失有效位，具体见 5.5。
5.2 峰值位置整体偏移几个数量级
现象：300K 曲线的峰值出现在 9.66e-6 μm 而不是 9.66 μm，或者 6000K 的峰值落在 483 nm 却画在 0.000483 μm 上。
原因：λ 输入单位是微米，但常数 c1、c2 仍然按米制代入，导致指数项和 λ⁵ 项里差了 1e6 的幂次。这是单位不统一造成的系统性偏移，各条曲线形状不变，但整体位置错位。
解决：统一采用函数内单位换算路径。所有输入以微米为准进入函数，函数内部先转米制计算，再转回微米制输出。不要在调用侧反复乘 1e6 或除 1e6，这种临时换算最容易漏一处。建议每次改完代码后，用 4.1 里那个 300K 10μm 的自检打印语句验证数值再继续。
5.3 线性坐标下低温曲线完全贴底
现象：同一张图里只能看到 3000K 和 6000K 两条曲线，300K 和 500K 的曲线紧贴 y=0，放大也看不到峰值。
原因：辐射亮度的最大值随 T 的五次方增长。6000K 峰值约 3e7，300K 峰值约 10，两者差距接近 7 个数量级，线性坐标下低温曲线被压缩到坐标轴底部。
解决：换用 semilogy 纵轴对数坐标。如果目的是展示各温度曲线的峰值位置而不是绝对量级，可以对每条曲线做峰值归一化再画线性图。归一化后的曲线能看清相对形状，但会牺牲绝对数值的比较能力。做地表温度反演演示时，建议保留对数版本作为主图，归一化版本作为辅助分析。
5.4 循环 plot 后图例失效或无法单独改线型
现象：循环里调用 plot 但不接收返回值，画完图后执行 legend，发现图例只能显示最后一条曲线，或者想把某条曲线加粗却找不到句柄。
原因：循环内如果没用 hold on，新 plot 会覆盖旧图；没存句柄时，后续操作只能靠 findobj 或重新绘图，非常被动。
解决：循环前执行 hold on，循环内用 h(i) = plot(...) 保存句柄，最后 legend(h, ...) 显式传句柄数组。要改某一温度的线型或颜色时，直接修改 h(i) 的 LineWidth 或 Color 属性，不需要重画整幅图。这也是把绘图代码当成可维护脚本而不是一次性命令来写的基本习惯。
5.5 exp(x)−1 在长波端出现数值抖动
现象：波长很大或温度很高时，长波端曲线出现轻微波动，数值上表现为尾部的有效位丢失，曲线不平滑。
原因：当 x = hc/λkT 较小（小于 1e-3）时，exp(x) 接近 1，直接执行 exp(x) − 1 会消除大量有效数字。双精度下 exp(1e-4) 约 1.000100005，减 1 后只保留约 6 位有效数字，而乘以极大的 λ⁵ 倒数后误差被放大。
解决：使用 expm1 替代 exp(x) − 1。expm1 是 Matlab 内建函数，专门计算 exp(x) − 1 且在小 x 时保持完整精度，R2016b 及以上版本直接可用。配合 5.1 的截断逻辑，这两个问题可以一并解决。
6. 进阶验证：维恩位移定律校验峰值与归一化叠置分析
6.1 峰值位置的自动提取与理论对比
画完多曲线图后，不要只停留在视觉确认。用代码提取每条曲线的峰值位置，跟维恩位移定律给出的理论值对比，能快速发现单位或常数错误。实现方式是用 max 函数沿波长方向找最大值索引，再根据索引取波长坐标。
% 提取峰值位置并与维恩位移定律对比 [maxB, idx] = max(B, [], 1); lambda_max = lambda(idx); lambda_theory = 2898 ./ T; % 打印对比结果 for i = 1:numel(T) fprintf('T = %4d K, 数值峰值 = %7.3f μm, 理论峰值 = %7.3f μm\n', ... T(i), lambda_max(i), lambda_theory(i)); end如果 lambda 在峰值附近采样足够密，数值峰值和理论值的偏差应当小于 1%。300K 的理论峰值 9.66 μm，2000 个点在 0.1 到 50 μm 范围内平均间隔约 0.025 μm，折合相对误差约 0.25%，完全在可接受范围内。如果偏差超过 2%，优先检查 planck 函数里的单位换算和 c2 常数。
6.2 归一化叠置：验证曲线形状的标度无关性
把所有温度曲线除以各自的峰值后叠在一张图上，会出现一个有价值的现象：曲线几乎完全重合。普朗克定律具有标度形式，B_λ / B_max 只是 λ/λ_max 的函数，温度只决定峰值位置和峰值大小，不改变相对曲线形状。利用这点，可以做一个通用曲线用于快速估算任意温度任意波长处的辐射亮度。
实际项目里，我曾经用归一化曲线做过红外探测器波段匹配的快速估算：先确认待测目标温度对应的峰值波长，再把探测器响应窗口换算成 λ/λ_max 区间，直接用通用曲线读响应比例，速度快且精度足够方案预研。这种做法比每改一个温度就跑一遍完整仿真省事得多，但前提是曲线经过了 6.1 的峰值校验，否则整条归一化曲线都会失真。
6.3 导出图形时的字体与格式收尾
图形要放进论文或汇报 PPT，需要在导出前检查几个视觉细节：字体设为 Times New Roman，字号不低于 10pt，曲线线宽保持 1.5 以上，图例不要遮挡曲线峰值。Matlab 里用 exportgraphics 导出矢量格式，效果比截图好得多。
% 导出前调整图形属性 set(gca, 'FontName', 'Times New Roman', 'FontSize', 11); exportgraphics(gcf, 'planck_curves.pdf', 'ContentType', 'vector');如果导出目标是演示文稿，直接存 600 dpi 的 PNG 更稳妥。矢量 PDF 插进 Word 有时会遇到字体替换问题，这时导出 PNG 反而少折腾。我的习惯是矢量 PDF 存底稿，600 dpi PNG 用于成品交付，两种格式都保留一份。后来我给某红外载荷项目做波段预研，第一件事就是跑这个脚本，先验证 λ_max 与 2898/T 的偏差再谈数据。自从那次把 300K 曲线误读成异常信号之后，我每次画多曲线图都强制走一遍峰值校验流程。希望帮到你。
本文还有配套的精品资源，点击获取
