def sigma_nfw(self):
def _centered_sigma(self):<EOL><INDENT>bigF = np.zeros_like(self._x)<EOL>f = np.zeros_like(self._x)<EOL>numerator_arg = ((<NUM_LIT:1.> / self._x[self._x_small]) +<EOL>np.sqrt((<NUM_LIT:1.> / (self._x[self._x_small]**<NUM_LIT:2>)) - <NUM_LIT:1.>))<EOL>denominator = np.sqrt(<NUM_LIT:1.> - (self._x[self._x_small]**<NUM_LIT:2>))<EOL>bigF[self._x_small] = np.log(numerator_arg) / denominator<EOL>bigF[self._x_big] = (np.arccos(<NUM_LIT:1.> / self._x[self._x_big]) /<EOL>np.sqrt(self._x[self._x_big]**<NUM_LIT:2> - <NUM_LIT:1.>))<EOL>f = (<NUM_LIT:1.> - bigF) / (self._x**<NUM_LIT:2> - <NUM_LIT:1.>)<EOL>f[self._x_one] = <NUM_LIT:1.> / <NUM_LIT><EOL>if np.isnan(np.sum(f)) or np.isinf(np.sum(f)):<EOL><INDENT>print('<STR_LIT>')<EOL><DEDENT>if f.ndim == <NUM_LIT:2>:<EOL><INDENT>sigma = <NUM_LIT> * self._rs_dc_rcrit * f<EOL><DEDENT>else:<EOL><INDENT>rs_dc_rcrit_4D = self._rs_dc_rcrit.T.reshape(<NUM_LIT:1>, <NUM_LIT:1>,<EOL>f.shape[<NUM_LIT:2>],<EOL>f.shape[<NUM_LIT:3>])<EOL>sigma = <NUM_LIT> * rs_dc_rcrit_4D * f<EOL><DEDENT>return sigma<EOL><DEDENT>def _offset_sigma(self):<EOL><INDENT>numRoff = self._numRoff<EOL>numTh = self._numTh<EOL>numRbins = self._nbins<EOL>maxsig = self._sigmaoffset.value.max()<EOL>roff_1D = np.linspace(<NUM_LIT:0.>, <NUM_LIT> * maxsig, numRoff)<EOL>theta_1D = np.linspace(<NUM_LIT:0.>, <NUM_LIT> * np.pi, numTh)<EOL>rMpc_1D = self._rbins.value<EOL>theta = theta_1D.reshape(numTh, <NUM_LIT:1>, <NUM_LIT:1>)<EOL>roff = roff_1D.reshape(<NUM_LIT:1>, numRoff, <NUM_LIT:1>)<EOL>rMpc = rMpc_1D.reshape(<NUM_LIT:1>, <NUM_LIT:1>, numRbins)<EOL>r_eq13 = np.sqrt(rMpc ** <NUM_LIT:2> + roff ** <NUM_LIT:2> -<EOL><NUM_LIT> * rMpc * roff * np.cos(theta))<EOL>_set_dimensionless_radius(self, radii=r_eq13, integration=True)<EOL>sigma = _centered_sigma(self)<EOL>inner_integrand = sigma.value / (<NUM_LIT> * np.pi)<EOL>sigma_of_RgivenRoff = simps(inner_integrand, x=theta_1D, axis=<NUM_LIT:0>,<EOL>even='<STR_LIT>')<EOL>sig_off_3D = self._sigmaoffset.value.reshape(<NUM_LIT:1>, <NUM_LIT:1>, self._nlens)<EOL>roff_v2 = roff_1D.reshape(numRoff, <NUM_LIT:1>, <NUM_LIT:1>)<EOL>PofRoff = (roff_v2 / (sig_off_3D**<NUM_LIT:2>) *<EOL>np.exp(-<NUM_LIT:0.5> * (roff_v2 / sig_off_3D)**<NUM_LIT:2>))<EOL>dbl_integrand = sigma_of_RgivenRoff * PofRoff<EOL>sigma_smoothed = simps(dbl_integrand, x=roff_1D, axis=<NUM_LIT:0>,<EOL>even='<STR_LIT>')<EOL>_set_dimensionless_radius(self)<EOL>sigma_sm = np.array(sigma_smoothed.T) * units.solMass / units.pc**<NUM_LIT:2><EOL>return sigma_sm<EOL><DEDENT>if self._sigmaoffset is None:<EOL><INDENT>finalsigma = _centered_sigma(self)<EOL><DEDENT>elif np.abs(self._sigmaoffset).sum() == <NUM_LIT:0>:<EOL><INDENT>finalsigma = _centered_sigma(self)<EOL><DEDENT>else:<EOL><INDENT>finalsigma = _offset_sigma(self)<EOL>self._sigma_sm = finalsigma<EOL><DEDENT>return finalsigma<EOL>
Calculate NFW surface mass density profile.

        Generate the 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
            Surface mass density profiles (ndarray, in astropy.units of
            Msun/pc/pc). Each row corresponds to a single cluster halo.