用MATLAB实现普朗克公式：黑体辐射计算与数值仿真全流程

简介一套面向物理模拟、红外仿真及科研计算的MATLAB实现资源围绕普朗克黑体辐射公式展开。包内共有5个m文件压缩包整体仅3KB涵盖普朗克公式主程序、辐射度计算、曲线拟合与映射函数等模块代码结构清晰便于直接修改复用。实际计算时需先定义普朗克常数、光速与玻尔兹曼常数再设定波长与温度范围通过嵌套循环逐点计算黑体辐射出射度并利用slice函数绘制温度—波长—辐射强度三维关系图能直观呈现不同温度、波长下的辐射分布规律计算得到的辐射强度矩阵亦可导出用于后续定量分析对红外系统设计、热辐射分析和相关科研实验有直接参考价值。目前已有4842人学习下载对于需要掌握黑体辐射模型、开展红外仿真或理解物理公式编程实现的工程师和科研人员这套资源可作为从理论到代码的快速桥梁有效提升相关课题的开发效率。 做光学仿真、辐射测温、光谱分析的工程师或者研究生大概率都绕不开普朗克公式。它把温度和光谱辐亮度之间的关系用一条曲线讲得明明白白黑体在不同温度下辐射能量怎么分配全看这一条公式。用MATLAB计算普朗克公式是我自己在做LED色温评估和红外探测器响应预估时最常用的基础操作把流程理顺之后半小时内就能从零写出一个能出图、能积分、能验证的完整脚本。这篇文章就从最基本的物理形式讲起把单位制、数组运算、数值积分这些坑一个个填平适合刚接触辐射度量学、需要快速上手MATLAB数值仿真的同学也适合想把黑体辐射这一块彻底弄懂后再去做更复杂光学仿真的人。1. 项目概述普朗克公式到底在算什么1.1 黑体辐射与普朗克公式的现实用途所谓黑体可以理解成一个理想化的物体它把照到身上的所有辐射全部吸收不反射也不透射然后纯粹按照自身的温度向外辐射能量。1900年普朗克用一个大胆的量子化假设把经典物理解释不了的黑体辐射光谱给拟合了出来这个公式后来成了量子力学的起点之一。它的作用很直白给定一个温度T可以计算出任意波长λ处的光谱辐亮度。放在工程里这意味着如果我想知道一个2200K的加热炉表面在3~5μm红外波段发出的能量有多大或者一个6000K光源在可见光波段能量占比是多少都可以直接用这个公式算不需要做实验。正因为如此普朗克公式在辐射测温、遥感反演、光学系统信噪比估算、LED与激光照明设计里头都是最底层的地基。1.2 为什么用MATLAB做这件事严格说普朗克公式本身是初等函数手算或者Excel也能算几个点但实际问题从来不会只算单个温度单条曲线。我要对比多个温度下的光谱分布要做波段积分要反解峰值波长还要把结果跟维恩位移定律、斯特藩—玻尔兹曼定律相互验证这时候MATLAB的优势就出来了。首先是向量化计算非常干净一个温度向量配上一个波长向量可以直接通过点除得到整个二维谱面。其次是数值积分工具成熟integral函数对这类带指数项的振荡不强但动态范围很大的被积函数处理得很稳定。第三是绘图方便光谱曲线的横纵坐标、对数坐标、归一化显示都很好控制这对快速判断计算结果是否合理非常关键。相比于Python,MATLAB在这类物理公式验证场景下的上手成本更低尤其适合本科生或者刚切换仿真工具的研究生。2. 物理模型与公式形式的选择2.1 波长形式与频率形式的区别普朗克公式有两种最常见的写法按波长表示的光谱辐亮度B_λ(λ,T)和按频率表示的光谱辐亮度B_ν(ν,T)。这两个形式物理上完全等价只是自变量不同。它们之间存在一个雅可比变换关系B_ν B_λ × |dλ/dν|因为λ c/ν所以dλ/dν -c/ν²取绝对值后公式形式就变了。这个坑特别隐蔽。如果我在MATLAB里先用波长形式算出结果再想换到频率坐标画图绝不能直接把λ替换成c/ν然后继续用原来的公式否则峰值位置、曲线面积、纵坐标数值全都会乱掉。实际写代码的时候我习惯明确声明变量到底是波长形式还是频率形式注释里写清楚单位这样过几个星期再回来看还能一眼判断自己当时算的是什么。对于绝大多数光学应用波长形式更直观因为探测器、滤光片一般都以波长来标定所以我也建议新手先从B_λ入手。波长形式的标准写法是B_λ(λ,T) (2hc²/λ⁵) × 1/(exp(hc/(λk_BT)) - 1)其中h是普朗克常数c是真空光速k_B是玻尔兹曼常数。这个式子的物理单位是W/(m²·m·sr)表示单位立体角、单位波长间隔内单位面积的辐射功率。初学者很容易忽略“sr”这个立体角单位导致理解上总差着一层其实它就对应一个方向锥角内的能量。2.2 量纲与单位制的统一做数值仿真的第一条铁律就是统一单位。我在辅导学生项目时见过太多结果算出来完全不对最后发现只是某个常数少打了一个数量级。普朗克公式里有h、c、k_B三个常数用国际单位制表示分别是参数符号数值普朗克常数h6.62607015×10⁻³⁴ J·s真空光速c2.99792458×10⁸ m/s玻尔兹曼常数k_B1.380649×10⁻²³ J/K公式里波长必须用米温度必须用开尔文这样算出来的B_λ单位是W/(m²·m·sr)。实际画图时如果横轴波长用μm微米很多人习惯直接把横坐标显示成微米但代入公式前必须先把微米换算成米否则整个指数项的量级就错了。还有个常被忽略的小细节纵轴单位从“每米”变成“每微米”时数值要差10⁶倍。因为1m的波长间隔等于10⁶个μm的波长间隔所以如果把B_λ的单位写成W/(m²·μm·sr)数值上等于B_λ(m单位)乘以10⁻⁶。这个换算在跟实验数据对比时尤其重要实验仪器给出的光谱辐亮度经常是W/(m²·μm·sr)不换算的话纵轴直接差六个数量级。3. MATLAB实现从零开始写普朗克公式脚本3.1 基础脚本与第一张光谱曲线我的建议是先做一个极简版本不搞函数封装不搞GUI把核心计算逻辑摊在脚本里方便逐步检查。下面这段代码就是计算一个5800K黑体在0.1μm到10μm波段的光谱辐亮度并绘制曲线%% 基本常数定义国际单位制 h 6.62607015e-34; % 普朗克常数 J·s c 2.99792458e8; % 光速 m/s kB 1.380649e-23; % 玻尔兹曼常数 J/K %% 计算参数 lambda_m (0.1:0.005:10) * 1e-6; % 波长范围 0.1~10 μm单位 m T 5800; % 温度 K %% 第一、第二辐射常数简化书写 c1 2 * h * c^2; % W·m^2 c2 h * c / kB; % m·K %% 普朗克公式注意点除 B_lambda_m c1 ./ (lambda_m.^5) ./ (exp(c2 ./ (lambda_m .* T)) - 1); %% 换算成以μm为单位的谱辐亮度 B_lambda_um B_lambda_m * 1e-6; % W/(m^2·μm·sr) %% 绘图 plot(lambda_m * 1e6, B_lambda_um, LineWidth, 1.5); xlabel(波长 / \mum); ylabel(光谱辐亮度 / (W·m^{-2}·\mum^{-1}·sr^{-1})); title(5800K 黑体辐射光谱); grid on;这里最关键的是点除符号./。因为lambda_m是一个行向量T是一个标量所以lambda_m.^5、lambda_m.*T、以及整个分数都必须是逐元素运算。如果少了一个点MATLAB会尝试把两个矩阵做矩阵乘法或者矩阵除法轻则报错重则算出一个完全没有物理意义的结果而且不报错——这种情况最坑。再解释一下为什么横轴起点选0.1μm而不是0。普朗克公式在λ趋于0时确实趋于0但数值上直接代入0会导致lambda_m.^5等于0c1/0 得到无穷大exp(c2/0)也是无穷大最终得到一个NaN。所以在数值计算中永远不要让波长向量包含0我习惯从0.1μm开始对绝大多数场景已经足够因为0.1μm以下的极紫外辐射在常规工程中占比很小。3.2 批量计算不同温度与维恩位移验证实际项目中很少只算一个温度比如我想看300K、500K、1000K、2000K、5800K这五条曲线在同一个坐标系下的变化趋势最简单的方式是写一个for循环或者直接构造温度向量。单个温度时可以用标量T但多个温度时更推荐循环里逐个计算然后hold on绘图因为这样每条曲线的峰值位置和幅度都容易控制。T_list [300, 500, 1000, 2000, 5800]; figure; hold on; colors lines(length(T_list)); for i 1:length(T_list) T_i T_list(i); B_i c1 ./ (lambda_m.^5) ./ (exp(c2 ./ (lambda_m .* T_i)) - 1) * 1e-6; plot(lambda_m * 1e6, B_i, Color, colors(i, :), LineWidth, 1.2); end hold off; xlabel(波长 / \mum); ylabel(光谱辐亮度 / (W·m^{-2}·\mum^{-1}·sr^{-1})); legend(cellstr(num2str(T_list, %d K))); grid on; xlim([0 20]);跑完之后你会发现一个很有趣的现象300K的曲线峰值在10μm附件肉眼几乎看不出它在短波方向有能量5800K的曲线峰值在0.5μm附近也就是可见光范围。这可以引出维恩位移定律的验证λ_peak · T 2.897771955×10⁻³ m·K。用数值方式验证非常简单[B_max, idx] max(B_lambda_m); lambda_peak lambda_m(idx); disp(lambda_peak * T);离散网格下求峰值误差取决于波长向量的间隔。我习惯把间隔设成0.001μm得到的结果和理论值能差到百分之几以内想要更高精度可以把idx附近的点拿来做二次插值也可以直接用更细的波长步长。这个方法作为第一道自检非常合适——如果算出的峰值波长乘温度偏差超过0.5%就要回头检查常数有没有打错或者单位是否统一。4. 进阶辐射出射度积分与谱段能量占比4.1 数值积分验证斯特藩—玻尔兹曼定律普朗克公式给出的是光谱辐亮度也就是某个方向上的能量分布。如果想知道黑体向半空间所有方向辐射的总能量需要对立体角积分对朗伯辐射体来说这个积分的结果是乘一个π。因此黑体的辐射出射度M可以写成M π ∫₀^∞ B_λ dλ理论上这个积分的结果就是斯特藩—玻尔兹曼定律M σT⁴其中σ ≈ 5.670374419×10⁻⁸ W/(m²·K⁴)。用MATLAB做数值积分时我可以直接用integral函数但要小心积分区间不能真的从0到无穷。sigma 2 * pi^5 * kB^4 / (15 * h^3 * c^2); % 理论斯特藩常数 M_integral pi * integral((lam) c1 ./ (lam.^5) ./ (exp(c2 ./ (lam .* T)) - 1), ... 1e-9, 5e-4, ArrayValued, true); M_theory sigma * T^4; disp([M_integral, M_theory]);这段代码里积分下限选1e-91nm上限选5e-40.5mm。为什么不能选0和Inf因为lambda0时公式会算出NaNlambda很大时exp参数接近0公式退化成c1/(λ⁵·c2/(λT))这个值确实会趋于0但直接积分到Inf会增加不必要的计算量而且在温度较高时指数项在最远处变得太平缓数值积分器可能会浪费大量步长。实际工程中我一般根据峰值波长来选区间下限取λ_peak/1000上限取λ_peak×1000基本能保证积分覆盖99.99%以上的能量。还有一个更隐蔽的坑如果被积函数无法做到严格光滑integral可能返回一个收敛警告。这时候先检查是不是指数项溢出了如果是可以把exp(c2/(λT))拆开处理或者把整个公式改成更稳定的数值形式。4.2 自定义波段的能量占比计算很多时候我不需要整个光谱的总能量而是想知道某个波段占了多少。比如判断一个光源的可见光占比、或者计算某个红外探测器工作窗口内能接收到多少辐射能量这就要做定积分lam1 0.38e-6; lam2 0.78e-6; % 可见光波段 380~780nm MyBlackbodyIntegral (a, b) integral((lam) ... c1 ./ (lam.^5) ./ (exp(c2 ./ (lam .* T)) - 1), a, b, ArrayValued, true); E_total pi * MyBlackbodyIntegral(1e-9, 5e-4); E_visible pi * MyBlackbodyIntegral(lam1, lam2); ratio E_visible / E_total * 100; disp([可见光占比: , num2str(ratio), %]);这个脚本用来分析不同色温光源非常直观。5800K左右的光源可见光占比大约能到40%上下而3000K的暖色温光源可见光占比就会明显下降大量能量跑到红外区。这类计算在LED照明设计和植物光源配比中很实用理解之后完全可以把这个函数封装成一个blackbody_fraction(T, lam1, lam2)函数后续调用非常方便。我建议大家在跑完积分后主动用斯特藩—玻尔兹曼定律做一次交叉验证如果同样的温度下积分得到的总出射度跟σT⁴对不上要么是单位出了问题要么是积分区间截断太狠。这一步虽然简单却是整个仿真闭环中最能树立信心的操作。5. 常见问题与排查技巧5.1 高频报错单位与数组运算的锅我把这几年遇到的高频错误整理成一张速查表每次脚本报错或者图像异常时先对照一遍现象可能原因解决方法图像是一条水平直线或全为零波长单位写成微米后没转成米把lambda向量乘以1e-6再代公式曲线出现NaN尖峰波长向量包含0或负值起始波长从0.1μm起步峰值波长明显不对用了频率公式却拿波长当横轴统一成波长形式B_λ纵轴数值小到10⁻²⁶单位还是W/(m²·m·sr)纵轴乘以1e-6换成μm单位指数溢出警告c2/(λT)过大检查λ是否被微米占位了或者T是否太小温度是向量时报矩阵维度错误用了/而没用./把所有除号改成点除积分结果比σT⁴大好几倍忘了乘π或者积分区间重复明确出射度要乘π其中单位错误是最难察觉的。我自己的经验是写完代码后第一件事不是看曲线形状而是手动取一个特殊点做笔算。比如T5800K、λ0.5μm时我大概能估出数值应该在10⁸量级W/(m²·μm·sr)如果算出来差特别多立刻能意识到单位折算出了问题而不是等图像出来后再瞎猜。5.2 结果异常排查从峰值到积分怎么确认没算错如果曲线形状正常但数值可疑我推荐一套三层校验法。第一层用维恩位移定律找峰值波长然后乘温度看是否接近2.898×10⁻³ m·K。第二层用斯特藩定律把全域积分结果除以T⁴看是否接近5.67×10⁻⁸。第三层是随机抽查单点比如查T300K时λ10μm处的光谱辐亮度跟文献或已有数据做对比。这三层校验独立性强任何一层出问题都能定位到具体环节。维恩位移不对大概率是公式形式或波长范围问题斯特藩积分不对优先看立体角积分时是否漏了π单点不对则要看常数代入是否有误。实际中我发现学生写错常数时往往曲线形状看起来依然正常因为指数项主导了形状常数只影响绝对幅度因此单一靠“看图像是否合理”来验收是不可靠的必须靠独立校验。还有一点经验温度过低时比如300K曲线峰值在10μm附近可见光区的数值几乎为零。如果此时用线性坐标绘图会被y轴起点压得很低。遇到这种情形可以用semilogy对数坐标绘制或者把兴趣波段截出来单独放大显示。这不算错误只是绘图策略要跟着物理量动态范围走。6. 从普朗克公式到更复杂的光学仿真写到这里其实最基础的MATLAB普朗克公式计算已经完整闭环了。但如果你的目标不只是验证公式而是要用它去推进一个更大的项目那还可以做很多扩展。最直接的扩展是把脚本改造成函数输入温度和波长范围自动返回光谱曲线和峰值信息然后塞进光学设计流程里作为光源模型。比如做成像系统仿真时已知目标温度可以用这个函数估算目标在探测器波段的辐射功率再配合光学系统的透过率和探测器响应度就能估算出信号电压。我那段时间做红外探测器指标分解时这套流程用了很多次。第二个方向是做实测光谱的等效色温拟合。实测一个光源的光谱数据后用lsqcurvefit去拟合普朗克公式中的T就能估计这个光源的等效黑体温度。拟合之前要注意把实测数据的单位统一成W/(m²·μm·sr)同时先把峰值波段附近的坏点去掉拟合结果才会稳定。另外很多光学仿真还关心光子数分布而非能量分布。普朗克公式的能量形式除以单个光子能量hc/λ就可以得到光子数谱分布这在量子效率计算、光通信噪声分析里会用到。改动虽然不大但物理概念要理清一不小心就会把能量和光子的峰值位置搞混。我个人实际使用中的体会是普朗克公式的MATLAB实现并不难难的是每一步都能明确自己在算什么、单位是什么、结果怎么自检。只要把单位制、点运算、积分区间这三个关键点抓牢这套脚本就可以从一次作业变成一件长期可用的工具后续往任何光学仿真方向延伸都会非常顺手。如果你做完这个基础版本后有具体的使用场景卡住了欢迎对照这里的方法逐项排查这套自检逻辑在别处同样适用。本文还有配套的精品资源点击获取