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
    x = r_use[None, :] / R_s[:, None]
    prof = self._fx_projected(x)
    norm = 2 * R_s * self._norm(M_use, R_s, c_M)
    prof = prof[:, :] * norm[:, None]
    if np.ndim(r) == 0:
        prof = np.squeeze(prof, axis=-1)
    if np.ndim(M) == 0:
        prof = np.squeeze(prof, axis=0)
    return prof