docs: def sample(self, swave, binned=True, fluxunits='counts'):
          """Sample the observation at the given wavelength.
          Also see :ref:`pysynphot-command-sample`.
  
          Parameters
          ----------
          swave : float
              Wavelength to sample.
  
          binned : bool
              Sample binned dataset (no interpolation).
              Otherwise, native (perform interpolation).
  
          fluxunits : {'counts'}
              Only the unit of counts is supported for now.
  
          Returns
          -------
          ans : float
              Sampled flux in given unit.
  
          Raises
          ------
          NotImplementedError
              Flux unit is not supported or non-scalar wavelength is given.
  
          ValueError
              Given wavelength out of range.
  
          """
          if self._binflux is None:
            self.initbinflux()
  
          if fluxunits != 'counts':
              s = "Sorry, only counts are supported at this time"
              raise NotImplementedError(s)
          else:
              #Save current fluxunits, in case they're different
              saveunits = None
              if not units.ismatch('counts', self.fluxunits):
                  saveunits = self.fluxunits
                  self.convert('counts')
  
  
          if binned:
              #Then we don't interpolate, just return the appropriate values
              #from binflux
              if np.isscalar(swave):
                  #Find the bin in which it belongs.
                  #_bin_edge[i] is the low edge of the bin centered
                  #at binwave[i].
  
                  idx = np.where(swave >= self._bin_edges)[0]
                  #idx[-1] is the largest edge that is still smaller
                  #than swave
                  try:
                      ans = self.binflux[idx[-1]]
                  except IndexError:
                      s = 'Value out of range: wavelength %g not contained in range [%g, %g]'
                      s = s % (swave, self.binwave[0], self.binwave[-1])
                      raise ValueError(s)
  
              else:
                  #The logic for this case doesn't yet work on arrays
                  s = "Sorry, only scalar samples are supported at this time"
                  raise NotImplementedError(s)
  
          else:
              #Then we do interpolate on wave/flux
              if np.isscalar(swave):
                  delta = 0.00001
                  wv = np.array([swave - delta, swave, swave + delta])
                  ans = np.interp(wv, self.wave, self.flux)[1]
              else:
                  # This raises UnboundLocalError -- needs to be fixed!
                  ans = np.interp(wv, self.wave, self.flux)
  
  
          #Change units back, if necessary, then return
          if saveunits is not None:
              self.convert(saveunits)
          return ans