import astropy.units as u
import numpy as np
import gala.potential as gp
from gala.units import galactic, solarsystem
import matplotlib.pyplot as plt

pot = gp.NFWPotential(m=1E11*u.Msun, r_s=20.*u.kpc, units=galactic)
pos = np.zeros((3,100)) * u.kpc
pos[0] = np.logspace(np.log10(20./100.), np.log10(20*100.), pos.shape[1]) * u.kpc
m_profile = pot.mass_enclosed(pos)

plt.figure()
plt.loglog(pos[0], m_profile, marker='') # doctest: +SKIP
plt.xlabel("$r$ [{}]".format(pos.unit.to_string(format='latex')))
plt.ylabel("$M(<r)$ [{}]".format(m_profile.unit.to_string(format='latex')))
plt.tight_layout()