Calculate NFW differential surface mass density profile.

        Generate the differential surface mass density profiles of each cluster
        halo, assuming a spherical NFW model. Optionally includes the effect of
        cluster miscentering offsets, if the parent object was initialized
        with offsets.

        Returns
        ----------
        Quantity
            Differential surface mass density profiles (ndarray, in
            astropy.units of Msun/pc/pc). Each row corresponds to a single
            cluster halo.

def deltasigma_nfw(self):
        """Calculate NFW differential surface mass density profile.

        Generate the differential surface mass density profiles of each cluster
        halo, assuming a spherical NFW model. Optionally includes the effect of
        cluster miscentering offsets, if the parent object was initialized
        with offsets.

        Returns
        ----------
        Quantity
            Differential surface mass density profiles (ndarray, in
            astropy.units of Msun/pc/pc). Each row corresponds to a single
            cluster halo.
        """
        def _centered_dsigma(self):
            # calculate g

            firstpart = np.zeros_like(self._x)
            secondpart = np.zeros_like(self._x)
            g = np.zeros_like(self._x)

            small_1a = 4. / self._x[self._x_small]**2
            small_1b = 2. / (self._x[self._x_small]**2 - 1.)
            small_1c = np.sqrt(1. - self._x[self._x_small]**2)
            firstpart[self._x_small] = (small_1a + small_1b) / small_1c

            big_1a = 8. / (self._x[self._x_big]**2 *
                           np.sqrt(self._x[self._x_big]**2 - 1.))
            big_1b = 4. / ((self._x[self._x_big]**2 - 1.)**1.5)
            firstpart[self._x_big] = big_1a + big_1b

            small_2a = np.sqrt((1. - self._x[self._x_small]) /
                               (1. + self._x[self._x_small]))
            secondpart[self._x_small] = np.log((1. + small_2a) /
                                               (1. - small_2a))

            big_2a = self._x[self._x_big] - 1.
            big_2b = 1. + self._x[self._x_big]
            secondpart[self._x_big] = np.arctan(np.sqrt(big_2a / big_2b))

            both_3a = (4. / (self._x**2)) * np.log(self._x / 2.)
            both_3b = 2. / (self._x**2 - 1.)
            g = firstpart * secondpart + both_3a - both_3b

            g[self._x_one] = (10. / 3.) + 4. * np.log(0.5)

            if np.isnan(np.sum(g)) or np.isinf(np.sum(g)):
                print('\nERROR: g is not all real\n', g)

            # calculate & return centered profile
            deltasigma = self._rs_dc_rcrit * g

            return deltasigma

        def _offset_dsigma(self):
            original_rbins = self._rbins.value

            # if offset sigma was already calculated, use it!
            try:
                sigma_sm_rbins = self._sigma_sm
            except AttributeError:
                sigma_sm_rbins = self.sigma_nfw()

            innermost_sampling = 1.e-10  # stable for anything below 1e-5
            inner_prec = self._numRinner
            r_inner = np.linspace(innermost_sampling,
                                  original_rbins.min(),
                                  endpoint=False, num=inner_prec)
            outer_prec = self._factorRouter * self._nbins
            r_outer = np.linspace(original_rbins.min(),
                                  original_rbins.max(),
                                  endpoint=False, num=outer_prec + 1)[1:]
            r_ext_unordered = np.hstack([r_inner, r_outer, original_rbins])
            r_extended = np.sort(r_ext_unordered)

            # set temporary extended rbins, nbins, x, rs_dc_rcrit array
            self._rbins = r_extended * units.Mpc
            self._nbins = self._rbins.shape[0]
            _set_dimensionless_radius(self)  # uses _rbins, _nlens
            rs_dc_rcrit = self._rs * self._delta_c * self._rho_crit
            self._rs_dc_rcrit = rs_dc_rcrit.reshape(self._nlens,
                                                    1).repeat(self._nbins, 1)

            sigma_sm_extended = self.sigma_nfw()
            mean_inside_sigma_sm = np.zeros([self._nlens,
                                             original_rbins.shape[0]])

            for i, r in enumerate(original_rbins):
                index_of_rbin = np.where(r_extended == r)[0][0]
                x = r_extended[0:index_of_rbin + 1]
                y = sigma_sm_extended[:, 0:index_of_rbin + 1] * x

                integral = simps(y, x=x, axis=-1, even='first')

                # average of sigma_sm at r < rbin
                mean_inside_sigma_sm[:, i] = (2. / r**2) * integral

            mean_inside_sigma_sm = mean_inside_sigma_sm * (units.Msun /
                                                           units.pc**2)

            # reset original rbins, nbins, x
            self._rbins = original_rbins * units.Mpc
            self._nbins = self._rbins.shape[0]
            _set_dimensionless_radius(self)
            rs_dc_rcrit = self._rs * self._delta_c * self._rho_crit
            self._rs_dc_rcrit = rs_dc_rcrit.reshape(self._nlens,
                                                    1).repeat(self._nbins, 1)
            self._sigma_sm = sigma_sm_rbins  # reset to original sigma_sm

            dsigma_sm = mean_inside_sigma_sm - sigma_sm_rbins

            return dsigma_sm

        if self._sigmaoffset is None:
            finaldeltasigma = _centered_dsigma(self)
        elif np.abs(self._sigmaoffset).sum() == 0:
            finaldeltasigma = _centered_dsigma(self)
        else:
            finaldeltasigma = _offset_dsigma(self)

        return finaldeltasigma