docs: def sample(self, wave, interp=True):
          """Sample the spectrum at given wavelength(s).
  
          This method has two behaviors:
  
          * When ``interp=True``, wavelength(s) must be provided
            as a Numpy array. Interpolation is done in internal units
            (Angstrom and ``photlam``).
          * When ``interp=False``, wavelength must be a scalar number.
            The flux that corresponds to the closest matching wavelength
            value is returned. This option should only be used for sampling
            binned data in `~pysynphot.observation.Observation`.
  
          Parameters
          ----------
          wave : number or array_like
              Wavelength(s) to sample, given in user unit.
  
          interp : bool
              Allow flux interpolation. Default is `True`.
  
          Returns
          -------
          ans : number or array_like
              Sampled flux in user unit.
  
          Raises
          ------
          NotImplementedError
              Non-scalar wavelength set provided when interpolation
              is not allowed.
  
          """
          if interp:
              # convert input wavelengths to Angstroms since the __call__ method
              # will be expecting that
              angwave = self.waveunits.ToAngstrom(wave)
  
              # First use the __call__ to get it in photlam
              flux = self(angwave)
  
              if hasattr(self, 'primary_area'):
                  area = self.primary_area
              else:
                  area = None
  
              # Then convert to the desired units
              ans = units.Photlam().Convert(angwave, flux,
                                            self.fluxunits.name, area=area)
  
          else:
              # Get the arrays in the proper units
              wave_array, flux_array = self.getArrays()
              if N.isscalar(wave):
                  # Find the correct index
                  diff = abs(wave-wave_array)
                  idx = diff.argmin()
  
                  ans = flux_array[idx]
  
              else:
                  raise NotImplementedError(
                      "Interp=False not yet supported for non-scalars")
  
          return ans