function blam = planck(temp, lambda)
% 
% Calculate the Planck function for the given temperature and wavelength.
% 
% Arguments:
% 
%       temp          (input)    temperature of blackbody (Kelvin)
%       
%       lambda        (input)    wavelength at which to evaluate (meters)
%       
%       blam         (output)    value of Planck function 
%                                         (Joules per sec per square meter
%                                          per unit wavelength)
%
% MWR 3/30/2001
%
% Improved the physical constants used in the calculations,
% so that calculated peak (and integrated power) will match more closely the
% book values.
% MWR 4/7/2003
% 
  % sanity checks
  if (temp <= 0) 
    error('temperature must be greater than zero');
  end
  if (lambda <= 0)
    error('wavelength must be greater than zero');
  end

  % some physical constants:
  %
  %   speed of light (meters/sec)
  c = 2.9979246e8;
  %   Planck's constant (Joules-sec)
  h = 6.626069e-34;
  %   Boltzmann constant (Joules/Kelvin)
  k = 1.3806503e-23;

  % check to see if we can use the Rayleigh-Jeans approximation
  tiny = 0.001;
  arg = (h*c)/(lambda*k*temp);
  if (arg < tiny) 
    % yes, can use the Rayleigh-Jeans approximation
    blam = (2*k*temp*c)/(lambda^4);

  else
    % nope, must use the full formula
    top = 2*h*c*c / (lambda^5);
    bot = exp((h*c)/(lambda*k*temp)) - 1;

    if (bot == 0)
      blam = 0;
    else
      blam = top/bot; 
    end

  end

