def compute_counts(self, skydir, fn, ebins=None):
        """Compute signal and background counts for a point source at
        position ``skydir`` with spectral parameterization ``fn``.

        Parameters
        ----------
        skydir : `~astropy.coordinates.SkyCoord`

        ebins : `~numpy.ndarray`

        Returns
        -------
        sig : `~numpy.ndarray`
            Signal counts array.  Dimensions are energy, angular
            separation, and event type.

        bkg : `~numpy.ndarray`
            Background counts array.  Dimensions are energy, angular
            separation, and event type.

        """

        if ebins is None:
            ebins = self.ebins
            ectr = self.ectr
        else:
            ectr = np.exp(utils.edge_to_center(np.log(ebins)))

        skydir_cel = skydir.transform_to('icrs')
        skydir_gal = skydir.transform_to('galactic')

        sig = []
        bkg = []
        bkg_fit = None
        if self._gdiff_fit is not None:
            bkg_fit = []


        for psf, exp in zip(self._psf, self._exp):

            coords0 = np.meshgrid(*[skydir_cel.ra.deg, ectr], indexing='ij')
            coords1 = np.meshgrid(*[skydir_cel.dec.deg, ectr], indexing='ij')

            # expv = exp.interpolate(skydir_cel.icrs.ra.deg,
            #                       skydir_cel.icrs.dec.deg,
            #                       ectr)

            expv = exp.interpolate(coords0[0], coords1[0], coords0[1])

            coords0 = np.meshgrid(*[skydir_gal.l.deg, ectr], indexing='ij')
            coords1 = np.meshgrid(*[skydir_gal.b.deg, ectr], indexing='ij')

            bkgv = self._gdiff.interpolate(np.ravel(coords0[0]),
                                           np.ravel(coords1[0]),
                                           np.ravel(coords0[1]))
            bkgv = bkgv.reshape(expv.shape)
            

            # bkgv = self._gdiff.interpolate(
            #    skydir_gal.l.deg, skydir_gal.b.deg, ectr)

            isov = np.exp(np.interp(np.log(ectr), np.log(self._iso[0]),
                                    np.log(self._iso[1])))
            bkgv += isov

            
            s0, b0 = irfs.compute_ps_counts(ebins, expv, psf, bkgv, fn,
                                            egy_dim=1,
                                            spatial_model=self.spatial_model,
                                            spatial_size=self.spatial_size)

            sig += [s0]
            bkg += [b0]

            if self._iso_fit is not None:
                isov_fit = np.exp(np.interp(np.log(ectr), np.log(self._iso_fit[0]),
                                            np.log(self._iso_fit[1])))
            else:
                isov_fit = isov

            if self._gdiff_fit is not None:
                bkgv_fit = self._gdiff_fit.interpolate(np.ravel(coords0[0]),
                                                       np.ravel(coords1[0]),
                                                       np.ravel(coords0[1]))
                bkgv_fit = bkgv_fit.reshape(expv.shape)
                bkgv_fit += isov_fit
                s0, b0 = irfs.compute_ps_counts(ebins, expv, psf,
                                                bkgv_fit, fn, egy_dim=1,
                                                spatial_model=self.spatial_model,
                                                spatial_size=self.spatial_size)
                bkg_fit += [b0]

        sig = np.concatenate([np.expand_dims(t, -1) for t in sig])
        bkg = np.concatenate([np.expand_dims(t, -1) for t in bkg])
        if self._gdiff_fit is not None:
            bkg_fit = np.concatenate([np.expand_dims(t, -1) for t in bkg_fit])

        return sig, bkg, bkg_fit