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 per steradian)
//
// 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

