code: def interp_bin(self, egy_bins, dtheta, scale_fn=None):
          
  
          npts = 4
          egy_bins = np.exp(utils.split_bin_edges(np.log(egy_bins), npts))
          egy = np.exp(utils.edge_to_center(np.log(egy_bins)))
          log_energies = np.log10(egy)
  
          vals = self.interp(egy[None, :], dtheta[:, None],
                             scale_fn=scale_fn)
          wts = np.exp(self._wts_fn((log_energies,)))
          wts = wts.reshape((1,) + wts.shape)
          vals = np.sum(
              (vals * wts).reshape((vals.shape[0], int(vals.shape[1] / npts), npts)), axis=2)
          vals /= np.sum(wts.reshape(wts.shape[0],
                                     int(wts.shape[1] / npts), npts), axis=2)
          return vals
docs: Evaluate the bin-averaged PSF model over the energy bins ``egy_bins``.
  
          Parameters
          ----------
          egy_bins : array_like
              Energy bin edges in MeV.
  
          dtheta : array_like
              Array of angular separations in degrees.
  
          scale_fn : callable        
              Function that evaluates the PSF scaling function.
              Argument is energy in MeV.