function L = blackbody(W,T)

%BLACKBODY Blackbody radiation
%   L = BLACKBODY(W,T) returns the radiance from a blackbody with the
%   temperature T (Kelvin) at the wavelengths given by the vector W
%   (micrometers) according to Planck's radiation law.
%
%   The unit of L is the somewhat obscure but practical microflick, i.e.,
%   microW/(sr * cm^2 * micrometers).

% W vavelength (micrometers)
% T temperature (Kelvin)

% h = 6.6256e-34; %Planck (Js)
% c = 2.9979e8;   %velocity of light (m/s)
% k = 1.3805e-23; %Boltzmann (J/K)
% lambda = W/1e6; %wavelength (m)
% L = (2*h*c^2) ./ ((lambda.^5)  .* (exp( (h*c) ./ (k*T*lambda) )-1) );

if T == 0
   L = zeros(size(W));
   return
end

A = 1.190938933177920e+010; % 2*h*c^2 * (1e6)^5 / 1e6 (microns) /1e4 (cm^2) *1e6 (microW)
B = 1.438818271640710e+004; % h*c/k * 1e6

L = A ./ ((repmat(W,1,numel(T)).^5)  .* (exp( B ./ (W*T) )-1) );