def METHOD_NAME(self, cosmo, r, M, a):
    r_use = np.atleast_1d(r)
    M_use = np.atleast_1d(M)
    # Comoving virial radius
    R_M = self.mass_def.get_radius(cosmo, M_use, a) / a
    c_M = self.concentration(cosmo, M_use, a)
    R_s = R_M / c_M
    norm = self._norm(M_use, R_s, c_M)
    x = r_use[None, :] / R_s[:, None]
    prof = norm[:, None] / (x * (1 + x)**3)
    if self.truncated:
        prof[r_use[None, :] > R_M[:, None]] = 0
    if np.ndim(r) == 0:
        prof = np.squeeze(prof, axis=-1)
    if np.ndim(M) == 0:
        prof = np.squeeze(prof, axis=0)
    return prof