code: def density(self, R, Rs, rho0, r_core):
          
  
          M0 = 4*np.pi*rho0 * Rs ** 3
          return (M0/4/np.pi) * ((r_core + R)*(R + Rs)**2) ** -1
docs: three dimenstional truncated NFW profile
  
          :param R: radius of interest
          :type R: float/numpy array
          :param Rs: scale radius
          :type Rs: float
          :param rho0: density normalization (central core density)
          :type rho0: float
          :return: rho(R) density