function []=Planck_test()
emis = 1.0;
wavelengths=0.4:0.01:4;
temperatures = 800:200:2400;
nwave = length(wavelengths);
ntemp = length(temperatures);
eB = zeros(ntemp,nwave);
Tnamda = zeros(1,ntemp);
for i=1:ntemp
    T0=temperatures(i);
    for j=1:nwave
        lamda = wavelengths(j);
        L1=planck(lamda,emis,T0);
        eB(i,j)=L1;
    end
end
subplot(1,2,1);
plot(wavelengths,eB);
xlabel('波长 (μm)');
ylabel('热辐射 (Wm-2μm-1)');
% legend('800K to 2400K');
for i=1:ntemp
   idx = find(eB(i,:)==max(eB(i,:))); 
   Tnamda(i) = wavelengths(idx(1));
end
subplot(1,2,2);
constant = 2897.8+zeros(size(Tnamda));
plot(Tnamda,Tnamda.*temperatures,'o',Tnamda,constant,'-r');
ylim([2800 3000]);
xlabel('最大温度对应波长(μm)');
ylabel('温度 (K)');
legend('模拟的维恩位移常数','维恩位移常数');
end

function r=planck(lambda,emis,T)
%lambda is um
lambda = lambda*1e-6;

pi=3.14159265359;
c=3e8;
h=(6.6256e-34)
k=(1.38054e-23)

c1=2*pi*h*c*c;
c2=h*c/k;

L = emis*c1/lambda^5/(exp(c2/(lambda*T))-1)
r= (L*1e-6);
end