简介面向热辐射与量子物理方向的学习者一份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.1910e8c2 乘上 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 μm2000 个点平均间隔约 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 μm6000K 对应峰值约 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) 在 λ 极小时会超过双精度浮点上限。当指数内部数值超过约 709exp 返回 Inf随后 Inf/Inf 产生 NaNplot 遇到 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 的曲线紧贴 y0放大也看不到峰值。原因辐射亮度的最大值随 T 的五次方增长。6000K 峰值约 3e7300K 峰值约 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 μm2000 个点在 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 曲线误读成异常信号之后我每次画多曲线图都强制走一遍峰值校验流程。希望帮到你。本文还有配套的精品资源点击获取
