code: def mass_3d(self, R, Rs, rho0, r_core):
          
  
          Rs = float(Rs)
          b = r_core * Rs ** -1
          c = R * Rs ** -1
  
          M0 = 4*np.pi*Rs**3 * rho0
  
          return M0 * (1+b**2) ** -1 * (0.5*np.log(1+c**2) + b**2*np.log(c*b**-1 + 1) - b*np.arctan(c))
docs: :param R: projected distance
          :param Rs: scale radius
          :param rho0: central core density
          :param r_core: core radius