def planck(wave, temp, wavelength=True):
    """The Planck radiation or Blackbody radiation as a function of wavelength 
    or wavenumber. SI units.
    _planck(wave, temperature, wavelength=True)
    wave = Wavelength/wavenumber or a sequence of wavelengths/wavenumbers (m or m^-1)
    temp = Temperature (scalar) or a sequence of temperatures (K)


    Output: Wavelength space: The spectral radiance per meter (not micron!)
            Unit = W/m^2 sr^-1 m^-1

            Wavenumber space: The spectral radiance in Watts per square meter
            per steradian per m-1:
            Unit = W/m^2 sr^-1 (m^-1)^-1 = W/m sr^-1

            Converting from SI units to mW/m^2 sr^-1 (cm^-1)^-1:
            1.0 W/m^2 sr^-1 (m^-1)^-1 = 0.1 mW/m^2 sr^-1 (cm^-1)^-1

    """
    units = ['wavelengths', 'wavenumbers']
    if wavelength:
        LOG.debug("Using {0} when calculating the Blackbody radiance".format(
            units[(wavelength == True) - 1]))

    if np.isscalar(temp):
        temperature = np.array([temp, ], dtype='float64')
    else:
        temperature = np.array(temp, dtype='float64')

    shape = temperature.shape
    if np.isscalar(wave):
        wln = np.array([wave, ], dtype='float64')
    else:
        wln = np.array(wave, dtype='float64')

    if wavelength:
        const = 2 * H_PLANCK * C_SPEED ** 2
        nom = const / wln ** 5
        arg1 = H_PLANCK * C_SPEED / (K_BOLTZMANN * wln)
    else:
        nom = 2 * H_PLANCK * (C_SPEED ** 2) * (wln ** 3)
        arg1 = H_PLANCK * C_SPEED * wln / K_BOLTZMANN

    arg2 = np.where(np.greater(np.abs(temperature), EPSILON),
                    np.array(1. / temperature), -9).reshape(-1, 1)
    arg2 = np.ma.masked_array(arg2, mask=arg2 == -9)
    LOG.debug("Max and min - arg1: %s  %s", str(arg1.max()), str(arg1.min()))
    LOG.debug("Max and min - arg2: %s  %s", str(arg2.max()), str(arg2.min()))
    try:
        exp_arg = np.multiply(arg1.astype('float32'), arg2.astype('float32'))
    except MemoryError:
        LOG.warning(("Dimensions used in numpy.multiply probably reached "
                     "limit!\n"
                     "Make sure the Radiance<->Tb table has been created "
                     "and try running again"))
        raise

    LOG.debug("Max and min before exp: %s  %s", str(exp_arg.max()),
              str(exp_arg.min()))
    if exp_arg.min() < 0:
        LOG.warning("Something is fishy: \n" +
                    "\tDenominator might be zero or negative in radiance derivation:")
        dubious = np.where(exp_arg < 0)[0]
        LOG.warning(
            "Number of items having dubious values: " + str(dubious.shape[0]))

    denom = np.exp(exp_arg) - 1
    rad = nom / denom
    rad = np.where(rad.mask, np.nan, rad.data)
    radshape = rad.shape
    if wln.shape[0] == 1:
        if temperature.shape[0] == 1:
            return rad[0, 0]
        else:
            return rad[:, 0].reshape(shape)
    else:
        if temperature.shape[0] == 1:
            return rad[0, :]
        else:
            if len(shape) == 1:
                return np.reshape(rad, (shape[0], radshape[1]))
            else:
                return np.reshape(rad, (shape[0], shape[1], radshape[1]))