%溫度不同的普朗克黑體單色輻射能力與波長的曲線
clear                                  %清除變量
k=1.38054e-23;                         %玻爾茲曼常數
h=6.626e-34;                           %普朗克常數
c=2.997925e8;                          %光速
sigma=5.6688e-008;                     %斯特潘常數
b=0.0029;                              %維恩常數
t=1400:100:2000;                       %熱力學溫度向量
n=length(t);                           %向量長度
lambda=[0:0.01:5]*1e-6;                %波長向量
lambda(1)=eps;                         %給零加一小量使分母不為零
[T,L]=meshgrid(t,lambda);              %波長和溫度矩陣
M=2*pi*h*c^2./(exp(h*c./(k*T.*L))-1)./L.^5;%單色輻射能力
figure                                 %創建圖形窗口
plot(lambda*1e6,M)                     %畫曲線
hl=legend([repmat('\itT\rm=',n,1),num2str(t'),repmat('K',n,1)]);%標記圖例
fs=16;                                 %字體大小
set(hl,'fontsize',fs)                  %設置圖例大小
grid on                                %加網格
title('普朗克黑體單色輻射能力與波長的關系','fontsize',fs)%標題
xlabel('波長\it\lambda\rm/\mum','fontsize',fs)%橫坐標
yl='單色輻射能力\itM\rm(\it\lambda\rm,\itT\rm)/(W\cdotm^-^3)';%縱坐標字符串
ylabel(yl,'fontsize',fs)               %縱坐標
txt=['\itb\rm=' num2str(b) 'm\cdotK']; %維恩常數文本
txt=[txt ',\it\sigma\rm=' num2str(sigma) 'W/(m^2\cdotK^4)'];%斯特潘常數文本
text(0,max(M(:))/10,txt,'fontsize',fs) %顯示常數
hold on                                %保持圖像
[mx,ix]=max(M);                        %找最大值和下標
stem(lambda(ix)*1e6,mx,'--','filled')  %畫直桿圖
text(lambda(ix)*1e6,mx,[num2str(lambda(ix)'*1e6)],'fontsize',fs)%顯示峰值波長
t=1300:2020;                           %較密的溫度向量
lambda=b./t;                           %波長向量
m=2*pi*h*c^2./(exp(h*c./(k*t.*lambda))-1)./lambda.^5;%單色輻射能力向量
plot(lambda*1e6,m)                     %畫峰值曲線