code: def density(self, R, Rs, rho0):
          
          return rho0/(R/Rs*(1+R/Rs)**2)
docs: three dimenstional NFW profile
  
          :param R: radius of interest
          :type R: float/numpy array
          :param Rs: scale radius
          :type Rs: float
          :param rho0: density normalization (characteristic density)
          :type rho0: float
          :return: rho(R) density