docs: def initbinflux(self):
          """Calculate binned flux and edges.
  
          Flux is computed by integrating the spectrum
          on the specified binned wavelength set, using
          information from the natural wavelength set.
  
          Native wave/flux arrays should be considered samples
          of a continuous function, but not their binned counterparts.
          Thus, it makes sense to interpolate ``(wave, flux)`` but not
          ``(binwave, binflux)``.
  
          .. note::
  
              Assumes that the wavelength values in the binned
              wavelength set are the *centers* of the bins.
  
              Uses ``pysynphot.pysynphot_utils.calcbinflux()`` C-extension,
              if available, for binned flux calculation.
  
          """
          endpoints = binning.calculate_bin_edges(self.binwave)
  
          # merge these endpoints in with the natural waveset
          spwave = spectrum.MergeWaveSets(self.wave, endpoints)
          spwave = spectrum.MergeWaveSets(spwave,self.binwave)
  
          # compute indices associated to each endpoint.
          indices = np.searchsorted(spwave, endpoints)
          self._indices = indices[:-1]
          self._indices_last = indices[1:]
  
          # prepare integration variables.
          flux = self(spwave)
          avflux = (flux[1:] + flux[:-1]) / 2.0
          self._deltaw = spwave[1:] - spwave[:-1]
  
          # sum over each bin.
          if utils_imported is True:
              self._binflux, self._intwave = \
                pysynphot_utils.calcbinflux(len(self.binwave),
                                            self._indices,
                                            self._indices_last,
                                            avflux,
                                            self._deltaw)
          else:
              #Note that, like all Python striding, the range over which
              #we integrate is [first:last).
              self._binflux = np.empty(shape=self.binwave.shape,dtype=np.float64)
              self._intwave = np.empty(shape=self.binwave.shape,dtype=np.float64)
              for i in range(len(self._indices)):
                  first = self._indices[i]
                  last = self._indices_last[i]
                  self._binflux[i]=(avflux[first:last]*self._deltaw[first:last]).sum()/self._deltaw[first:last].sum()
                  self._intwave[i]=self._deltaw[first:last].sum()
  
          #Save the endpoints for future use
          self._bin_edges = endpoints