docs: 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]))