docs: def get_reflectance(self, sun_zenith, sat_zenith, azidiff, bandname, redband=None):
          """Get the reflectance from the three sun-sat angles"""
          # Get wavelength in nm for band:
          if isinstance(bandname, float):
              LOG.warning('A wavelength is provided instead of band name - ' +
                          'disregard the relative spectral responses and assume ' +
                          'it is the effective wavelength: %f (micro meter)', bandname)
              wvl = bandname * 1000.0
          else:
              wvl = self.get_effective_wavelength(bandname)
              wvl = wvl * 1000.0
  
          rayl, wvl_coord, azid_coord, satz_sec_coord, sunz_sec_coord = self.get_reflectance_lut()
  
          # force dask arrays
          compute = False
          if HAVE_DASK and not isinstance(sun_zenith, Array):
              compute = True
              sun_zenith = from_array(sun_zenith, chunks=sun_zenith.shape)
              sat_zenith = from_array(sat_zenith, chunks=sat_zenith.shape)
              azidiff = from_array(azidiff, chunks=azidiff.shape)
              if redband is not None:
                  redband = from_array(redband, chunks=redband.shape)
  
          clip_angle = rad2deg(arccos(1. / sunz_sec_coord.max()))
          sun_zenith = clip(sun_zenith, 0, clip_angle)
          sunzsec = 1. / cos(deg2rad(sun_zenith))
          clip_angle = rad2deg(arccos(1. / satz_sec_coord.max()))
          sat_zenith = clip(sat_zenith, 0, clip_angle)
          satzsec = 1. / cos(deg2rad(sat_zenith))
          shape = sun_zenith.shape
  
          if not(wvl_coord.min() < wvl < wvl_coord.max()):
              LOG.warning(
                  "Effective wavelength for band %s outside 400-800 nm range!",
                  str(bandname))
              LOG.info(
                  "Set the rayleigh/aerosol reflectance contribution to zero!")
              if HAVE_DASK:
                  chunks = sun_zenith.chunks if redband is None else redband.chunks
                  res = zeros(shape, chunks=chunks)
                  return res.compute() if compute else res
              else:
                  return zeros(shape)
  
          idx = np.searchsorted(wvl_coord, wvl)
          wvl1 = wvl_coord[idx - 1]
          wvl2 = wvl_coord[idx]
  
          fac = (wvl2 - wvl) / (wvl2 - wvl1)
          raylwvl = fac * rayl[idx - 1, :, :, :] + (1 - fac) * rayl[idx, :, :, :]
          tic = time.time()
  
          smin = [sunz_sec_coord[0], azid_coord[0], satz_sec_coord[0]]
          smax = [sunz_sec_coord[-1], azid_coord[-1], satz_sec_coord[-1]]
          orders = [
              len(sunz_sec_coord), len(azid_coord), len(satz_sec_coord)]
          f_3d_grid = atleast_2d(raylwvl.ravel())
  
          if HAVE_DASK and isinstance(smin[0], Array):
              # compute all of these at the same time before passing to the interpolator
              # otherwise they are computed separately
              smin, smax, orders, f_3d_grid = da.compute(smin, smax, orders, f_3d_grid)
          minterp = MultilinearInterpolator(smin, smax, orders)
          minterp.set_values(f_3d_grid)
  
          if HAVE_DASK:
              ipn = map_blocks(self._do_interp, minterp, sunzsec, azidiff,
                               satzsec, dtype=raylwvl.dtype, chunks=azidiff.chunks)
          else:
              ipn = self._do_interp(minterp, sunzsec, azidiff, satzsec)
  
          LOG.debug("Time - Interpolation: {0:f}".format(time.time() - tic))
  
          ipn *= 100
          res = ipn
          if redband is not None:
              res = where(redband < 20., res,
                          (1 - (redband - 20) / 80) * res)
  
          res = clip(res, 0, 100)
          if compute:
              res = res.compute()
          return res