docs: def reflectance_from_tbs(self, sun_zenith, tb_near_ir, tb_thermal, **kwargs):
          """
          The relfectance calculated is without units and should be between 0 and 1.
  
          Inputs:
  
            sun_zenith: Sun zenith angle for every pixel - in degrees
  
            tb_near_ir: The 3.7 (or 3.9 or equivalent) IR Tb's at every pixel
                        (Kelvin)
  
            tb_thermal: The 10.8 (or 11 or 12 or equivalent) IR Tb's at every
                        pixel (Kelvin)
  
            tb_ir_co2: The 13.4 micron channel (or similar - co2 absorption band)
                       brightness temperatures at every pixel. If None, no CO2
                       absorption correction will be applied.
  
          """
          # Check for dask arrays
          if hasattr(tb_near_ir, 'compute') or hasattr(tb_thermal, 'compute'):
              compute = False
          else:
              compute = True
          if hasattr(tb_near_ir, 'mask') or hasattr(tb_thermal, 'mask'):
              is_masked = True
          else:
              is_masked = False
  
          if np.isscalar(tb_near_ir):
              tb_nir = np.array([tb_near_ir, ])
          else:
              tb_nir = np.asanyarray(tb_near_ir)
  
          if np.isscalar(tb_thermal):
              tb_therm = np.array([tb_thermal, ])
          else:
              tb_therm = np.asanyarray(tb_thermal)
  
          if tb_therm.shape != tb_nir.shape:
              errmsg = 'Dimensions do not match! {0} and {1}'.format(
                  str(tb_therm.shape), str(tb_nir.shape))
              raise ValueError(errmsg)
  
          tb_ir_co2 = kwargs.get('tb_ir_co2')
          lut = kwargs.get('lut', self.lut)
  
          if tb_ir_co2 is None:
              co2corr = False
              tbco2 = None
          else:
              co2corr = True
              if np.isscalar(tb_ir_co2):
                  tbco2 = np.array([tb_ir_co2, ])
              else:
                  tbco2 = np.asanyarray(tb_ir_co2)
  
          if not self.rsr:
              raise NotImplementedError("Reflectance calculations without "
                                        "rsr not yet supported!")
  
          # Assume rsr is in microns!!!
          # FIXME!
          self._rad3x_t11 = self.tb2radiance(tb_therm, lut=lut)['radiance']
          thermal_emiss_one = self._rad3x_t11 * self.rsr_integral
  
          l_nir = self.tb2radiance(tb_nir, lut=lut)['radiance'] * self.rsr_integral
  
          if thermal_emiss_one.ravel().shape[0] < 10:
              LOG.info('thermal_emiss_one = %s', str(thermal_emiss_one))
          if l_nir.ravel().shape[0] < 10:
              LOG.info('l_nir = %s', str(l_nir))
  
          sunzmask = (sun_zenith < 0.0) | (sun_zenith > 88.0)
          sunz = where(sunzmask, 88.0, sun_zenith)
  
          mu0 = np.cos(np.deg2rad(sunz))
          # mu0 = np.where(np.less(mu0, 0.1), 0.1, mu0)
          self._rad3x = l_nir
          self._solar_radiance = self.solar_flux * mu0 / np.pi
  
          # CO2 correction to the 3.9 radiance, only if tbs of a co2 band around
          # 13.4 micron is provided:
          if co2corr:
              self.derive_rad39_corr(tb_therm, tbco2)
              LOG.info("CO2 correction applied...")
          else:
              self._rad3x_correction = 1.0
  
          nomin = l_nir - thermal_emiss_one * self._rad3x_correction
          denom = self._solar_radiance - thermal_emiss_one * self._rad3x_correction
          data = nomin / denom
          mask = (self._solar_radiance - thermal_emiss_one *
                  self._rad3x_correction) < EPSILON
  
          logical_or(sunzmask, mask, out=mask)
          logical_or(mask, np.isnan(tb_nir), out=mask)
  
          self._r3x = where(mask, np.nan, data)
  
          # Reflectances should be between 0 and 1, but values above 1 is
          # perfectly possible and okay! (Multiply by 100 to get reflectances
          # in percent)
          if hasattr(self._r3x, 'compute') and compute:
              res = self._r3x.compute()
          else:
              res = self._r3x
          if is_masked:
              res = np.ma.masked_array(res, mask=np.isnan(res))
          return res