% 518 Ps6_1 example.
% Brian Hornbuckle, February 20, 2008.

clear;

f = logspace(5,17,1000);         % frequency, Hz
[dummy,length_f] = size(f);      % find out length of f vector

k = 1.3807e-23;				        % Boltzmann's constant, J/K
h = 6.6261e-34;				        % Plank's constant, J*s
c = 1/sqrt(4*pi*1e-7*8.854e-12);	% speed of light in vacuum, m/s

T = [100 3000];                       % temperatures, K
B = zeros(max(length(T)),length_f);   % initialize array to hold spectral brightness, W m^-2 sr^-1 Hz^-1

for i = 1:max(length(T))
  B(i,:) = 2.*h.*f.^3./c.^2./(exp(h.*f./k./T(i)) - 1);
end;

figure(1)

loglog(f,B(2,:),'m--',f,B(1,:),'b-.');  axis([1e7,1e14,1e-21,1e-8]);
xlabel('frequency, Hz');  ylabel('spectral brightness, W m^{-2} sr^{-1} Hz^{-1}');
legend([num2str(T(2)),' K'],[num2str(T(1)),' K'],2);
title('Planck Law');  
text(5e8,1e-20,'microwave region');
set(gca,'XTick',[1e7,1e9,1e11,1e13]);
grid on;